A Hybrid Finite Element—Pseudospectral Approach for Simulating Nonlinear Pollutant Transport in Watercourses

1. Introduction

Water pollution is one of the most critical environmental challenges of the 21st century, particularly in natural watercourses that serve as essential reservoirs for biodiversity, drinking water, agriculture, and industry. The propagation of pollutants in rivers, lakes, and estuaries involves complex physical, chemical, and biological processes governed by fluid flow, dispersion, and reaction mechanisms [2,6,14]. Understanding and predicting these transport phenomena is crucial for environmental risk assessment, water quality management, and the design of remediation strategies [16].

Over the past decades, mathematical modeling and numerical simulation have become indispensable tools in environmental hydrodynamics. They offer a cost-effective alternative to experimental monitoring, which is often limited by spatial coverage, temporal resolution, or financial constraints [3,8]. Various numerical approaches have been proposed for simulating pollutant transport, including finite difference methods [13], finite volume methods [10], and finite element methods (FEM) [4]. While FEM provides flexibility in handling complex geometries, spectral and pseudospectral methods are known for their superior accuracy in resolving smooth solutions [1,15].

Hybrid numerical schemes that combine the geometric adaptability of FEM with the precision of spectral methods have gained attention in recent years, especially for advection-diffusion-reaction equations [9,12]. These methods can yield high-order accuracy while maintaining stability under stiff nonlinearities or irregular domains. However, their application to environmental flow problems-particularly those involving coupled Navier-Stokes and reaction-diffusion systems-remains under explored.

In this work, we propose a hybrid numerical approach that couples the finite element method with a pseudospectral discretization scheme to simulate the nonlinear transport of pollutants in a two-dimensional watercourse. The pollutant concentration is governed by a convection-diffusion-reaction equation, while the flow velocity is computed from the steady-state incompressible Navier-Stokes equations. The nonlinear reaction term models degradation or transformation processes that may occur during pollutant dispersion.

Our main contributions are as follows:

• We formulate a coupled nonlinear system for pollutant transport in incompressible flows, incorporating both diffusion and nonlinear reaction dynamics.
• We develop a hybrid FEM-pseudospectral scheme that ensures stability and convergence even under limited solution regularity.
• We conduct a set of numerical experiments to illustrate the effects of physical parameters (e.g., Reynolds number, reaction nonlinearity) on pollutant distribution.

The remainder of the paper is organized as follows: Section 2 presents the mathematical model and governing equations. Section 3 details the spatial discretization and hybrid numerical scheme. Section 4 describes the algorithmic implementation and simulations. Section 5 analyzes the numerical results. Finally, Section 6 concludes the paper with perspectives for future work.

We consider the transport and diffusion of a chemical pollutant in a two-dimensional, bounded aquatic domain $\mathsf{\Omega }\subset {\mathbb{R}}^2$, over a finite time interval $t\in [0,T]$. The pollutant concentration, denoted by $C(t,x,y)$, is governed by a convection?diffusion-reaction equation:

$${\displaystyle \frac{\partial C}{\partial t}}-\lambda \mathrm{\Delta }C+\mathbf{v}·\nabla C+F\left(C\right)=0,\mathrm{in}(0,T]\times \mathsf{\Omega },$$

(1)

subject to homogeneous Neumann boundary conditions:

$${\displaystyle \frac{\partial C}{\partial \mathbf{n}}}=0,\mathrm{on}(0,T]\times \partial \mathsf{\Omega },$$

(2)

and an initial condition:

$$C(0,x,y)=C_0(x,y),\mathrm{in}\mathsf{\Omega }.$$

(3)

Here, $\mathbf{v}=(v_1,v_2)$ is the steady, incompressible velocity field of the ambient fluid, and $F\left(C\right)$ represents a nonlinear reaction term describing chemical transformation, biological decay, or pollutant degradation.

The nonlinear function $F\left(C\right)$ is defined as:

$$F\left(C\right)={\mu C\left|C\right|}^{\beta -1},$$

(4)

where $\mu \ge 0$ is a reaction rate constant and $\beta >0$ is a nonlinearity parameter.

This form is widely used to represent autocatalytic reactions, nonlinear degradation, or sorption-desorption phenomena in environmental models [11]. For example:

• When $\beta =1$, the reaction is linear and corresponds to classical first-order kinetics.
• For $\beta >1$, the nonlinearity introduces reaction acceleration with increasing concentration, modeling stronger pollutant self-interaction or saturation effects.

The use of $F\left(C\right)={\mu C\left|C\right|}^{\beta -1}$ also ensures continuous differentiability and generality across a wide class of power-law reaction models.

The velocity field $\mathbf{v}$ is assumed to satisfy the steady incompressible Navier-Stokes equations:

$$\begin{array}{cc}\hfill -\mathrm{\Delta }{v}_1+F_1(v_1,v_2)& =f_1,\mathrm{in}\mathsf{\Omega },\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill -\mathrm{\Delta }{v}_2+F_2(v_1,v_2)& =f_2,\mathrm{in}\mathsf{\Omega },\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial v_i}{\partial n}}& =0,\mathrm{on}\partial \mathsf{\Omega },i=1,2,\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial v_1}{\partial x}}+{\displaystyle \frac{\partial v_2}{\partial y}}& =0,\mathrm{in}\mathsf{\Omega },\hfill \end{array}$$

(8)

with the nonlinear convective terms defined by:

$$\begin{array}{cc}\hfill F_1(v_1,v_2)& ={\displaystyle \frac{1}{\mathrm{Re}}}\left(v_1{\displaystyle \frac{\partial v_1}{\partial x}}+v_2{\displaystyle \frac{\partial v_1}{\partial y}}\right),\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill F_2(v_1,v_2)& ={\displaystyle \frac{1}{\mathrm{Re}}}\left(v_1{\displaystyle \frac{\partial v_2}{\partial x}}+v_2{\displaystyle \frac{\partial v_2}{\partial y}}\right).\hfill \end{array}$$

(10)

Here, $\mathrm{Re}$ is the Reynolds number, a dimensionless parameter characterizing the ratio between inertial and viscous forces in the fluid:

$$\mathrm{Re}={\displaystyle \frac{\rho UL}{\nu }},$$

(11)

where $\rho$ is the fluid density, U a characteristic velocity, L a characteristic length, and $\nu$ the dynamic viscosity. Low Reynolds numbers indicate diffusion-dominated flows, whereas higher values suggest convection-dominated regimes with stronger mixing effects.

The model is based on the following:

• The fluid is Newtonian, incompressible, and its velocity field is stationary (i.e., independent of time).
• The pollutant is passive: it does not affect the flow field (one-way coupling).
• There is no pollutant exchange across the domain boundaries (homogeneous Neumann conditions).
• The domain $\mathsf{\Omega }$ may be geometrically complex, but is fixed in time and two-dimensional.

These simplifications make the model tractable while retaining the essential features of pollutant transport dynamics in realistic environmental systems such as shallow rivers or estuarine channels.

2. Discretized Model and Numerical Scheme

To numerically solve the coupled system described in Section 2, we employ a hybrid spatial discretization strategy that combines the Finite Element Method (FEM) for geometric flexibility and a pseudospectral method for high accuracy on each element. The resulting semi-discrete system enables efficient simulation of nonlinear convection-diffusion-reaction phenomena.

The computational implementation is carried out in MATLAB^®, which provides a flexible environment for assembling local-to-global transformations, handling sparse matrices, and integrating nonlinear dynamical systems.

2.1. Spatial Discretization Framework

Let ${\mathcal{T}}_h$ be a conforming triangulation of the domain $\mathsf{\Omega }$, consisting of $n_t$ triangular elements $T^k$. Each physical triangle $T^k$ is associated with an affine mapping ${\mathcal{H}}^k$ from the reference triangle

$$T_{\mathrm{ref}}=\left\{(r,s),-1\le r,s\le 1;r+s\le 0\right\}.$$

(12)

On each element, we define a local spectral basis using Dubiner polynomials [5], which are orthogonal and well-suited for spectral element methods on triangles.

Given N a local approximation parameter, representing the degree of the interpolation polynomial on each triangle, the global basis functions $g_{mn}^k$ where $m,n\ge 0$ and $m+n\le N$ are constructed by mapping the local basis functions $g_{mn}$ from $T_{\mathrm{ref}}$ to $T^k$ and extending them by zero outside the element:

$$g_{mn}^k(x,y)=\left\{\begin{array}{cc}{g}_{mn}\left({\left({\mathcal{H}}^k\right)}^{-1}(x,y)\right),\hfill & \mathrm{if}(x,y)\in T^k,\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(13)

The spectral Dubiner basis on the reference triangle is defined by:

$$g_{ij}(r,s)=J_i^{0,0}\left({\displaystyle \frac{2(r+1)}{1-s}}-1\right)J_j^{2i+1,0}\left(s\right){(1-s)}^i,$$

(14)

where $J_n^{\alpha ,\beta }$ is the Jacobi polynomial of degree n and order $(\alpha ,\beta )$.

A schematic representation of the local-to-global mapping process and the assembly of the global matrices is shown in Figure 1.

2.2. Steady-State Velocity Field Computation via Spectral-Galerkin Method

To compute the steady-state velocity field $\mathbf{v}=(v_1,v_2)$, we discretize the stationary incompressible Navier-Stokes equations using a Galerkin approach. For each triangular element $T^k$, we formulate the weak form of the momentum equations (5) and (6), and approximate them by integrating against the local basis functions $g_{mn}^k$. This yields the following variational expression:

$${\int }_{T^k}\nabla v_i·\nabla g_{mn}^{k}dxdy+{\int }_{T^k}{F}_i(v_1,v_2)g_{mn}^{k}dxdy={\int }_{T^k}{f}_i{g}_{mn}^{k}dxdy,i=1,2.$$

(15)

By employing pseudo-spectral differentiation, we define local differentiation matrices ${\mathbf{D}}_x^k$ and ${\mathbf{D}}_y^k$ over the element $T^k$, such that for any function f defined on $T^k$, the derivatives at the collocation points are approximated as:

$${\mathbf{f}}_x={\mathbf{D}}_x^k\mathbf{f},{\mathbf{f}}_y={\mathbf{D}}_y^k\mathbf{f},$$

(16)

where $\mathbf{f}$, ${\mathbf{f}}_x$, and ${\mathbf{f}}_y$ are vectors containing the values of f and its partial derivatives with respect to x and y at the collocation points of $T^k$.

Additionally, we introduce a quadrature matrix $\mathbf{Q}$, which is a diagonal matrix whose entries correspond to quadrature weights associated with the collocation points on the reference triangle $T_{\mathrm{ref}}$. This enables the approximation of integrals over $T^k$ as:

$${\int }_{T^k}f(x,y)g(x,y)dxdy\approx det\left(J^k\right)·{\mathbf{f}}^T\mathbf{Q}\mathbf{g},$$

(17)

where $J^k$ is the Jacobian matrix associated with the affine transformation from the reference triangle $T_{\mathrm{ref}}$ to the physical element $T^k$.

By applying the pseudo-spectral differentiation formula (16) and the quadrature approximation (17) to the local weak formulation (15), we obtain the following discrete system for each velocity component:

$${\mathbf{A}}^k{\mathbf{v}}_i^k+\mathbf{Q}{\mathbf{F}}_i({\mathbf{v}}_1^k,{\mathbf{v}}_2^k)=\mathbf{Q}{\mathbf{f}}_i^k,i=1,2,$$

(18)

where the local stiffness matrix is defined as:

$${\mathbf{A}}^k={\left({\mathbf{D}}_x^k\right)}^T\mathbf{Q}{\mathbf{D}}_x^k+{\left({\mathbf{D}}_y^k\right)}^T\mathbf{Q}{\mathbf{D}}_y^k,$$

(19)

and the nonlinear convective term is discretized as:

$${\mathbf{F}}_i({\mathbf{v}}_1^k,{\mathbf{v}}_2^k)={\displaystyle \frac{1}{\mathrm{Re}}}\left(\mathrm{diag}\left({\mathbf{v}}_1^k\right){\mathbf{D}}_x^k{\mathbf{v}}_i^k+\mathrm{diag}\left({\mathbf{v}}_2^k\right){\mathbf{D}}_y^k{\mathbf{v}}_i^k\right).$$

(20)

Here, the vectors ${\mathbf{v}}_i^k$ and ${\mathbf{f}}_i^k$ contain the discrete values of the velocity components $v_i$ and source terms $f_i$ at the collocation points of $T^k$. The notation $\mathrm{diag}\left({\mathbf{v}}_i^k\right)$ refers to a diagonal matrix whose entries are those of the vector ${\mathbf{v}}_i^k$.

By assembling the local equations (18) over all triangular elements in the mesh, we obtain an assembled system in local numbering:

$$\mathbf{A}{\mathbf{v}}_i^{\mathrm{loc}}+{\mathbf{Q}}^{\mathrm{loc}}{\mathbf{F}}_i^{\mathrm{loc}}({\mathbf{v}}_1^{\mathrm{loc}},{\mathbf{v}}_2^{\mathrm{loc}})={\mathbf{Q}}^{\mathrm{loc}}{\mathbf{f}}_i^{\mathrm{loc}},i=1,2,$$

(21)

where:

• ${\mathbf{v}}_i^{\mathrm{loc}}$ is the concatenated vector obtained by stacking all local velocity vectors ${\mathbf{v}}_i^k$ from each triangle $T^k$;
• $\mathbf{A}$ is a block-diagonal stiffness matrix whose diagonal blocks are the local stiffness matrices ${\mathbf{A}}^k$;
• ${\mathbf{Q}}^{\mathrm{loc}}$ is a block-diagonal matrix with each diagonal block equal to the quadrature matrix $\mathbf{Q}$;
• ${\mathbf{F}}_i^{\mathrm{loc}}({\mathbf{v}}_1^{\mathrm{loc}},{\mathbf{v}}_2^{\mathrm{loc}})$ and ${\mathbf{f}}_i^{\mathrm{loc}}$ are the concatenated vectors of the local nonlinear convective terms ${\mathbf{F}}_i({\mathbf{v}}_1^k,{\mathbf{v}}_2^k)$ and the source terms ${\mathbf{f}}_i^k$, respectively.

All local contributions are related to the global degrees of freedom via a global transformation matrix $\mathbf{Z}$, such that:

$${\mathbf{v}}_i^{\mathrm{loc}}=\mathbf{Z}{\mathbf{v}}_i^{\mathrm{glob}},i=1,2,$$

(22)

where ${\mathbf{v}}_i^{\mathrm{glob}}$ is the vector of global unknowns for the i-th velocity component. Substituting the above transformation into the assembled local system (21), we obtain the global nonlinear system:

$$\mathbf{A}\mathbf{Z}{\mathbf{v}}_i^{\mathrm{glob}}+{\mathbf{Q}}^{\mathrm{loc}}{\mathbf{F}}_i^{\mathrm{loc}}(\mathbf{Z}{\mathbf{v}}_1^{\mathrm{glob}},\mathbf{Z}{\mathbf{v}}_2^{\mathrm{glob}})={\mathbf{Q}}^{\mathrm{loc}}{\mathbf{f}}_i^{\mathrm{loc}},i=1,2.$$

(23)

This system represents a large-scale nonlinear problem, where the unknown global velocity vectors ${\mathbf{v}}_1^{\mathrm{glob}}$ and ${\mathbf{v}}_2^{\mathrm{glob}}$ must also satisfy the discrete incompressibility condition (8) which is discretized using the same spectral differentiation framework. This results in the following global discrete constraint:

$${\mathbf{D}}_x\mathbf{Z}{\mathbf{v}}_1^{\mathrm{glob}}+{\mathbf{D}}_y\mathbf{Z}{\mathbf{v}}_2^{\mathrm{glob}}=\mathbf{0},$$

(24)

where ${\mathbf{D}}_x$ and ${\mathbf{D}}_y$ are block-diagonal matrices that assemble all local spectral differentiation matrices ${\mathbf{D}}_x^k$ and ${\mathbf{D}}_y^k$, respectively.

2.3. Pollutant Concentration Update

Once the steady-state velocity field $\mathbf{v}=(v_1,v_2)$ is computed, the pollutant transport equation is discretized in space using the same local spectral basis and quadrature framework as for the velocity field.

By multiplying the governing equation (1) by the basis function $g_{mn}^k$ and integrating over a triangle $T^k$, we obtain the weak formulation:

$${\displaystyle \frac{d}{dt}}{\int }_{T^k}C{g}_{mn}^{k}dxdy+\lambda {\int }_{T^k}\nabla C·\nabla g_{mn}^{k}dxdy+{\int }_{T^k}(\mathbf{v}·\nabla C)g_{mn}^{k}dxdy+{\int }_{T^k}F\left(C\right)g_{mn}^{k}dxdy=0.$$

(25)

Proceeding as in the velocity formulation, the advective term is approximated using spectral differentiation and quadrature:

$${\int }_{T^k}(\mathbf{v}·\nabla C)g_{mn}^{k}dxdy\approx det\left(J^k\right)·{\left({\mathbf{g}}_{mn}^k\right)}^T\mathbf{Q}\left[\mathrm{diag}\left({\mathbf{v}}_1^k\right){\mathbf{D}}_x^k+\mathrm{diag}\left({\mathbf{v}}_2^k\right){\mathbf{D}}_y^k\right]{\mathbf{C}}^k,$$

(26)

where ${\mathbf{C}}^k$ is the vector of pollutant concentration values at the collocation points of $T^k$.

The resulting semi-discrete approximation of (25) reads:

$${\left({\mathbf{g}}_{mn}^k\right)}^T\mathbf{Q}{\displaystyle \frac{d}{dt}}{\mathbf{C}}^k+\lambda {\left({\mathbf{g}}_{mn}^k\right)}^T{\mathbf{A}}^k{\mathbf{C}}^k+{\left({\mathbf{g}}_{mn}^k\right)}^T\mathbf{Q}{\mathbf{B}}^k({\mathbf{v}}_1^k,{\mathbf{v}}_2^k){\mathbf{C}}^k+{\left({\mathbf{g}}_{mn}^k\right)}^T\mathbf{Q}\mathbf{F}\left({\mathbf{C}}^k\right)=0,$$

(27)

where we define:

$${\mathbf{B}}^k({\mathbf{v}}_1^k,{\mathbf{v}}_2^k)=\mathrm{diag}\left({\mathbf{v}}_1^k\right){\mathbf{D}}_x^k+\mathrm{diag}\left({\mathbf{v}}_2^k\right){\mathbf{D}}_y^k,$$

(28)

and where $\mathbf{F}\left({\mathbf{C}}^k\right)$ denotes the evaluation of the nonlinear reaction term $F\left(C\right)$ at the collocation points of $T^k$.

Simplifying (27) and isolating the time derivative yields the local system of ODEs:

$${\displaystyle \frac{d}{dt}}{\mathbf{C}}^k=-\lambda {\mathbf{Q}}^{-1}{\mathbf{A}}^k{\mathbf{C}}^k-{\mathbf{B}}^k({\mathbf{v}}_1^k,{\mathbf{v}}_2^k){\mathbf{C}}^k-\mathbf{F}\left({\mathbf{C}}^k\right).$$

(29)

Aggregating these local systems for all $k=1,\dots ,n_t$, and using the same block matrix notation introduced earlier, we obtain the global system in local numbering:

$${\displaystyle \frac{d}{dt}}{\mathbf{C}}^{\mathrm{loc}}=-\lambda {\mathbf{A}}^{\mathrm{diff},\mathrm{loc}}{\mathbf{C}}^{\mathrm{loc}}-\mathbf{B}({\mathbf{v}}_1^{\mathrm{loc}},{\mathbf{v}}_2^{\mathrm{loc}}){\mathbf{C}}^{\mathrm{loc}}-\mathbf{F}\left({\mathbf{C}}^{\mathrm{loc}}\right),$$

(30)

where:

• ${\mathbf{A}}^{\mathrm{diff},\mathrm{loc}}$ is a block-diagonal matrix whose k-th diagonal block is ${\mathbf{Q}}^{-1}{\mathbf{A}}^k$;
• $\mathbf{B}({\mathbf{v}}_1^{\mathrm{loc}},{\mathbf{v}}_2^{\mathrm{loc}})$ is a block-diagonal matrix with blocks ${\mathbf{B}}^k({\mathbf{v}}_1^k,{\mathbf{v}}_2^k)$.

Finally, transforming this system into global numbering using the global-to-local mapping matrix $\mathbf{Z}$ and the local-to-global projection matrix $\mathbf{W}$, we arrive at the global semi-discrete formulation:

$${\displaystyle \frac{d}{dt}}{\mathbf{C}}^{\mathrm{glob}}=-\lambda \mathbf{W}{\mathbf{A}}^{\mathrm{diff},\mathrm{loc}}\mathbf{Z}{\mathbf{C}}^{\mathrm{glob}}-\mathbf{W}\mathbf{B}({\mathbf{v}}_1^{\mathrm{loc}},{\mathbf{v}}_2^{\mathrm{loc}})\mathbf{Z}{\mathbf{C}}^{\mathrm{glob}}-\mathbf{W}\mathbf{F}\left(\mathbf{Z}{\mathbf{C}}^{\mathrm{glob}}\right).$$

(31)

The semi-discrete problem in pollutant concentration thus consists of computing the global time-dependent solution ${\mathbf{C}}^{\mathrm{glob}}\left(t\right)$ over the time interval $[0,T]$ subject to the initial condition:

$${\mathbf{C}}^{\mathrm{glob}}\left(0\right)={\mathbf{C}}_0^{\mathrm{glob}},$$

(32)

where ${\mathbf{C}}_0^{\mathrm{glob}}$ is the vector of initial concentration values evaluated at the global collocation points.

2.4. Convergence, Stability, and Algorithmic Implementation

2.4.1. Convergence of the Discretized Velocity Field Problem

The discretization of the stationary incompressible Navier-Stokes equations using a spectral-Galerkin approach on triangular elements results in a nonlinear saddle-point problem, constrained by the discrete incompressibility condition. The convergence of the velocity field approximation ${\mathbf{v}}^{\mathrm{glob}}$ hinges on two main aspects: the spectral approximation properties of the Dubiner basis functions and the regularity of the exact solution. Assuming that the exact velocity field $\mathbf{v}$ lies in a Sobolev space ${\left[H^s\left(\mathsf{\Omega }\right)\right]}^2$ for $s>1$, the spectral element discretization achieves exponential convergence in space with respect to the polynomial degree N on each triangle, provided the mesh ${\mathcal{T}}_h$ is quasi-uniform and the geometry is smooth or polygonal [7,9].

More precisely, let ${\mathbf{v}}_N$ denote the discrete spectral solution obtained from system (23) and $\mathbf{v}$ the exact solution of the continuous weak problem. we can state the following result:

Theorem 1.

Let $\mathsf{\Omega }\subset {\mathbb{R}}^2$ be a polygonal domain, and suppose the exact velocity field $\mathbf{v}\in {\left[H^s\left(\mathsf{\Omega }\right)\right]}^2$ for some $s>1$. Let ${\mathbf{v}}_N$ denote the discrete velocity field obtained from the spectral-Galerkin system (23) using polynomial degree N on each triangle. Then, there exists a constant $C>0$, independent of N, such that:

$$\parallel \mathbf{v}-{\mathbf{v}}_N{\parallel }_{H^1\left(\mathsf{\Omega }\right)}\le C{N}^{-s+1}{\parallel \mathbf{v}\parallel }_{H^s\left(\mathsf{\Omega }\right)}.$$

(33)

If $\mathbf{v}$ is analytic on Ω, then the convergence is spectral, i.e., exponentially fast:

$$\parallel \mathbf{v}-{\mathbf{v}}_N{\parallel }_{H^1\left(\mathsf{\Omega }\right)}\le C{e}^{-bN},$$

(34)

for some $b>0$.

Proof. 

The proof follows classical approximation theory for spectral element methods [7,9]. On each triangle $T^k$, the spectral basis formed by Dubiner polynomials yields an approximation error of order $N^{-s}$ in $L^2$ and $N^{-s+1}$ in $H^1$, due to inverse inequalities. The assembly of the global error estimate then follows from a Céa-type lemma applied to the discrete weak formulation of (15), using the inf-sup stability of the discrete spaces under suitable projection or stabilization. The exponential convergence in the analytic case comes from the classical result on orthogonal polynomial approximation of analytic functions on triangles.    □

2.4.2. Convergence, Stability, and Algorithmic Implementation of the Semi Discretized Pollutant Problem

The semi-discrete formulation derived above leads to a system of nonlinear ordinary differential equations (ODEs) for the pollutant concentration ${\mathbf{C}}^{\mathrm{glob}}\left(t\right)$. To guarantee the well-posedness and convergence of the numerical solution, we must ensure stability and boundedness of the solution over time.

2.4.2.1. Energy Boundedness Condition

We assume that the nonlinear reaction term $F\left(C\right)$ satisfies the following dissipative condition:

$$F\left(C\right)·C\ge \gamma C^2,\mathrm{for}\mathrm{some}\gamma \ge 0.$$

(35)

Under this assumption, one can derive an energy estimate by multiplying equation (31) by ${\mathbf{C}}^{\mathrm{glob},T}$ and applying symmetry properties of the matrices involved. This yields:

$${\displaystyle \frac{1}{2}}{\displaystyle \frac{d}{dt}}\parallel {\mathbf{C}}^{\mathrm{glob}}{\parallel }_{{\mathbf{Q}}^{\mathrm{glob}}}^2\le -\lambda \parallel \nabla {\mathbf{C}}^{\mathrm{glob}}{\parallel }^2+\gamma {\parallel {\mathbf{C}}^{\mathrm{glob}}\parallel }^2,$$

(36)

where ${\parallel ·\parallel }_{{\mathbf{Q}}^{\mathrm{glob}}}$ is the discrete weighted $L^2$ norm. By applying Grönwall’s inequality, this implies boundedness of the solution for all $t\in [0,T]$ provided that $\gamma$ remains small compared to the diffusion and advection damping rates.

2.4.2.2. Time Discretization and Solver Strategy

Rather than designing an ad-hoc time integrator, we advance the semi-discrete system obtained after spatial discretization by calling MATLAB’s built-in adaptive ODE solver. Denote by

$${\displaystyle \frac{d}{dt}}{\mathbf{C}}^{\mathrm{glob}}\left(t\right)=\mathcal{G}\left(t,{\mathbf{C}}^{\mathrm{glob}}\left(t\right)\right):=-\lambda \mathbf{W}{\mathbf{A}}^{\mathrm{diff},\mathrm{loc}}\mathbf{Z}{\mathbf{C}}^{\mathrm{glob}}-\mathbf{W}\mathbf{B}({\mathbf{v}}_1^{\mathrm{loc}},{\mathbf{v}}_2^{\mathrm{loc}})\mathbf{Z}{\mathbf{C}}^{\mathrm{glob}}-\mathbf{W}\mathbf{F}\left(\mathbf{Z}{\mathbf{C}}^{\mathrm{glob}}\right),$$

with initial condition ${\mathbf{C}}^{\mathrm{glob}}\left(0\right)={\mathbf{C}}_0^{\mathrm{glob}}$; the steady velocity field ${\mathbf{v}}^{\mathrm{loc}}=({\mathbf{v}}_1^{\mathrm{loc}},{\mathbf{v}}_2^{\mathrm{loc}})$ is precomputed once and kept fixed during the concentration update. We then solve the IVP

$${\displaystyle \frac{d}{dt}}{\mathbf{C}}^{\mathrm{glob}}\left(t\right)=\mathcal{G}\left(t,{\mathbf{C}}^{\mathrm{glob}}\left(t\right)\right),t\in [0,T],{\mathbf{C}}^{\mathrm{glob}}\left(0\right)={\mathbf{C}}_0^{\mathrm{glob}},$$

by invoking MATLAB’s adaptive solver ode15s. This black-box approach automatically controls the time step to meet user-specified tolerances, and it avoids any custom IMEX design.

2.4.2.3. MATLAB Solver-Based Workflow

The implementation proceeds as follows:

• Precompute constant operators: assemble sparse matrices $\mathbf{L}:=\mathbf{W}{\mathbf{A}}^{\mathrm{diff},\mathrm{loc}}\mathbf{Z}$ and $\mathbf{B}\left({\mathbf{v}}^{\mathrm{loc}}\right)$ (kept fixed if the velocity is steady), and prepare a function for the nonlinear reaction vector $\mathbf{F}(·)$ evaluated at $\mathbf{Z}{\mathbf{C}}^{\mathrm{glob}}$.
• Define the RHS callback for the ODE solver:

  $$\mathrm{RHS}(t,\mathbf{C})=-\lambda \mathbf{L}\mathbf{C}-\mathbf{W}\mathbf{B}\left({\mathbf{v}}^{\mathrm{loc}}\right)\mathbf{Z}\mathbf{C}-\mathbf{W}\mathbf{F}\left(\mathbf{Z}\mathbf{C}\right).$$
• Call the solver with the desired tolerances:

  $$[t,\mathbf{C}\left(t\right)]=\mathtt{ode}\mathtt{15}\mathtt{s}(\mathrm{RHS},[0,T],{\mathbf{C}}_0^{\mathrm{glob}},\mathtt{odeset}(\text{’}\mathrm{RelTol}\text{’},{10}^{-6},\text{’}\mathrm{AbsTol}\text{’},{10}^{-9})).$$
• Post-process$\mathbf{C}\left(t\right)$ at requested output times, and map to local or physical fields via $\mathbf{Z}$ when needed.

3. Numerical Experiments

This section presents numerical experiments to validate and illustrate the performance of the proposed spectral-Galerkin framework for solving the coupled Navier–Stokes and pollutant transport equations. The simulations are implemented in MATLAB R2024a, leveraging sparse matrix operations and customized local-to-global assembly routines. Two main components are addressed: the steady-state velocity field computation and the transient pollutant concentration evolution.

The simulations are implemented in a modular MATLAB framework. Key features include:

• Local-to-global transformation via sparse matrix $\mathbf{Z}$.
• Assembly of spectral differentiation matrices ${\mathbf{D}}_x^k$, ${\mathbf{D}}_y^k$ using Jacobi or Dubiner basis.
• Block-diagonal structures for stiffness $\mathbf{A}$ and convection matrices $\mathbf{B}$ for efficient matrix- vector operations.
• Divergence enforcement via projection on $ker({\mathbf{D}}_x\mathbf{Z},{\mathbf{D}}_y\mathbf{Z})$.
• Time integration for concentration handled with vectorized Runge–Kutta schemes.

3.1. Numerical Resolution Strategy for the Discrete Velocity Field

The discrete system governing the steady-state velocity field, given in (23) together with the incompressibility constraint (24), results in a large-scale nonlinear algebraic problem. Direct solvers quickly become computationally prohibitive as the number of elements and polynomial degree increase, due to the size and sparsity pattern of the global matrices. To overcome this limitation, we adopt an efficient iterative scheme based on gradient methods, implemented in MATLAB.

3.1.0.4. Linearization

At each nonlinear iteration, the convective terms ${\mathbf{F}}_i^{\mathrm{loc}}(\mathbf{Z}{\mathbf{v}}_1^{\mathrm{glob}},\mathbf{Z}{\mathbf{v}}_2^{\mathrm{glob}})$ are linearized around the current velocity approximation ${\mathbf{v}}^{\left(m\right)}$. This yields a Jacobian-based correction system of the form:

$$\mathbf{J}\left({\mathbf{v}}^{\left(m\right)}\right)\delta \mathbf{v}=-\mathbf{R}\left({\mathbf{v}}^{\left(m\right)}\right),$$

(37)

where $\mathbf{J}$ is the Jacobian matrix of the nonlinear operator, $\delta \mathbf{v}$ the correction, and $\mathbf{R}$ the nonlinear residual.

3.1.0.5. Gradient-Based Solver

The correction system is solved iteratively using a preconditioned conjugate gradient (PCG) method. This choice exploits the symmetry and sparsity of the discretized operators, significantly reducing memory requirements and computation time. In MATLAB, this is implemented via the built-in operator solve (equivalent to mldivide) applied to the sparse Jacobian system, with the gradient iterations handled internally. For large problems, additional preconditioners such as incomplete LU (ILU) or Jacobi can be specified.

3.1.0.6. Algorithm

The overall resolution strategy is summarized as follows:

• Initialize with an admissible velocity field ${\mathbf{v}}^{\left(0\right)}$ satisfying the discrete incompressibility constraint (24).
• At iteration m, assemble the nonlinear residual $\mathbf{R}\left({\mathbf{v}}^{\left(m\right)}\right)$ and the Jacobian $\mathbf{J}\left({\mathbf{v}}^{\left(m\right)}\right)$.
• Solve the linear correction system

  $$\mathbf{J}\left({\mathbf{v}}^{\left(m\right)}\right)\delta \mathbf{v}=-\mathbf{R}\left({\mathbf{v}}^{\left(m\right)}\right)$$

  using MATLAB’s solve with a gradient-based iterative scheme.
• Update the velocity field: ${\mathbf{v}}^{(m+1)}={\mathbf{v}}^{\left(m\right)}+\delta \mathbf{v}$.
• Enforce the incompressibility condition (24) by projecting the updated velocity field onto the divergence-free subspace.
• Repeat until convergence: $\parallel \delta \mathbf{v}\parallel /\parallel {\mathbf{v}}^{\left(m\right)}\parallel <\epsilon$, with tolerance $\epsilon \sim {10}^{-8}$.

3.1.0.7. Implementation in MATLAB

The following schematic illustrates the key commands:

J = assembleJacobian(v);

R = assembleResidual(v);

delta_v = J \ (-R);    % equivalent to solve(J, -R) in MATLAB

v_new = v + delta_v;

Thanks to MATLAB’s efficient sparse linear algebra routines, the above approach enables the rapid computation of the velocity field even on fine meshes and for higher-order polynomial approximations.

This gradient-based resolution strategy ensures both computational efficiency and numerical stability, and is integrated into the global simulation framework presented in Section 3.2.

3.2. Manufactured Velocity Field for Convergence Verification

To validate the convergence properties of the method, we consider a manufactured solution defined on $\mathsf{\Omega }={(-1,1)}^2$:

$$v_{1\mathrm{ex}}(x,y)=cos\left(\frac{\pi x}{2}\right)sin\left(\frac{\pi y}{2}\right),v_{2\mathrm{ex}}(x,y)=-sin\left(\frac{\pi x}{2}\right)cos\left(\frac{\pi y}{2}\right),$$

which is divergence-free and satisfies homogeneous Neumann conditions on the boundary. The corresponding source term $\mathbf{f}=(f_1,f_2)$ is computed by substituting ${\mathbf{v}}_{\mathrm{ex}}=(v_{1\mathrm{ex}},v_{2\mathrm{ex}})$ into the steady Navier–Stokes equations.

The numerical velocity field ${\mathbf{v}}_N$ is computed on pseudo spectal collocation points on triangular meshes using polynomial degree $N=4,6,8,\dots$. The $L^2$ errors between ${\mathbf{v}}_N$ and $\mathbf{v}$ the vector values of the exact solution at pseudo spectral collocation points are evaluated. Table 1 reports the exponential decay of the errors, confirming spectral convergence.

A convergence plot (Figure 2) confirms the exponential accuracy of the scheme with respect to the polynomial order N.

3.3. Numerical Simulation of the Coupled Velocity–Concentration System

To assess the performance of the proposed spectral–Galerkin discretization in a realistic setting, we simulate the coupled velocity and pollutant transport problem described in Section 2. The computational domain $\mathsf{\Omega }$ corresponds to a river channel containing a central islet, which introduces geometric complexity and induces heterogeneous flow patterns. The domain is discretized into 119 triangular spectral elements using a mesh generator (see Figure 3). Unless otherwise specified, we fixe $\mu =1$, while the influence of the nonlinear reaction parameter $\beta$ and the advection strength parameter $\lambda$ are systematically investigated. The polynomial order of the pseudospectral approximation is set to $N=4$ throughout all computations. The Reynolds number is prescribed as $\mathrm{Re}=100$, and the forcing term is identical to that employed in the steady-state velocity computation described earlier.

Simulation Setup

The simulation proceeds in two stages. First, the steady-state velocity field is computed by solving the nonlinear system (23) subject to the incompressibility constraint (24). The resulting velocity profiles are then used as input data for the transient pollutant transport problem governed by the semi-discrete system (31). Time integration of the concentration field is performed with a second-order explicit Runge–Kutta method, which provides a good balance between accuracy and computational efficiency. Homogeneous Neumann boundary conditions are prescribed for the concentration, enforcing a no-flux condition at the riverbanks and islet boundary. The pollutant is initially localized in the vicinity of a prescribed interior point of $\mathsf{\Omega }$, mimicking the instantaneous release of a contaminant source.

Test Scenarios

Several test cases are conducted to quantify the effect of nonlinear reaction dynamics. In particular, the role of the parameter $\beta$ in shaping the pollutant distribution is examined by varying its value over a wide range, while $\lambda$ controls the intensity of pollutant diffusion. Unless otherwise noted, a reference configuration with $\beta =1.5$ and $\lambda =0.001$ is considered.

Results: Velocity Field

Figure 4 illustrates the steady-state velocity field for $\mathrm{Re}=100$. The computed streamlines highlight the strong influence of the islet on the overall flow structure, with the formation of recirculation zones and regions of locally accelerated flow. For higher Reynolds numbers (e.g., $\mathrm{Re}=1000$), the velocity gradients become sharper, leading to enhanced advective transport in the downstream region.

Results: Pollutant Transport

The pollutant transport simulations reveal the interplay between advection, diffusion, and nonlinear reaction. Figure 5 shows the temporal evolution of the concentration field for $\beta =1.5$. The pollutant plume initially spreads due to local diffusion, while advection progressively transports the contaminant along streamlines. Over time, nonlinear reactions amplify concentration gradients, leading to heterogeneous accumulation patterns within the domain.

The influence of the nonlinear reaction parameter $\beta$ is further analyzed in Figure 6 and Figure 7, which display the pollutant concentration at $t=1.3$ for $\beta \in \{2,3,4,5\}$. Increasing $\beta$ results in pronounced pollutant accumulation near the riverbanks and around the islet boundary, demonstrating the critical role of reaction terms in shaping long-term pollutant patterns.

Discussion

The numerical experiments confirm the robustness and accuracy of the proposed spectral–Galerkin discretization for coupled velocity–concentration problems. The method remains stable for complex geometries and polynomial degrees up to $N=4$. The results highlight the following major physical effect:

• Increasing $\beta$ enhances pollutant accumulation near boundaries and reactive hotspots, reflecting the dominant role of nonlinear reaction kinetics.

Overall, the spectral–Galerkin approach proves effective for capturing intricate nonlinear interactions between fluid flow and pollutant dynamics, making it a promising tool for large-scale environmental flow simulations.

4. Conclusions

In this work, we developed a hybrid numerical framework that combines the Finite Element Method (FEM) and Pseudospectral techniques to simulate the nonlinear diffusion and advection of pollutants in a two-dimensional aquatic environment. The proposed method addresses the challenges posed by complex geometries and nonlinear reaction terms by leveraging the geometric flexibility of FEM and the high resolution of spectral discretizations.

We formulated a coupled system involving the stationary incompressible Navier-Stokes equations for fluid velocity and a nonlinear convection-diffusion-reaction equation for pollutant concentration. A detailed numerical scheme was presented, along with two modular algorithms for solving the velocity and concentration fields efficiently. Our simulations demonstrated the robustness and accuracy of the method under varying Reynolds numbers and nonlinear reaction parameters.

Key observations include the method’s numerical stability across a range of discretization parameters and the pronounced sensitivity of pollutant accumulation patterns to the nonlinearity parameter $\beta$. These insights underscore the importance of accurately modeling reaction kinetics in water quality simulations.

Future work may extend this framework to three-dimensional domains, incorporate time-dependent velocity fields, or account for multi-species transport and reactive boundary conditions. Furthermore, integrating real-world data or coupling with watershed models could enhance the method’s applicability in environmental decision-making and pollution mitigation strategies.