# A Novel Multivariate Spectral Local Quasilinearization Method (MV-SLQLM) for Modelling Flow, Moisture, Heat, and Solute Transport in Soil

Authors: Elias Mwakilama; Vusi Magagula; Duncan Gathungu
Journal: Journal of Applied Mathematics
 (2023)

---

## Abstract

Conventionally, the problem of studying the transport of water, heat, and solute in soil or groundwater systems has been numerically solved using finite difference (FD) or finite element (FE) methods. FE methods are attractive over FD methods because they are geometrically flexible. However, recent studies demonstrate that spectral collocation (SC) methods converge exponentially faster than FD or FE methods using a few grid points or on coarse grids. This work proposes and applies a multivariate spectral local quasilinearization method (MV-SLQLM) to model the transportation and interaction of soil moisture, heat, and solute concentration in a nonbare soil ridge. The MV-SLQLM uses a quasilinearization method (QLM) to linearize any nonlinear equations and then employs a local linearization method (LLM) to decouple the linearized system of PDEs to form a sequence of equations that are solved in a computationally efficient manner. The MV-SLQLM is thus an extension of the bivariate spectral local linearization method (BI-SLLM) that fails to deal with a 2D problem and is a modification of the MV-SQLM whose efficiency is compromised when operating on high dense solution matrices. We use the residual error norms of the difference between successive iterations to affirm convergence to the expected solution. To illustrate the application and check the solution accuracy, we conduct systematic analyses of the effect of model parameters on distribution profiles. Findings are in good agreement with theory and literature, thereby revealing suitability of the MV-SLQLM to solve coupled nonlinear PDEs with environmental fluid dynamics applications.


---

## Body

## 1. Introduction

The problem of modelling the transport and interaction of heat and mass in porous media is of interest to the extent that there is a considerable body of work on analytical, semianalytical, and numerical solutions. Because of their ease of use and low costs at implementation, analytical and semianalytical solutions are considered efficient tools for modelling the fate and transport of heat, water, or contaminants in soil or groundwater systems. Among others [1–3], Shao et al. [4] derived analytical solutions for heat transport equations in soil media. Including [5, 6], Chen et al. [7] proposed analytical solutions to linearized Richards equations for describing infiltration. Sanskrityayn et al. [8] and other studies [9–18] provided analytical solutions to modelling the transport of contaminants in soil and groundwater systems. Other studies [19–21] proposed analytical solutions to solute and water transport equations for verifying numerical solvers. Gao et al. [22] and others [23–25] employed semianalytical solutions to model the transport of solutes and water in soil and aquifers. However, the main limitation with analytical and semianalytical solutions is their inability to handle complex nonlinear models. Moreover, besides a few citable examples, many analytical tools are not multidimensional and do not consider interaction. Numerical models can handle such complexity but at the expense of defining several parameters that are often unknown, leading to model outputs that are highly uncertain [9]. In addition, numerical solutions have natural approximation errors, can be computationally demanding, and may generate very different results using the same input parameters [9].Nonetheless, several numerical solutions have been proposed to model heat, water, and solute transfer in porous media using finite difference (FD), finite volume (FV), or finite element (FE) methods. The solutions are problem-specific or as general commercial software [26]. Ehrhardt and Mickens [27] and Wang and Roeger [28] proposed nonstandard FD schemes to solve convection-diffusion equations with constant coefficients and convection-diffusion-reaction equations. Soko et al. [29] studied the inverse problem of identifying multiple pollution sources in groundwater systems using a FV method. Noborio et al. [30, 31] developed a two-dimensional (2D) Galerkin FE model for water, heat, and solute transport in soil ridge. Benjamin et al. [32] studied the effect of ridge size and geometry on simulated water and heat transport in soil ridge using a FE model. Baghlani and Nikzad [33] proposed a FD-differential quadrature method for solving 2D unsaturated water flow. Juncu et al. [34] proposed a nonlinear multigrid FD method (FDM) for the numerical solution of the unsaturated flow equation in 2D. Abbas et al. [35] studied a Darcy-Brinkman-Forchheimer model for the laminar natural convection in porous media using a FE method (FEM). Kelleners et al. [36] proposed a FE numerical model of coupled water flow and heat transport in soil and snow. Bennett [37] proposed a FD scheme to solve a Richards’ equation. Reichardt and Timm [38] described water redistribution after infiltration into the soil using a FE tool. Ren et al. [39] investigated the effects of temperature gradient on soil water-salt transfer with evaporation using an FD model. Cimolin and Discacciati [40] studied the Navier-Stokes/Forchheimer model for infiltration through porous media using FE. Fujisawa and Murakami [41] employed FEM to numerically analyze coupled flows in porous and fluid domains using the Darcy-Brinkman equations. Celia et al. [42] proposed a general mass-conservative FE numerical solution coupled with Picard iteration for the unsaturated flow equation. Using FEM, Li et al. [43] proposed a numerical homogenization of the Richards equation for unsaturated water flow through heterogeneous soils. Šimunek et al. [44] proposed an HYDRUS model for simulating the movement of water, heat, and multiple solutes in variably saturated media. Fujimaki [45] and Watanabe [46] applied the Hydrus to study irrigation amounts and model water and heat flow in a directionally frozen silt soil, respectively. Haile [47] proposed a VS2DRR FE-based tool to model the transport of reactive contaminants in the vadose zone.In spite of the aforementioned FD, FV, or FE numerical applications, one major setback with the FD methods is their demand for several or fine grid points to generate approximate solutions of quadratic convergence [48]. Although the FE methods are geometrically flexible and therefore attractive over FD or FV methods [30], classical FE mesh needs to be fine to model inhomogeneity surface material accurately [49]. Consequently, in both FD and FE methods, demand for fine meshes implies excessive use of high number of elements and costly computational means. Despite some advances in the performance of FE methods, through use of h-p or Galerkin versions [50, 51], the development of alternative numerical schemes for efficient and accurate heat and mass transfer analysis remains a subject of interest.Spectral methods have emerged as attractive numerical techniques in Engineering and Science [52–55]. As demonstrated on various theoretical diffusion-reaction or advection-diffusion problems [56–58], spectral collocation methods (SCM) are more efficient than FD or FE methods. Unlike the FDM and FEM, the SCM converge exponentially fast using a few grid points and on coarse grids at a minimal computational cost [59]. Motsa [60] developed a new spectral local linearization method (SLLM) for solving nonlinear boundary layer flow problems. Later, the SLLM was extended to a BI-SLLM [61]. Thumma and Magagula [62] studied the transient electromagnetohydrodynamic radiative squeezing flow between two parallel Riga plates using the SLLM. Motsa [63] presented and compared the spectral relaxation method (SRM) and spectral local quasilinearization method (SQLM) for solving unsteady boundary layer flow problems. Ansari et al. [64] studied the Jeffrey nanofluid flow near a Riga plate using the SQLM. Motsa [63] and Magagula et al. [65] proposed a bivariate-SQLM (BI-SQLM) to solve nonlinear PDEs. Jamil et al. [66] studied the unsteady MHD micropolar-nanofluids with homogeneous-heterogeneous chemical reactions over a stretching surface using the BI-SLLM. Since both BI-SQLM and BI-SLLM face challenges with problems of large time scales, Motsa et al. [67] and Magagula et al. [68] proposed the nonoverlapping multidomain spectral collocation methods, i.e., multidomain BI-SLLM [69] and multidomain BI-SQLM [65, 70]. Later, Mkhatshwa et al. [71] proposed an overlapping multidomain spectral method to study conjugate problems of conduction and MHD free convection flow of nanofluids over flat plates. Dlamini and Magagula [72] proposed a multivariate SQLM (MV-SQLM) to solve 2D Burger’s equations. Lekoko et al. [73] and Goqo et al. [74] applied the MV-SQLM to analyze pressure and heat distribution in a filter chamber and to study heat and mass transfer in a porous media saturated with nanofluid, respectively. Related to the MV-SQLM is the trivariate spectral collocation method (TV-SCM) proposed by Samuel and Motsa [75]. The MV-SQLM is an extension of the BI-SQLM as the latter does not handle 2D problems. However, the BI-SQLM and MV-SQLM generate more dense solution matrices when solving a coupled system of ODEs or PDEs, thereby compromising numerical efficiency. The BI-SLLM, however, employs the LLM to decouple the system and iteratively solve the problem on low-dense solution matrices, but for 1D problems.Despite the case, performance analysis of the previous spectral variants has depended heavily upon theoretical or practical nanofluid type of problems [59], whose medium of transmission differs from that of soil. Soil structure presents complex modelling processes that lead to highly coupled flow governing equations. In this work, we propose and apply a multivariate spectral local quasilinearization method (MV-SLQLM), extending the BI-SLLM and modifying the MV-SQLM, to model the transport and interaction of heat and mass in a 2D homogeneous porous medium. The MV-SLQLM uses the quasilinearization method (QLM) of Bellman et al. [76] to linearize the nonlinear equations and the local linearization method (LLM) [60, 61] to decouple the system and efficiently solve the problem on low-dense matrices. The proposed solution method is motivated by a practical problem of understanding transport and interaction of soil moisture, soil temperature, and solute concentration in a soil ridge with water in the furrows. Previously, Benjamin et al. [32] investigated the bare ridge tillage (height) effect on the simulated water and heat transport, disregarding presence of solute. Noborio et al. [30] studied the interaction of water, heat, and solute in furrow-irrigated soil, but with a focus on the effect of bare soil ridge surface and furrow. Our work assumes nonbare soil ridge to include the effect of plant-root uptake activity. Using the proposed MV-SLQLM, this study generates knowledge of the simultaneous distribution and interaction of soil water, moisture, heat, and solute in a nonbare soil ridge, useful for the proper management of ridge-tillage soil field to enhance crop production.The rest of the paper is organized as follows: in Section2, we present the investigated domain and the governing equations. In Section 3, we describe the solution method, error, and stability analysis and illustrate how the method is employed to solve the mentioned problem. Simulation results are presented and discussed in Section 4, along with convergence rate of the proposed method. Finally, Section 5 summarizes and concludes the current study.


## 2. Investigated Domain and Governing Equations

Guided by the characteristics of natural convection flow in a homogeneous porous medium, this study considers a problem of environmental application that describes the transport and interaction of soil moisture, heat, and solute concentration in a nonbare soil ridge when liquid water is present in the furrows (Figure1). Knowledge of interactive transport mechanisms inside a nonbare soil ridge is vital to agronomists and environmental researchers for upholding good environmental conditions (SDG 13), promoting life on land (SDG 15), and ensuring food secure nations (SDG 2). Figure 1(a) describes the study domain.Figure 1
Schematic configuration of the homogeneous porous media.
(a)
Study domain(b)
Porous pack with boundariesNaturally, a soil ridge resembles an experimental soil column (porous pack), heated at the top and cooled at the bottom so that it also contains a vapor zone near top surface, a liquid zone at the bottom, and a two-phase zone in between (Figure1(b), adapted from [77]). Hence, the model under consideration is that of an unsteady natural convection flow in a porous medium of nonbare ridge such that equations describing transport of moisture and solute concentration are coupled with sink terms. Further, it is assumed that buoyancy forces are present due to temperature and concentration variation of the fluid flow. The continuity equation (1) and equations governing flow, momentum, energy, moisture, and concentration are (2)–(6) [78]
(1)∂u∂x+∂w∂z=0,(2)ρ∂u∂t+u∂u∂x+w∂u∂z=μ∂2u∂x2+∂2u∂z2−μku−ρFu2k,(3)ρ∂w∂t+u∂w∂x+w∂w∂z=−∂P∂z+μ∂2w∂z2+∂2w∂z2−μkw−ρFw2k+ρgβT−T0+β∗C−C0,(4)C^A∂T∂t+C^wqw+C^vqv∂T∂x+C^wqw+C^vqv∂T∂z=λ^∂2T∂x2+∂2T∂z2−L0∂qv∂x+∂qv∂z−μ∂qw∂x2+ςT−T0,(5)∂θ∂t−∂∂xDθ∂θ∂x+DT∂T∂x−∂∂zDθ∂θ∂z+DT∂T∂z=∂Kθ∂z−Sz,t,(6)ρbκd+θw∂C∂t+qx∂C∂x+qz∂C∂z+ηw′θwC+ρbCa−ζsCs=∂∂xθwD∂C∂x+DTM∂2T∂x+∂∂zθwD∂C∂z+DTM∂2T∂z,where u and w are velocity components, and P is the pressure of the fluid; ρ and μ are density and viscosity of the fluid, respectively; and F and k are, respectively, the inertia resistance coefficient and the permeability coefficient. Following the Darcy-Boussinesq approximation assumption [79] on the vertical flow, g is the force due to gravity, and β and β∗ are, respectively, volumetric coefficients of expansion of heat and contaminants. The terms T0 and C0 in equation (3) represent medium surface temperature and concentration, respectively. The presence of sink terms, Sz,t in equation (5) and Cs in equation (6), mimics the plant-root uptake activity due to the nonbare ridge surface assumption.In equation (4), C^A=C^s+L0∂θv/∂t is the apparent heat capacity or the heat flow retardation factor [46], where C^s=ρcps, C^w=ρcpl, and C^v=ρcpv are, respectively, the volumetric heat capacities of the porous medium, the liquid phase, and the vapor phase; qw and qv are the liquid flux and vapor flux densities, respectively; λ^ is the apparent thermal conductivity of soil; and L0 is the volumetric latent heat of vaporization of liquid water. Otherwise, in equation (5), Dθ, DT, and Kθ are, respectively, unsaturated moisture diffusivity, coefficient of Soret effect of temperature in moisture flux, and unsaturated hydraulic conductivity.In equation (6), θw is volumetric water content −; D is the hydrodynamic dispersion coefficient m2s−1; q is the Darcian fluid flux density m3m−2s−1; ρb is the soil bulk density Kgm−3; Ca is the concentration of contaminants adsorbed onto the solid surface mgKg−1; ηw′ is the first-order reaction rate coefficient s−1; and Sc is the contaminant sink term Mgm−3s−1. Additionally, kd is the distribution coefficient for linear adsorption process m3Mg−1, ζs is the volumetric flux of water per unit volume of soil matrix s−1; Cs is the concentration of contaminants at sink; Rm=ρbκd+θw is the soil matrix retardation factor −; DTM is the thermal gradient diffusion coefficient m2s∘C−1; and qx and qz are, respectively, seepage velocity components in the horizontal and vertical directions.The initial and boundary conditions to (2)–(6) are [78]
(7a)u⟶x,z,0=0,Tx,z,0=Ti,θx,z,t=θ0Cx,z,0=Ci,(7b)ua,z,t=Uwm,wa,z,t=0,ub,z,t=Uwm,wb,z,t=0,(7c)ux,c,t=0,wx,c,t=0,ux,d,t=0,wx,d,t=Vwm,(7d)∂T∂xa=q′,∂T∂xb=q′,Tx,c,t=0,Tx,d,t=T0−Ti,(7e)∂θ∂xa=Ex,∂θ∂xb=Ex,θx,c,t=0,θx,d,t=θr,(7f)∂C∂xa=0,∂C∂xb=0,Cx,c,t=0,Cx,d,t=0.Applying the following nondimensional variables(8)x∗=xL,z∗=zL,ϕ=dL,u∗=uU0ϕ,w∗=wU0ϕ,q∗=qU0ϕ,P∗=PρU02,(9)T∗=T−TiT0−Ti,C∗=C−CiC0−Ci,t∗=tL/U0,we obtain (upon dropping the asterisk) the dimensionless equations [78]
(10)∂u∂t+ϕu∂u∂x+ϕw∂u∂z=1Re∂2u∂x2+∂2u∂z2−uRe·Λ−F∗Reu2,(11)∂w∂t+ϕu∂w∂x+ϕw∂w∂z=−1ϕ∂P∂z+1Re∂2w∂x2+∂2w∂z2−wRe∙Λ−F∗Rew2+1ϕGrT+GrCRe2,(12)∂T∂t+ϱ1ϕqx∂T∂x+qz∂T∂z=1Pr∙Re∂2T∂x2+∂2T∂z2−ϕEcReϱ2L0∂qv∂x+∂qv∂z+ϕ2EcRe∂qw∂x2+δ1Pr∙ReT,(13)Rm∂C∂t+ϕqx∂C∂x+qz∂C∂z=θwPe∂2C∂x2+∂2C∂z2+ζ¯∂2T∂x2+∂2T∂z2+ϱ3C+ϱ4Ca+ϱ5,where Pr=μcp/k is the Prandtl number; Re=ρU0L/μ is the Reynolds number; Ec=U02/cpT0−Ti is the Eckert number; F∗=ϕReF/Λ is the Forchheimer number, where Λ is the Darcy number; GrT=gβT0−TiL3ρ2/μ2 is the Grashof number due to soil temperature; GrC=gβC0−CiL3ρ2/μ2 is the Grashof number due to soil concentration; and Pe=U0L/D is the Peclet number.In addition,ρ1=ρcpw/ρbcps is the heat capacity ratio parameter, δ1=ς/λ^ is the heat generation/absorption parameter, ζ¯=DTMT0−Ti/U0L is the thermophoresis parameter, ρ2=L/U0 is the characteristic time, ρ3=ηw′θwL/U0 is the chemical reaction parameter, ρ4=ηw′ρbL/U0C0−Ci is the sorption parameter, ρ5=ζsLCs/U0C0−Ci is the chemical sink parameter, and Rm=1+ρbκd/θw is the soil matrix retardation factor.Equation (5) is highly nonlinear and degenerate because of Kθ and Dθ that also depend on the field variable θ. For simplicity, we first rewrite equation (5) as
(14)∂θ∂t=DT∂2T∂x2+∂2T∂z2+Dθ∂2θ∂x2+∂2θ∂z2−∂Kθ∂z+∂θz∂t,where ∂θz/∂t replaces the sink term Sz,t to describe the “rate of moisture uptake by the plant root in the z-direction” only. Then, using the relation [20]
(15)∂Kθ∂z=Ksθs−θr∂θ∂z,equation (14) is rewritten as
(16)∂θ∂t=DT∂2T∂x2+∂2T∂z2+Dθ∂2θ∂x2+∂2θ∂z2−Ksθs−θr∂θ∂z+∂θz∂t.Using the expression [34]
(17)θ∗=θ−θrθs−θr,where θ∗ is the “relative water/moisture content,” and the following nondimensional variables
(18)x∗=xd,z∗=zL,T∗=T−TiT0−Ti,t∗=DθtL2,lead to a simplified and nondimensional moisture (normalized) equation [78]
(19)Rn∂θ∂t=ΔTΔθLe′∂2T∂x2+∂2T∂z2+∂2θ∂x2+∂2θ∂z2−δ2∂θ∂z,where Rn=1+∂θz/∂u¯ is the moisture retardation factor; Le′=DT/Dθ is the Lewis number or “reciprocal” Luikov number; δ2=KsL/Dθθs−θr is the ratio of “apparent” diffusion coefficients; and ΔT=T0−Ti.The dimensionless initial and boundary conditions are [78]
(20a)u⟶x,z,0=0,Tx,z,0=Ti,θx,z,t=θ0,Cx,z,0=Ci,(20b)ua,z,t=UwmU0ϕ,wa,z,t=0,ub,z,t=UwmU0ϕ,wb,z,t=0,(20c)ux,c,t=0,wx,c,t=0,ux,d,t=0,wx,d,t=VwmU02ϕ,(20d)∂T∂xa=q′dT0−Ti,∂T∂xb=q′dT0−Ti,Tx,c,t=−TiT0−Ti,Tx,d,t=T0−2TiT0−Ti,(20e)∂C∂xa=0,∂C∂xb=0,Cx,c,t=0,Cx,d,t=1,(20f)∂θ∂xa=Exdθs−θr,∂θ∂xb=Exdθs−θr,θx,c,t=0,θx,d,t=θ0−θrθs−θr.Since it is generally difficult to obtain estimates of Darcy flux densities for the investigated domain, the soil heat (12) and contaminant (13) equations are instead coupled to the flow equations ((10) and (11)). To achieve this, we first consider the expression of Darcy flux density q⟶ [80]
(21)q⟶=−Khϕedhdl,where h is the hydraulic head m, ϕe is the effective porosity, and l is the single direction distance over which the measured change in the hydraulic head is measured m. Using the definitions of moisture specific capacity Ch=dθ/dh and unsaturated moisture diffusivity Dθ=Kh/Ch yields a new expression
(22)Dθ=Khdhdθ.Combining equations (21) and (22) yields
(23)q⟶=−Dθϕedθdl,such that
(24)∂q⟶∂l=−Dθϕe∂2θ∂l2.Using the relationship between the Darcy flux densityq⟶ and the average seepage velocity u⟶ [80],
(25)u⟶=q⟶ϕe,together with equation (24), we deduce the following expressions;
(26)qx=uϕe;qz=wϕe;∂q∂x=−Dθϕe∂2θ∂x2;∂q∂z=−Dθϕe∂2θ∂z2.Then, using the first two terms of (26), we replace the Darcy flux density components qx and qz in equations (12) and (13) with the “average” seepage velocity components u and w, respectively; the solutions of equations (10) and (11), respectively. Consequently, equations (12) and (13) are coupled with equations (10) and (11). In a similar manner, the last two expressions of (26) replace the term ∂q/∂x+∂q/∂z in equation (12) to be eloquently coupled with equation (19).Lastly, the parametersϱ3=ηw′θwL/U0, ϱ4=ηw′ρbL/U0C0−Ci, and ϱ5=ζsLCs/U0C0−Ci in equation (13) are further rearranged. Introducing a “one,” d/d to the RHS of ϱ3, results into
(27)ϱ3´=ηw′θwd2D˘,where D˘≡U0d is the apparent diffusivity of the porous matrix. Hence, equation (27) is expressed in the form of the Damko¨hler (Da) number [23]
(28)Da=kjb22D∗,when kj≡ηw′, b2≡d, and D∗≡D˘, where b, kj, and D are, respectively, half-fracture aperture length, reaction rate constant, and matrix diffusion coefficient. Hence, the nondimensional parameter ϱ′3 (equation (27)) is termed “apparent” Damko¨hler number. In a similar manner, the rearranged forms of ϱ4 and ϱ5 are, respectively,
(29)ϱ4´=ρbκdC0−CiDaandϱ5´=ζsd2CsD˘C0−Ci,where ϱ4´ in (29) is now called the “weighted” Damko¨hler number.


## 3. Derivation of Solution Method

### 3.1. MV-SLQLM Approximation and Construction in (2+1) Dimensions

Suppose we have a general nonlinear PDE of the form(30)∂∂t=F∂2u∂x2,∂2u∂z2,∂u∂x,∂u∂z,u,for x,z,t∈a,b×c,d×0,T, where F is a nonlinear operator.To solve equation (30), first, we require the QLM to linearize any nonlinear equations [67]. Then, we approximate the unknown solution ux,z,t of each linear equation using the triple Lagrange interpolation functions Lix, Ljz, and Lkt at the Chebyshev-Gauss-Lobatto (CGL) collocation points x~,z~,t~ located in the local support domain −1,13. Like in the generalized differential quadrature (GDQ) method [81], the use of Gauss-Lobatto (GL) rule permits grid points to be located on the boundaries, appropriate for the application of boundary conditions. Moreover, grid points with nonuniform intervals provide more accurate results in the calculation of numerical integrals and derivatives.As with other related meshless methods [49, 82], accuracy and efficiency of any spectral method depend on the number of nodes (collocation points) in the local support domain. We thus obtain the spectral domain −1,13 from the original domain x,z,t∈a,b×c,d×0,T via the linear transformation,
(31)x=12b−ax~+12a+b,z=12d−cz~+12c+d,andt=12t~+1,where [55]
(32)x~i=cosπiNxi=0Nx,z~j=cosπjNzj=0Nz,t~k=cosπkNtk=0Nt,is the collocation points described by the Chebyshev polynomial TNv of order N∈ℕ,
(33)TNv=cosNarccosv.Then, we approximate each functionux,z,t by the Chebyshev interpolation uNx,z,t, where
(34)UNx,z,t=∑k=0Nt∑j=0Nz∑i=0NxukjiLktLjzLix,together with the corresponding approximate derivatives at the collocation points,
(35a)∂n∂xnux~i,z~j,t~k=∑r=0Nt∑q=0Nz∑p=0NxukjiLrt~kLqz~jDpin,(35b)∂n∂znux~i,z~j,t~k=∑r=0Nt∑q=0Nz∑p=0NxukjiLrt~kDqjnLpx~i,(35c)∂n∂tnux~i,z~j,t~k=∑r=0Nt∑q=0Nz∑p=0NxukjiDrknLqz~jLpx~i,where Dpin=Linx~i, Dqjn=Ljnz~j, and Drkn=Lknt~k are, respectively, Nx+1×Nx+1, Nz+1×Nz+1, and Nt+1×Nt+1 Chebyshev differentiation matrices for i,j,k=0,1,⋯,N.Further, a Lagrange polynomial of orderN, for example, at the Chebyshev-Gauss-Lobatto points Liv, is defined as [72]
(36)Liv=−1i+11−v2TN′vσ¯iN2v−vi,i=0,1,2,⋯,N,where σ¯i=1 for i=1,2,⋯,N−1, and σ¯0=σ¯N=2. In addition, the first order Chebyshev differentiation matrix at the collocation points is given by [52, 55]
(37)Dij1=σ¯j−1j+iσ¯ivj−vi,i≠j,−vi21−vi2,i=j≠0,N,2N2+16,i,j=0,2N2+16,i,j=N.For simplicity and efficiency, we express equations (34), (35a), (35b), and (35c) in matrix form, using the Kronecker product, as [72]
(38a)ux~i,z~j,t~k=Lt~k⊗Lz~j⊗Lx~iU,(38b)∂n∂xnux~i,z~j,t~k=Lt~k⊗Lz~j⊗D⟶xnU,(38c)∂n∂znux~i,z~j,t~k=Lt~k⊗D⟶zn⊗Lx~iU,(38d)∂n∂tnux~i,z~j,t~k=D⟶tn⊗Lz~j⊗Lx~iU,where U=u000,⋯,u00Nx,u010,⋯,u01Nx,⋯,u0NxNz,u100,⋯,uNtNzNxT,
(39a)Lt~k=L0t~k,L1t~k,L2t~k,⋯,LNtt~k,(39b)Lz~j=L0z~j,L1z~j,L2z~j,⋯,LNzz~j,(39c)Lx~i=L0x~i,L1x~i,L2x~i,⋯,LNxx~i,and
(40)Lt~k⊗Lz~j⊗Lx~i=L0t~kL0z~jLx~i,⋯,L0t~kLNzz~jLx~i,L1t~kL0z~jLx~i,⋯,L1t~kLNzz~jLx~i,⋯,LNtt~kLNzz~jLx~i.Third, we formulate the spectral collocated coupled system, defined for each equation, as(41a)U~s+1=Us,Ws,Uxs,Uxxs,Uzs,Uzzs,Ws+1,(41b)W~s+1=Ws,Us,Wxs,Wxxs,Wzs,Wzzs,Us+1,(41c)T~s+1=Ts,Us,Ws,Txs,Txxs,Tzs,Tzzs,Θs+1,(41d)Θ~s+1=Θs,Θxs,Θxxs,Θzs,Θzzs,Ts+1,(41e)C~s+1=Cs,Us,Ws,Cxs,Cxxs,Czs,Czzs,Ts+1.We apply the LLM algorithm to decouple the spectral collocated system (41a), (41b), (41c), (41d), and (41e) and solve each equation as follows:
(a)
Solve forU~s+1 in equation (41a), assuming Ws and Ws+1 are known from the previous iteration(b)
WithU~s+1 now known, solve for W~s+1 in equation (41b)(c)
Solve forT~s+1 in (41c), using U~s+1 and W~s+1, assuming Θs+1 is known from the previous iteration(d)
WithT~s+1 now known, proceed to solve for Θ~s+1 in (41d)(e)
Finally, usingU~s+1, W~s+1, and T~s+1, solve for C~s+1 in (41e)Beginning with the initial guesses, steps (a) to (e) are repeated, in turn, fors=0,1,2,⋯ to derive approximate solutions to the given model.


### 3.2. Error, Stability, and Convergence Analysis

#### 3.2.1. Error Bound

Following the approach in Bhrawy [83] and Samuel and Motsa [84], in this subsection, we present an analytical expression for the error norm and an upper bound on the expected error of the approximation, equation (34), to the smooth function ux,z,t∈Ω≡a,b×c,d×0,T. Let
(42)∏=SpanLixLjzLkt,∀i=0,1,2,⋯,Nx;j=0,1,2,⋯,Nz;k=0,1,2,⋯,Nt.Given a smooth functionux,z,t∈Ω≡a,b×c,d×0,T. Let UNx,z,t, defined in (34), be the best approximation of ux,z,t so that ∀VN∈∏,
(43)ux,z,t−UNx,z,t∞≤ux,z,t−VNx,z,t∞.Then, there exists arbitrary pointsγ1,γ1′∈a,b, γ2,γ2′∈c,d, and γ3,γ3′∈0,T, so that
(44)ux,z,t−VNx,z,t=∂Nx+1uγ1,z,t∂xNx+1Nx+1!∏i=0Nxx−xi+∂Nz+1ux,γ2,t∂zNz+1Nz+1!∏j=0Nzz−zj+∂Nt+1ux,z,γ3∂tNt+1Nt+1!∏k=0Ntt−tk−∂Nx+Nz+Nt+3uγ1,γ2,γ3∂xNx+1∂zNz+1∂tNt+1Nt+1!Nz+1!Nx+1!×∏i=0Nxx−xi×∏j=0Nzz−zj×∏k=0Ntt−tk,where VNx,z,t is the multivariate chebyshev interpolating polynomial of ux,z,t at xi, zj, and tk.Taking the infinite norm of equation (44) yields
(45)ux,z,t−VNx,z,t∞≤maxx,z,t∈Ω∂Nx+1uγ1,z,t∂xNx+1∏i=0Nxx−xiNx+1!+maxx,z,t∈Ω∂Nz+1ux,γ2,t∂zNz+1∏j=0Nxz−zjNz+1!+maxx,z,t∈Ω∂Nt+1ux,z,γ3∂tNt+1∏k=0Ntt−tkNt+1!−maxx,z,t∈Ω∂Nx+Nz+Nt+3uγ1′,γ2′,γ3′∂xNx+1∂zNz+1∂tNt+1×∏i=0Nxx−xi×∏j=0Nxz−zj×∏k=0Ntt−tkNx+1!Nz+1!Nt+1!.Sinceux,z,t is a smooth function on Ω; then, its Nx+1, Nz+1, Nt+1, and Nx+Nz+Nt+3-th partial derivatives are both continuous and upper-bounded, so that
(46)maxx,z,t∈Ω∂Nx+1uγ1,z,t∂xNx+1≤K1,maxx,z,t∈Ω∂Nz+1ux,γ2,t∂zNz+1≤K2,maxx,z,t∈Ω∂Nt+1ux,z,γ3∂tNt+1≤K3,maxx,z,t∈Ω∂Nx+Nz+Nt+3uγ1′,γ2′,γ3′∂xNx+1∂zNz+1∂tNt+1≤K4,where K1,K2,K3,andK4 are arbitrary constants in ℝ+.Substituting (46) into (45) yields the norm error bound estimate
(47)ux,z,t−VNx,z,t∞≤K1∏i=0Nxx−xi∞Nx+1!+K2∏j=0Nzz−zj∞Nz+1!+K3∏k=0Ntt−tk∞Nt+1!−K4×∏i=0Nxx−xi∞×∏j=0Nzz−zj∞×∏k=0Ntt−tk∞Nt+1!Nz+1!Nx+1!,where VN∈∏Nx,Nz,Nt.However, (47) can be further simplified. From the linear relationships defined in (31), we have
(48)x−xi=12b−ax~−x~i,z−zj=12d−cz~−z~j,andt−tk=12t~−t~k,so that from (47),
(49)∏i=0Nxx−xi∞=∏i=0Nxb−a2x~−x~i∞≤b−a2Nxmaxx~,xi∈−1,1∏i=0Nxx~−x~i≤b−aNx+12Nx+1maxx~,xi∈−1,1∏i=0Nxx~−x~i.Using equation (36),
(50)Lix~=−1i+11−x~2TNxx~′σ¯iNx2x~−x~i,i=0,1,2,⋯,Nx,so that the Nx+1-th degree Lagrange polynomial is
(51)LNx+1x~=−1Nx+21−x~2TNxx~′σ¯iNx2x~−x~i.Taking the modulus of equation (51) and rearranging the terms yield
(52)x~−x~i=1−x~2TNxx~′Nx2LNx+1x~.Since the Chebyshev interpolating polynomial, equation (33) is the solution to the Chebyshev differential equation [85]
(53)1−x2d2ydx2−xdydx+n2y=0,for nonnegative integers n,n∈ℤ0 and x<1, it follows that
(54)1−x2TNxx′=−xTNxx−n2.Using equation (54) and applying the triangular inequality, equation (52) can then be expressed as
(55)x~−x~i≤−x~TNxx~+−n2Nx2LNx+1x~.Employing the inequality relation stated in [84] and noting that TNxx~≤1,∀xi∈−1,1, then
(56)1−x~2TNx′x~=−x~TNxx~+−n2≤2Nx,so that (55) becomes
(57)x~−x~i≤2Nx.Therefore, from (57), maxx~,xi∈−1,1∏i=0Nxx~−x~i≤2Nx+2/Nx. It then also follows from (49) that
(58)∏i=0Nxx−xi∞≤b−aNx+12Nx+1maxx~,xi∈−1,1∏i=0Nxx~−x~i≤2Nx+2b−aNx+1Nx2Nx+1.In a similar manner,(59)∏j=0Nzz−zj∞≤d−cNz+12Nz+1maxz~,zj∈−1,1∏j=0Nzz~−z~i≤2Nz+2d−cNz+1Nz2Nz+1,(60)∏k=0Ntt−tk∞≤12Nt+1maxt~,tk∈−1,1∏k=0Ntt~−t~k.Substituting (58), (59), and (60) into (47) yields
(61)ux,z,t−VNx,z,t∞≤K12Nx+2b−aNx+1Nx2Nx+1Nx+1!+K22Nz+2d−cNz+1Nz2Nz+1Nz+1!+K32Nt+2Nt2Nt+1Nt+1!−K4×2Nx+2b−aNx+1×2Nz+2d−cNz+1×2Nt+2NxNzNt2Nx+Nz+Nt+3Nt+1!Nz+1!Nx+1!.Hence, by combining (61) and (43), an upper bound of the error for the approximate solutions is obtained. Since the model, equations (10)–(13) and (19), does not have the exact solution, it is inferred that the convergence of the proposed method depends basically on the above error bound.


#### 3.2.2. Stability

According to Noborio et al. [30], the display of dispersions or oscillatory errors in the numerical solutions is completely unavoidable regardless of whether one uses FD or FE methods to simulate transport due to the characteristics of the solute transport equation. Hence, to ensure a generally stable MV-SLQLM, two complementary measures are employed on the basis of solute transport equation, equation (13). First, in reference [86], we investigated the existence and uniqueness of the classical solution to equation (13) using the proposed “VF-EM” technique. Defining an error e~≔ct−c~t and making inference of Ilati and Dehghan [82], the result (see equation 3.10 in [86]) means e~ is strongly stable. Moreover, the well-posed analysis of equation (13) in time integral also demonstrates that the numerical technique is stable in L2-norm. Second, we examine numerically the stability of a discretized scheme for convection-diffusion terms in equation (13). Among the existing numerical methods for analyzing scheme’s stability on the basis of convection-diffusion terms, Yu et al. [87] observed that majority are usually based on five assumptions, i.e., one-dimensional, linear, first kind boundary condition, free source term, and uniform grid system.However, from [87], any deviation from the aforementioned five assumptions may attribute to enhanced convective stability. In this study, except for the linear nature, equation (13) is of (2+1) dimensions, together with a nonfree source term, and boundary conditions of first and second types. It is worth also mentioning that the convective terms of equations (12) and (13) use flow velocity functions derived from the pair of highly nonlinear momentum equations, equations (10) and (11). Moreover, the MV-SLQLM is built on Lagrange interpolation bases that do not depend on the order of nodes (nonuniform) and are optimally stable [88, 89]. In spite of the fact that these properties naturally contribute to a more stable numerical scheme, we further check the stability of the convective-diffusion solution profiles following the study of Abdullah [90]. Using a one-dimension steady convective-diffusion equation, Abdullah [90] described a linear relationship between grid size and local Pe number;
(62)Nmin=mPe+b,3.125<Pe≤100,3,0≤Pe≤3.125,where Nmin is the minimum grid size, sufficient for the flow profile to behave physically correctly; Pe is the grid or local Peclet number; and m and b are the slope and intercept constants, respectively. Since we have a low Pe, 0≤Pe≤100, advection-diffusion-reaction equation (13), we investigate the effect of varying Pe on the solute concentration profiles at the minimum grid number, Nmin.


#### 3.2.3. Convergence

Making inference of (61) and (43), convergence to the expected solution is affirmed by considering the L∞-norm of the difference between successive iterations, using the residual error analysis
(63)EU~=max0≤i≤NtU~s+1,i−U~s,i∞,EW~=max0≤i≤NtW~s+1,i−W~s,i∞,ET~=max0≤i≤NtT~s+1,i−T~s,i∞,EΘ~=max0≤i≤NtΘ~s+1,i−Θ~s,i∞,EC~=max0≤i≤NtC~s+1,i−C~s,i∞.,defined for the difference between values of successive iterations.


### 3.3. Numerical Application

Now, we apply the above described procedure, the MV-SLQLM, to the model under consideration, equations (10)–(13) and (19), together with (20a), (20b), (20c), (20d), (20e), and (20f). Beginning with the nonlinear Darcy-Brinkman-Forchheimer equation (equation (10)) in the x−direction
(64)∂u∂t+ϕuux+ϕwuz=a1uxx+uzz−a2u−a3u2,where a1=1/Re, a2=1/Re·Λ, and a3=F∗/Re, we employ the QLM to look for the linear expression
(65)α2,sx,z,tus"+1+α1,sx,z,tus′+1+α0,sx,z,tus+1−uts+1=Rsx,z,t,where
(66a)Rsx,z,t=∑m=02αm,sx,z,tusm−Fus,us′,us",(66b)Fus,us′,us"=a1uxxs+uzzs−ϕuux−ϕwuz−a2us−a3us2,(66c)αm,sx,z,t=∂Fus,us′,us"∂usm,m=0,1,2.Using equations (66b) and (66c), we have
(67)α0,sx,z,t=−ϕuxs−a2−2a3us;α1,sx,z,t=−ϕus−ϕws;α2,sx,z,t=a1,so that(68)Rsx,z,t=−ϕuxsus+ϕwsuzs+ϕws+1uzs−a3us2.Hence, equation (65) becomes
(69)a1uxxs+1+uzzs+1−ϕusuxs+1−ϕwsuzs+1−ϕuxs+a2+2a3usus+1−uts+1−ϕuzsws+1=−ϕuxsus+wsuzs−a3us2.Applying the Chebyshev spectral collocation method to equation (69) yields
(70)a1I⟶Nt+1⊗I⟶Nz+1⊗1hx2D⟶x2+I⟶Nt+1⊗1hz2D⟶z2⊗I⟶Nx+1−ϕdiagUsI⟶Nt+1⊗I⟶Nz+1⊗1hxD⟶x−ϕdiagWsI⟶Nt+1⊗1hzD⟶z⊗I⟶Nx+1−diagϕUxs+a2+2a3Us−1htD⟶t⊗I⟶Nz+1⊗I⟶Nx+1Us+1−ϕdiagUzsWs+1=−ϕUxsUs+UzsWs−a3UsUs,where hx=xl−x0/2, hz=zl−z0/2, and ht=T/2; and I⟶Nt+1, I⟶Nx+1, and I⟶Nz+1 are, respectively, identity matrices of sizes Nt+1×Nt+1, Nx+1×Nx+1, and Nz+1×Nz+1.SinceU in equation (70) is known at the boundaries of the domain, rewriting equation (70) for the interior points of the domain, yields the new expression
(71)a1I⟶Nt⊗I⟶Nz−1⊗1hx2D⟶x22:Nx,:+I⟶Nt⊗1hz2D⟶z22:Nz,:⊗I⟶Nx−1−diagϕUxs+a2+2a3Us−ϕdiagUsI⟶Nt⊗I⟶Nz−1⊗1hxD⟶x2:Nx,:−ϕdiagWsI⟶Nt⊗1hzD⟶z2:Nz,:⊗I⟶Nx−1−1htD⟶t1:Nt,:⊗I⟶Nz−1⊗I⟶Nx−1U~s+1−ϕdiagUzsWs+1=−ϕUxsUs+UzsWs−a3U~sU~s−ψijk1+ψ~ijk1−σijk1+σ~ijk1+δijk1,where U~=U2:Nx,2:Nz,1:Nt, and
(72a)ψ1=a1D⟶x22:Nx,1⊗U1,2:Nz,1:Nt−UsD⟶x2:Nx,1⊗U1,2:Nz,1:Nt,(72b)ψ~ijk1=a1D⟶x22:Nx,Nx+1⊗UNx+1,2:Nz,1:Nt−UsD⟶x2:Nx,Nx+1⊗UNx+1,2:Nz,1:Nt,(72c)σijk1=a1D⟶z22:Nz,1⊗U2:Nx,1,1:Nt−WsD⟶z2:Nz,1⊗U2:Nx,1,1:Nt,(72d)σ~ijk1=a1D⟶z22:Nz,Nz+1⊗U2:Nx,Nz+1,1:Nt−WsD⟶z2:Nz,1⊗U2:Nx,1,1:Nt,(72e)δijk1=−D⟶t1:Nt,Nt+1⊗U2:Nx,2:Nz,Nt+1.Thus, at step (a) of the LLM, we need to solve forU~s+1 in the matrix equation (equation (71)) defined as
(73)A⟶1U~s+1=B⟶1+K⟶1,where
(74)A⟶1=a1I⟶Nt⊗I⟶Nz−1⊗1hx2D⟶x22:Nx,:+I⟶Nt⊗1hz2D⟶z22:Nz,:⊗I⟶Nx−1−diagϕUxs+a2+2a3Us−ϕdiagUsI⟶Nt⊗I⟶Nz−1⊗1hxD⟶x2:Nx,:−ϕdiagWsI⟶Nt⊗1hzD⟶z2:Nz,:⊗I⟶Nx−1−1htD⟶t1:Nt,:⊗I⟶Nz−1⊗I⟶Nx−1,and
(75a)B⟶1=−ϕUxsUs+UzsWs−a3U~sU~s+ϕdiagUzsWs+1,(75b)K⟶1=−ψijk1+ψ~ijk1−σijk1+σ~ijk1+δijk1.In a similar manner, the Darcy-Brinkman-Forchheimer equation (equation (11)) in z−direction is expressed as
(76)∂w∂t=−b0Pz−ϕuwx−ϕwwz+b1wxx+wzz−b2w−b3w2+b4,where b0=1/ϕ, b1=1/Re, b2=1/Re·Λ, b3=F∗/Re, and b4=1/ϕGrT+GrC/Re2.Applying the QLM, we obtain(77)b1wxxs+1+wzzs+1−ϕuswxs+1−ϕwswzs+1−ϕwzs+b2+2b3wsws+1−b0Pzs+1−wts+1−ϕwxsus+1=−ϕwxus+wzsws−b3wsws−b4.Then, by applying the Chebyshev spectral collocation method to equation (77), we get
(78)b1I⟶Nt+1⊗I⟶Nz+1⊗1hx2D⟶x2+I⟶Nt+1⊗1hz2D⟶z2⊗I⟶Nx+1−ϕdiagUsI⟶Nt+1⊗I⟶Nz+1⊗1hxD⟶x−ϕdiagWsI⟶Nt+1⊗1hzD⟶z⊗I⟶Nx+1−diagϕWzs+b2+2b3Ws−1htD⟶t⊗I⟶Nz+1⊗I⟶Nx+1Ws+1−b0Pzs+1−ϕdiagWxsUs+1=−ϕWxsUs+WzsWs−b3WsWs−b4,where hx, hz, ht, I⟶Nt+1, I⟶Nx+1, and I⟶Nz+1 are all defined as before.Considering the known values ofW at the initial and boundary conditions, equation (78) can be rearranged as
(79)b1I⟶Nt⊗I⟶Nz−1⊗1hx2D⟶x22:Nx,:+I⟶Nt⊗1hz2D⟶z22:Nz,:⊗I⟶Nx−1−ϕdiagUsI⟶Nt⊗I⟶Nz−1⊗1hxD⟶x2:Nx,:−ϕdiagWsI⟶Nt⊗1hzD⟶z2:Nz,:⊗I⟶Nx−1−diagϕWzs+b2+2b3Ws−1htD⟶t1:Nt,:⊗I⟶Nz−1⊗I⟶Nx−1W~s+1=b0P~zs+1+ϕdiagWxsU~s+1−ϕW~xsU~s+W~zsW~s−b3W~sW~s−ψijk2+ψ~ijk2−σijk2+σ~ijk2+δijk2−b4,where W~=W2:Nx,2:Nz,1:Nt, and
(80a)ψ2=b1D⟶x22:Nx,1⊗W1,2:Nz,1:Nt−UsD⟶x2:Nx,1⊗W1,2:Nz,1:Nt,(80b)ψ~ijk2=b1D⟶x22:Nx,Nx+1⊗WNx+1,2:Nz,1:Nt−UsD⟶x2:Nx,Nx+1⊗WNx+1,2:Nz,1:Nt,(80c)σijk2=b1D⟶z22:Nz,Nz+1⊗U2:Nx,Nz+1,1:Nt−WsD⟶z2:Nz,1⊗W2:Nx,1,1:Nt,(80d)σ~ijk2=b1D⟶z22:Nz,Nz+1⊗U2:Nx,Nz+1,1:Nt−WsD⟶z2:Nz,Nz+1⊗W2:Nx,Nz+1,1:Nt,(80e)δijk2=−D⟶t1:Nt,Nt+1⊗W2:Nx,2:Nz,Nt+1.Thus, in step (b) of the LLM, we solve forW~s+1 in the matrix equation
(81)A⟶2W~s+1=B⟶2+K⟶2,where
(82)A⟶2=b1I⟶Nt⊗I⟶Nz−1⊗1hx2D⟶x22:Nx,:+I⟶Nt⊗1hz2D⟶z22:Nz,:⊗I⟶Nx−1−ϕdiagUsI⟶Nt⊗I⟶Nz−1⊗1hxD⟶x2:Nx,:−ϕdiagWsI⟶Nt⊗1hzD⟶z2:Nz,:⊗I⟶Nx−1−diagϕWzs+b2+2b3Ws−1htD⟶t1:Nt,:⊗I⟶Nz−1⊗I⟶Nx−1,and
(83a)B⟶2=b0P~zs+1+ϕdiagWxsU~s+1−ϕW~xsU~s+W~zsW~s−b3W~sW~s,(83b)K⟶2=−ψijk2+ψ~ijk2−σijk2+σ~ijk2+δijk2−b4.The pressure termP~zs+1 is iteratively estimated from an initial guess
(84)Px,z=x4−4π2sinh2πsinhπxcosπz,deduced from the analytical solution [91]
(85)Px,z=x4−4∑n=1,odd∞1nπ2sinh2nπsinhnπxcosnπz,of the pressure equation ∇2P=0, assuming an initial state of P=0 everywhere and the zero-flux top and bottom boundary conditions.Next, we express the heat equation as(86)∂T∂t=c¯1∂2T∂x2+∂2T∂z2−c¯2u∂T∂x+w∂T∂z+c¯3∂2θ∂x2+∂2θ∂z2+c¯4T+ϕ2EcReϕ∂u∂x2,where c¯1=1/PrRe, c¯2=ϱ1ϕ2, c¯3=ϕEcϱ2L0Dθ/Reϕe, c¯4=δ1/PrRe, and c¯5=ϕ2Ec/Re.In terms of the QLM, (86) can be expressed as
(87)c¯1Txxs+1+Tzzs+1−c¯2usTxs+1−c¯2wsTzs+1+c¯3θxxs+1+θzzs+1+c¯4Ts+1−Tts+1=−ϕ2c¯5uxs2.Applying the Chebyshev spectral collocation method to equation (87), we get
(88)c¯1I⟶Nt+1⊗I⟶Nz+1⊗1hx2D⟶x2+I⟶Nt+1⊗1hz2D⟶z2⊗I⟶Nx+1−diagc¯2UsI⟶Nt+1⊗I⟶Nz+1⊗1hxD⟶x−diagc¯2WsI⟶Nt+1⊗1hzD⟶z⊗I⟶Nx+1+diagc4−1htD⟶t⊗I⟶Nz+1⊗I⟶Nx+1Ts+1+c¯3I⟶Nt+1⊗I⟶Nz+1⊗1hx2D⟶x2+I⟶Nz+1⊗1hz2D⟶z2⊗I⟶Nx+1Θs+1=−ϕ2c¯5Uxs2,where hx, hz, ht, I⟶Nt+1, I⟶Nx+1, and I⟶Nz+1 are all as previously described.In a similar manner, with the known values at the boundaries and initial points, we solve(89)c¯1I⟶Nt⊗I⟶Nz−1⊗1hx2D⟶x22:Nx,:+I⟶Nt⊗1hz2D⟶z22:Nz,:⊗I⟶Nx−1−diagc¯2UsI⟶Nt⊗I⟶Nz−1⊗1hxD⟶x2:Nx,:−diagc¯2WsI⟶Nt⊗1hzD⟶z2:Nz,:⊗I⟶Nx−1+diagc4−1htD⟶t1:Nt,:⊗I⟶Nz−1⊗I⟶Nx−1T~s+1+c¯3Θ~xxs+1+Θ~zzs+1=−ϕ2c¯5Uxs2−ψijk3+ψ~ijk3−σijk3+σ~ijk3+δijk3,where T~=T2:Nx,2:Nz,1:Nt, U~s≈U~s+1, W~s≈W~s+1, Θ~=Θ2:Nx,2:Nz,1:Nt, and
(90a)ψ3=c¯1D⟶x22:Nx,1⊗T1,2:Nz,1:Nt−c¯2U~sD⟶x2:Nx,1⊗T1,2:Nz,1:Nt,(90b)ψ~ijk3=c¯1D⟶x22:Nx,Nx+1⊗TNx+1,2:Nz,1:Nt−c¯2U~sD⟶x2:Nx,Nx+1⊗TNx+1,2:Nz,1:Nt,(90c)σijk3=c¯1D⟶z22:Nz,1⊗T2:Nx,1,1:Nt−c¯2W~sD⟶z2:Nz,1⊗T2:Nx,1,1:Nt,(90d)σ~ijk3=c¯1D⟶z22:Nz,Nz+1⊗T2:Nz,Nz+1,1:Nt−c¯2W~sD⟶z2:Nz,Nz+1⊗T2:Nx,Nz+1,1:Nt,(90e)δijk3=−D⟶t1:Nt,Nt+1⊗T2:Nx,2:Nz,Nt+1.Henceforth, at step (c) of the LLM, we solve the matrix equation,(91)A⟶3T~s+1=B⟶3+K⟶3,where
(92)A⟶3=c¯1I⟶Nt⊗I⟶Nz−1⊗1hx2D⟶x22:Nx,:+I⟶Nt⊗1hz2D⟶z22:Nz,:⊗I⟶Nx−1−diagc¯2UsI⟶Nt⊗I⟶Nz−1⊗1hxD⟶x2:Nx,:−diagc¯2WsI⟶Nt⊗1hzD⟶z2:Nz,:⊗I⟶Nx−1+diagc41htD⟶t1:Nt,:⊗I⟶Nz−1⊗I⟶Nx−1,and,
(93a)K⟶3=−ψijk3+ψ~ijk3−σijk3+σ~ijk3+δijk3,(93b)B⟶3=−c¯3Θ~xxs+Θ~zzs−ϕ2c¯5U~xs2.Since solving the system (91) requires knowledge of θ~s+1, we then consider the moisture equation, which is also coupled with temperature T,
(94)Rn∂θ∂t=θxx+θzz−δ2θz−b¯Txx+Tzz,where b¯=LeT1−T0. Applying the QLM, we need to solve
(95)θxxs+1+θzzs+1−δ2θzs+1−b¯Txxs+Tzzs=Rnθts+1.Then, applying the Chebyshev spectral collocation method to equation (95), we obtain
(96)I⟶Nt+1⊗I⟶Nz+1⊗1hx2D⟶x2+I⟶Nt+1⊗1hz2D⟶z2⊗I⟶Nx+1−diagδ2I⟶Nt+1⊗1hzD⟶z⊗I⟶Nx+1−Rn1htD⟶t⊗I⟶Nz+1⊗I⟶Nx+1Θs+1−b¯I⟶Nt+1⊗I⟶Nz+1⊗1hx2D⟶x2+I⟶Nt+1⊗1hz2D⟶z2⊗I⟶Nx+1Ts+1=0,where hx, hz, ht, I⟶Nt+1, I⟶Nx+1, and I⟶Nz+1 are all as previously described.Taking all the known initial and boundary values ofθs+1 and Ts+1 to the right-hand side of equation (96), we now solve the equation
(97)I⟶Nt⊗I⟶Nz−1⊗1hx2D⟶x22:Nx,:+I⟶Nt⊗1hz2D⟶z22:Nz,:⊗I⟶Nx+1−diagδ2I⟶Nt⊗1hzD⟶z2:Nz,:⊗I⟶Nx−1−Rn1htD⟶t1:Nt,:⊗I⟶Nz−1⊗I⟶Nx−1Θ~s+1−b¯Txxs+1+Tzzs+1=−ψijk4+ψ~ijk4−σijk4+σ~ijk4+δijk4,where
(98a)ψ4=D⟶x22:Nx,1⊗Θ1,2:Nz,1:Nt,(98b)ψ~ijk4=D⟶x22:Nx,Nx+1⊗ΘNx+1,2:Nz,1:Nt,(98c)σijk4=D⟶z22:Nz,1⊗Θ2:Nx,1,1:Nt−δ2D⟶z2:Nz,1⊗θ2:Nx,1,1:Nt,(98d)σ~ijk4=D⟶z22:Nz,Nz+1⊗Θ2:Nz,Nz+1,1:Nt−δ2D⟶z2:Nz,Nz+1⊗Θ2:Nx,Nz+1,1:Nt,(98e)δijk4=−D⟶t1:Nt,Nt+1⊗θ2:Nx,2:Nz,Nt+1,so that, at step (d) of the LLM, we solve the matrix equation
(99)A⟶4Θ~s+1=B⟶4+K⟶4,where
(100)A⟶4=I⟶Nt⊗I⟶Nz−1⊗1hx2D⟶x22:Nx,:+I⟶Nt⊗1hz2D⟶z22:Nz,:⊗I⟶Nx−1−diagδ2I⟶Nt⊗1hzD⟶z2:Nz,:⊗I⟶Nx−1−Rn1htD⟶t1:Nt,:⊗I⟶Nz−1⊗I⟶Nx−1,and
(101a)B⟶4=b¯Txxs+Tzzs,(101b)K⟶4=−ψijk4+ψ~ijk4−σijk4+σ~ijk4+δijk4.Finally, applying the QLM, the solute equation(102)Rm∂C∂t+e¯1u∂C∂x+w∂C∂z=e¯1∂2C∂x2+∂2C∂z2+ζ¯∂2T∂x2+∂2T∂z2+ϱ3C+ϱ4Ca+ϱ5,where e¯1=ϕ2 and e¯2=θw/Pe become
(103)e¯2Cxxs+1+Czzs+1−e¯1usCxs+1−e¯1wsCzs+1+ϱ3+ϱ4Cs+1−Cts+1=−ζ¯Txxs+Tzzs−ϱ5.Now, applying the Chebyshev spectral collocation to equation (103), we obtain
(104)e¯2I⟶Nt+1⊗I⟶Nz+1⊗1hx2D⟶x2+I⟶Nt+1⊗1hz2D⟶z2⊗I⟶Nx+1−diage¯1UsI⟶Nt+1⊗I⟶Nz+1⊗1hxD⟶x−diage¯1WsI⟶Nt+1⊗1hzD⟶z⊗I⟶Nx−1+diagϱ3+ϱ4−1htD⟶t⊗I⟶Nz+1⊗I⟶Nx+1Cs+1=ζ¯I⟶Nt+1⊗I⟶Nz+1⊗1hx2D⟶x2+I⟶Nt+1⊗1hz2D⟶z2⊗I⟶Nx+1Ts−diagϱ5,where hx, hz, ht, I⟶Nt+1, I⟶Nx+1, and I⟶Nz+1 are all as previously described. Taking out known values of initial and boundary conditions, we obtain
(105)e¯2I⟶Nt⊗I⟶Nz−1⊗1hx2D⟶x22:Nx,:+I⟶Nt⊗1hz2D⟶z22:Nz,:⊗I⟶Nx−1−diage¯1UsI⟶Nt⊗I⟶Nz−1⊗1hxD⟶x2:Nx,:−diage¯1WsI⟶Nt⊗1hzD⟶z2:Nz,:⊗I⟶Nx−1+diagϱ3+ϱ4−1htD⟶t1:Nt,:⊗I⟶Nz−1⊗I⟶Nx−1C~s+1=−ζ¯T~xxs+T~zzs−ψijk5+ψ~ijk5−σijk5+σ~ijk5+δijk5−diagϱ5,where C~=C2:Nx,2:Nz,1:Nt, U~s≈U~s+1, W~s≈W~s+1, and
(106a)ψijk5=e¯2D⟶x22:Nx,1⊗C1,2:Nz,1:Nt−e¯1U~sD⟶x2:Nx,1⊗C1,2:Nz,1:Nt,(106b)ψ~ijk5=e¯2D⟶x22:Nx,Nx+1⊗CNx+1,2:Nz,1:Nt−e¯1U~sD⟶x2:Nx,Nx+1⊗CNx+1,2:Nz,1:Nt,(106c)σijk5=e¯2D⟶z22:Nz,1⊗C2:Nx,1,1:Nt−e¯1W~sD⟶z2:Nz,1⊗C2:Nx,1,1:Nt,(106d)σ~ijk5=e¯2D⟶z22:Nz,Nz+1⊗C2:Nz,Nz+1,1:Nt−e¯1W~sD⟶z2:Nz,Nz+1⊗C2:Nx,Nz+1,1:Nt,(106e)δijk5=−D⟶t1:Nt,Nt+1⊗C2:Nx,2:Nz,Nt+1.Thus, at step (e), we solve forC~s+1 in the matrix equation
(107)A⟶5C~s+1=B⟶5+K⟶5,where
(108)A⟶5=e¯2I⟶Nt⊗I⟶Nz−1⊗1hx2D⟶x22:Nx,:+I⟶Nt⊗1hz2D⟶z22:Nz,:⊗I⟶Nx−1−diage¯1UsI⟶Nt⊗I⟶Nz−1⊗1hxD⟶x2:Nx,:−diage¯1WsI⟶Nt⊗1hzD⟶z2:Nz,:⊗I⟶Nx−1+diagϱ3+ϱ4−1htD⟶t1:Nt,:⊗I⟶Nz−1⊗I⟶Nx−1,and
(109a)K⟶5=−ψijk5+ψ~ijk5−σijk7+σ~ijk7+δijk5−diagρ5,(109b)B⟶5=−ζ¯T~xxs+T~zzs.The spectral-collocated field variables at the initial and boundary points are then given by(110a)U~⟶s+1:,:,1=0,T~s+1:,:,1=Ti,Θ~s+1:,:,1=θ0,C~s+1:,:,1=Ci,(110b)U~s+1Nx+1,:,:=UwmU0ϕ,W~s+11,:,:=0,W~s+1Nx+1,:,:=0,U~s+11,:,:=UwmU0ϕ,(110c)U~s+1:,Nz+1,:=0,U~s+1:,1,:=0,W~s+1:,Nz+1,:=0,W~s+1:,1,:=VwmU02ϕ,(110d)I⟶Nt+1⊗I⟶Nz+1⊗D⟶Nx+1,:,:1T~s+1=q′dT0−Ti,I⟶Nt+1⊗I⟶Nz+1⊗D⟶1,:,:1T~s+1=q′dT0−Ti,(110e)T~s+1:,Nz+1,:=−TiT0−Ti,T~s+1:,1,:=T0−2TiT0−Ti,C~s+1:,Nz+1,:=0,(110f)I⟶Nt+1⊗I⟶Nz+1⊗D⟶Nx+1,:,:1C~s+1=0,I⟶Nt+1⊗I⟶Nz+1⊗D⟶1,:,:1C~s+1=0,C~s+1:,1,:=1,(110g)I⟶Nt+1⊗I⟶Nz+1⊗D⟶Nx+1,:,:1Θ~s+1=Exdθs−θr,I⟶Nt+1⊗I⟶Nz+1⊗D⟶1,:,:1Θ~s+1=Exdθs−θr,(110h)Θ~s+1:,Nz+1,:=0,Θ~s+1:,1,:=θ0−θrθs−θr.Then, starting from initial guesses, we obtain approximate solutions forux,z,t, wx,z,t, Tx,z,t, θx,z,t, and Cx,z,t by iteratively solving the matrix equations (73), (81), (91), (99), and (107), in turn, for s=0,1,2,⋯. We only confirm convergence to the expected solutions by measuring the norm residual errors of the difference between successive iterations as defined in (64).


## 3.1. MV-SLQLM Approximation and Construction in (2+1) Dimensions

Suppose we have a general nonlinear PDE of the form(30)∂∂t=F∂2u∂x2,∂2u∂z2,∂u∂x,∂u∂z,u,for x,z,t∈a,b×c,d×0,T, where F is a nonlinear operator.To solve equation (30), first, we require the QLM to linearize any nonlinear equations [67]. Then, we approximate the unknown solution ux,z,t of each linear equation using the triple Lagrange interpolation functions Lix, Ljz, and Lkt at the Chebyshev-Gauss-Lobatto (CGL) collocation points x~,z~,t~ located in the local support domain −1,13. Like in the generalized differential quadrature (GDQ) method [81], the use of Gauss-Lobatto (GL) rule permits grid points to be located on the boundaries, appropriate for the application of boundary conditions. Moreover, grid points with nonuniform intervals provide more accurate results in the calculation of numerical integrals and derivatives.As with other related meshless methods [49, 82], accuracy and efficiency of any spectral method depend on the number of nodes (collocation points) in the local support domain. We thus obtain the spectral domain −1,13 from the original domain x,z,t∈a,b×c,d×0,T via the linear transformation,
(31)x=12b−ax~+12a+b,z=12d−cz~+12c+d,andt=12t~+1,where [55]
(32)x~i=cosπiNxi=0Nx,z~j=cosπjNzj=0Nz,t~k=cosπkNtk=0Nt,is the collocation points described by the Chebyshev polynomial TNv of order N∈ℕ,
(33)TNv=cosNarccosv.Then, we approximate each functionux,z,t by the Chebyshev interpolation uNx,z,t, where
(34)UNx,z,t=∑k=0Nt∑j=0Nz∑i=0NxukjiLktLjzLix,together with the corresponding approximate derivatives at the collocation points,
(35a)∂n∂xnux~i,z~j,t~k=∑r=0Nt∑q=0Nz∑p=0NxukjiLrt~kLqz~jDpin,(35b)∂n∂znux~i,z~j,t~k=∑r=0Nt∑q=0Nz∑p=0NxukjiLrt~kDqjnLpx~i,(35c)∂n∂tnux~i,z~j,t~k=∑r=0Nt∑q=0Nz∑p=0NxukjiDrknLqz~jLpx~i,where Dpin=Linx~i, Dqjn=Ljnz~j, and Drkn=Lknt~k are, respectively, Nx+1×Nx+1, Nz+1×Nz+1, and Nt+1×Nt+1 Chebyshev differentiation matrices for i,j,k=0,1,⋯,N.Further, a Lagrange polynomial of orderN, for example, at the Chebyshev-Gauss-Lobatto points Liv, is defined as [72]
(36)Liv=−1i+11−v2TN′vσ¯iN2v−vi,i=0,1,2,⋯,N,where σ¯i=1 for i=1,2,⋯,N−1, and σ¯0=σ¯N=2. In addition, the first order Chebyshev differentiation matrix at the collocation points is given by [52, 55]
(37)Dij1=σ¯j−1j+iσ¯ivj−vi,i≠j,−vi21−vi2,i=j≠0,N,2N2+16,i,j=0,2N2+16,i,j=N.For simplicity and efficiency, we express equations (34), (35a), (35b), and (35c) in matrix form, using the Kronecker product, as [72]
(38a)ux~i,z~j,t~k=Lt~k⊗Lz~j⊗Lx~iU,(38b)∂n∂xnux~i,z~j,t~k=Lt~k⊗Lz~j⊗D⟶xnU,(38c)∂n∂znux~i,z~j,t~k=Lt~k⊗D⟶zn⊗Lx~iU,(38d)∂n∂tnux~i,z~j,t~k=D⟶tn⊗Lz~j⊗Lx~iU,where U=u000,⋯,u00Nx,u010,⋯,u01Nx,⋯,u0NxNz,u100,⋯,uNtNzNxT,
(39a)Lt~k=L0t~k,L1t~k,L2t~k,⋯,LNtt~k,(39b)Lz~j=L0z~j,L1z~j,L2z~j,⋯,LNzz~j,(39c)Lx~i=L0x~i,L1x~i,L2x~i,⋯,LNxx~i,and
(40)Lt~k⊗Lz~j⊗Lx~i=L0t~kL0z~jLx~i,⋯,L0t~kLNzz~jLx~i,L1t~kL0z~jLx~i,⋯,L1t~kLNzz~jLx~i,⋯,LNtt~kLNzz~jLx~i.Third, we formulate the spectral collocated coupled system, defined for each equation, as(41a)U~s+1=Us,Ws,Uxs,Uxxs,Uzs,Uzzs,Ws+1,(41b)W~s+1=Ws,Us,Wxs,Wxxs,Wzs,Wzzs,Us+1,(41c)T~s+1=Ts,Us,Ws,Txs,Txxs,Tzs,Tzzs,Θs+1,(41d)Θ~s+1=Θs,Θxs,Θxxs,Θzs,Θzzs,Ts+1,(41e)C~s+1=Cs,Us,Ws,Cxs,Cxxs,Czs,Czzs,Ts+1.We apply the LLM algorithm to decouple the spectral collocated system (41a), (41b), (41c), (41d), and (41e) and solve each equation as follows:
(a)
Solve forU~s+1 in equation (41a), assuming Ws and Ws+1 are known from the previous iteration(b)
WithU~s+1 now known, solve for W~s+1 in equation (41b)(c)
Solve forT~s+1 in (41c), using U~s+1 and W~s+1, assuming Θs+1 is known from the previous iteration(d)
WithT~s+1 now known, proceed to solve for Θ~s+1 in (41d)(e)
Finally, usingU~s+1, W~s+1, and T~s+1, solve for C~s+1 in (41e)Beginning with the initial guesses, steps (a) to (e) are repeated, in turn, fors=0,1,2,⋯ to derive approximate solutions to the given model.


## 3.2. Error, Stability, and Convergence Analysis

### 3.2.1. Error Bound

Following the approach in Bhrawy [83] and Samuel and Motsa [84], in this subsection, we present an analytical expression for the error norm and an upper bound on the expected error of the approximation, equation (34), to the smooth function ux,z,t∈Ω≡a,b×c,d×0,T. Let
(42)∏=SpanLixLjzLkt,∀i=0,1,2,⋯,Nx;j=0,1,2,⋯,Nz;k=0,1,2,⋯,Nt.Given a smooth functionux,z,t∈Ω≡a,b×c,d×0,T. Let UNx,z,t, defined in (34), be the best approximation of ux,z,t so that ∀VN∈∏,
(43)ux,z,t−UNx,z,t∞≤ux,z,t−VNx,z,t∞.Then, there exists arbitrary pointsγ1,γ1′∈a,b, γ2,γ2′∈c,d, and γ3,γ3′∈0,T, so that
(44)ux,z,t−VNx,z,t=∂Nx+1uγ1,z,t∂xNx+1Nx+1!∏i=0Nxx−xi+∂Nz+1ux,γ2,t∂zNz+1Nz+1!∏j=0Nzz−zj+∂Nt+1ux,z,γ3∂tNt+1Nt+1!∏k=0Ntt−tk−∂Nx+Nz+Nt+3uγ1,γ2,γ3∂xNx+1∂zNz+1∂tNt+1Nt+1!Nz+1!Nx+1!×∏i=0Nxx−xi×∏j=0Nzz−zj×∏k=0Ntt−tk,where VNx,z,t is the multivariate chebyshev interpolating polynomial of ux,z,t at xi, zj, and tk.Taking the infinite norm of equation (44) yields
(45)ux,z,t−VNx,z,t∞≤maxx,z,t∈Ω∂Nx+1uγ1,z,t∂xNx+1∏i=0Nxx−xiNx+1!+maxx,z,t∈Ω∂Nz+1ux,γ2,t∂zNz+1∏j=0Nxz−zjNz+1!+maxx,z,t∈Ω∂Nt+1ux,z,γ3∂tNt+1∏k=0Ntt−tkNt+1!−maxx,z,t∈Ω∂Nx+Nz+Nt+3uγ1′,γ2′,γ3′∂xNx+1∂zNz+1∂tNt+1×∏i=0Nxx−xi×∏j=0Nxz−zj×∏k=0Ntt−tkNx+1!Nz+1!Nt+1!.Sinceux,z,t is a smooth function on Ω; then, its Nx+1, Nz+1, Nt+1, and Nx+Nz+Nt+3-th partial derivatives are both continuous and upper-bounded, so that
(46)maxx,z,t∈Ω∂Nx+1uγ1,z,t∂xNx+1≤K1,maxx,z,t∈Ω∂Nz+1ux,γ2,t∂zNz+1≤K2,maxx,z,t∈Ω∂Nt+1ux,z,γ3∂tNt+1≤K3,maxx,z,t∈Ω∂Nx+Nz+Nt+3uγ1′,γ2′,γ3′∂xNx+1∂zNz+1∂tNt+1≤K4,where K1,K2,K3,andK4 are arbitrary constants in ℝ+.Substituting (46) into (45) yields the norm error bound estimate
(47)ux,z,t−VNx,z,t∞≤K1∏i=0Nxx−xi∞Nx+1!+K2∏j=0Nzz−zj∞Nz+1!+K3∏k=0Ntt−tk∞Nt+1!−K4×∏i=0Nxx−xi∞×∏j=0Nzz−zj∞×∏k=0Ntt−tk∞Nt+1!Nz+1!Nx+1!,where VN∈∏Nx,Nz,Nt.However, (47) can be further simplified. From the linear relationships defined in (31), we have
(48)x−xi=12b−ax~−x~i,z−zj=12d−cz~−z~j,andt−tk=12t~−t~k,so that from (47),
(49)∏i=0Nxx−xi∞=∏i=0Nxb−a2x~−x~i∞≤b−a2Nxmaxx~,xi∈−1,1∏i=0Nxx~−x~i≤b−aNx+12Nx+1maxx~,xi∈−1,1∏i=0Nxx~−x~i.Using equation (36),
(50)Lix~=−1i+11−x~2TNxx~′σ¯iNx2x~−x~i,i=0,1,2,⋯,Nx,so that the Nx+1-th degree Lagrange polynomial is
(51)LNx+1x~=−1Nx+21−x~2TNxx~′σ¯iNx2x~−x~i.Taking the modulus of equation (51) and rearranging the terms yield
(52)x~−x~i=1−x~2TNxx~′Nx2LNx+1x~.Since the Chebyshev interpolating polynomial, equation (33) is the solution to the Chebyshev differential equation [85]
(53)1−x2d2ydx2−xdydx+n2y=0,for nonnegative integers n,n∈ℤ0 and x<1, it follows that
(54)1−x2TNxx′=−xTNxx−n2.Using equation (54) and applying the triangular inequality, equation (52) can then be expressed as
(55)x~−x~i≤−x~TNxx~+−n2Nx2LNx+1x~.Employing the inequality relation stated in [84] and noting that TNxx~≤1,∀xi∈−1,1, then
(56)1−x~2TNx′x~=−x~TNxx~+−n2≤2Nx,so that (55) becomes
(57)x~−x~i≤2Nx.Therefore, from (57), maxx~,xi∈−1,1∏i=0Nxx~−x~i≤2Nx+2/Nx. It then also follows from (49) that
(58)∏i=0Nxx−xi∞≤b−aNx+12Nx+1maxx~,xi∈−1,1∏i=0Nxx~−x~i≤2Nx+2b−aNx+1Nx2Nx+1.In a similar manner,(59)∏j=0Nzz−zj∞≤d−cNz+12Nz+1maxz~,zj∈−1,1∏j=0Nzz~−z~i≤2Nz+2d−cNz+1Nz2Nz+1,(60)∏k=0Ntt−tk∞≤12Nt+1maxt~,tk∈−1,1∏k=0Ntt~−t~k.Substituting (58), (59), and (60) into (47) yields
(61)ux,z,t−VNx,z,t∞≤K12Nx+2b−aNx+1Nx2Nx+1Nx+1!+K22Nz+2d−cNz+1Nz2Nz+1Nz+1!+K32Nt+2Nt2Nt+1Nt+1!−K4×2Nx+2b−aNx+1×2Nz+2d−cNz+1×2Nt+2NxNzNt2Nx+Nz+Nt+3Nt+1!Nz+1!Nx+1!.Hence, by combining (61) and (43), an upper bound of the error for the approximate solutions is obtained. Since the model, equations (10)–(13) and (19), does not have the exact solution, it is inferred that the convergence of the proposed method depends basically on the above error bound.


### 3.2.2. Stability

According to Noborio et al. [30], the display of dispersions or oscillatory errors in the numerical solutions is completely unavoidable regardless of whether one uses FD or FE methods to simulate transport due to the characteristics of the solute transport equation. Hence, to ensure a generally stable MV-SLQLM, two complementary measures are employed on the basis of solute transport equation, equation (13). First, in reference [86], we investigated the existence and uniqueness of the classical solution to equation (13) using the proposed “VF-EM” technique. Defining an error e~≔ct−c~t and making inference of Ilati and Dehghan [82], the result (see equation 3.10 in [86]) means e~ is strongly stable. Moreover, the well-posed analysis of equation (13) in time integral also demonstrates that the numerical technique is stable in L2-norm. Second, we examine numerically the stability of a discretized scheme for convection-diffusion terms in equation (13). Among the existing numerical methods for analyzing scheme’s stability on the basis of convection-diffusion terms, Yu et al. [87] observed that majority are usually based on five assumptions, i.e., one-dimensional, linear, first kind boundary condition, free source term, and uniform grid system.However, from [87], any deviation from the aforementioned five assumptions may attribute to enhanced convective stability. In this study, except for the linear nature, equation (13) is of (2+1) dimensions, together with a nonfree source term, and boundary conditions of first and second types. It is worth also mentioning that the convective terms of equations (12) and (13) use flow velocity functions derived from the pair of highly nonlinear momentum equations, equations (10) and (11). Moreover, the MV-SLQLM is built on Lagrange interpolation bases that do not depend on the order of nodes (nonuniform) and are optimally stable [88, 89]. In spite of the fact that these properties naturally contribute to a more stable numerical scheme, we further check the stability of the convective-diffusion solution profiles following the study of Abdullah [90]. Using a one-dimension steady convective-diffusion equation, Abdullah [90] described a linear relationship between grid size and local Pe number;
(62)Nmin=mPe+b,3.125<Pe≤100,3,0≤Pe≤3.125,where Nmin is the minimum grid size, sufficient for the flow profile to behave physically correctly; Pe is the grid or local Peclet number; and m and b are the slope and intercept constants, respectively. Since we have a low Pe, 0≤Pe≤100, advection-diffusion-reaction equation (13), we investigate the effect of varying Pe on the solute concentration profiles at the minimum grid number, Nmin.


### 3.2.3. Convergence

Making inference of (61) and (43), convergence to the expected solution is affirmed by considering the L∞-norm of the difference between successive iterations, using the residual error analysis
(63)EU~=max0≤i≤NtU~s+1,i−U~s,i∞,EW~=max0≤i≤NtW~s+1,i−W~s,i∞,ET~=max0≤i≤NtT~s+1,i−T~s,i∞,EΘ~=max0≤i≤NtΘ~s+1,i−Θ~s,i∞,EC~=max0≤i≤NtC~s+1,i−C~s,i∞.,defined for the difference between values of successive iterations.


## 3.2.1. Error Bound

Following the approach in Bhrawy [83] and Samuel and Motsa [84], in this subsection, we present an analytical expression for the error norm and an upper bound on the expected error of the approximation, equation (34), to the smooth function ux,z,t∈Ω≡a,b×c,d×0,T. Let
(42)∏=SpanLixLjzLkt,∀i=0,1,2,⋯,Nx;j=0,1,2,⋯,Nz;k=0,1,2,⋯,Nt.Given a smooth functionux,z,t∈Ω≡a,b×c,d×0,T. Let UNx,z,t, defined in (34), be the best approximation of ux,z,t so that ∀VN∈∏,
(43)ux,z,t−UNx,z,t∞≤ux,z,t−VNx,z,t∞.Then, there exists arbitrary pointsγ1,γ1′∈a,b, γ2,γ2′∈c,d, and γ3,γ3′∈0,T, so that
(44)ux,z,t−VNx,z,t=∂Nx+1uγ1,z,t∂xNx+1Nx+1!∏i=0Nxx−xi+∂Nz+1ux,γ2,t∂zNz+1Nz+1!∏j=0Nzz−zj+∂Nt+1ux,z,γ3∂tNt+1Nt+1!∏k=0Ntt−tk−∂Nx+Nz+Nt+3uγ1,γ2,γ3∂xNx+1∂zNz+1∂tNt+1Nt+1!Nz+1!Nx+1!×∏i=0Nxx−xi×∏j=0Nzz−zj×∏k=0Ntt−tk,where VNx,z,t is the multivariate chebyshev interpolating polynomial of ux,z,t at xi, zj, and tk.Taking the infinite norm of equation (44) yields
(45)ux,z,t−VNx,z,t∞≤maxx,z,t∈Ω∂Nx+1uγ1,z,t∂xNx+1∏i=0Nxx−xiNx+1!+maxx,z,t∈Ω∂Nz+1ux,γ2,t∂zNz+1∏j=0Nxz−zjNz+1!+maxx,z,t∈Ω∂Nt+1ux,z,γ3∂tNt+1∏k=0Ntt−tkNt+1!−maxx,z,t∈Ω∂Nx+Nz+Nt+3uγ1′,γ2′,γ3′∂xNx+1∂zNz+1∂tNt+1×∏i=0Nxx−xi×∏j=0Nxz−zj×∏k=0Ntt−tkNx+1!Nz+1!Nt+1!.Sinceux,z,t is a smooth function on Ω; then, its Nx+1, Nz+1, Nt+1, and Nx+Nz+Nt+3-th partial derivatives are both continuous and upper-bounded, so that
(46)maxx,z,t∈Ω∂Nx+1uγ1,z,t∂xNx+1≤K1,maxx,z,t∈Ω∂Nz+1ux,γ2,t∂zNz+1≤K2,maxx,z,t∈Ω∂Nt+1ux,z,γ3∂tNt+1≤K3,maxx,z,t∈Ω∂Nx+Nz+Nt+3uγ1′,γ2′,γ3′∂xNx+1∂zNz+1∂tNt+1≤K4,where K1,K2,K3,andK4 are arbitrary constants in ℝ+.Substituting (46) into (45) yields the norm error bound estimate
(47)ux,z,t−VNx,z,t∞≤K1∏i=0Nxx−xi∞Nx+1!+K2∏j=0Nzz−zj∞Nz+1!+K3∏k=0Ntt−tk∞Nt+1!−K4×∏i=0Nxx−xi∞×∏j=0Nzz−zj∞×∏k=0Ntt−tk∞Nt+1!Nz+1!Nx+1!,where VN∈∏Nx,Nz,Nt.However, (47) can be further simplified. From the linear relationships defined in (31), we have
(48)x−xi=12b−ax~−x~i,z−zj=12d−cz~−z~j,andt−tk=12t~−t~k,so that from (47),
(49)∏i=0Nxx−xi∞=∏i=0Nxb−a2x~−x~i∞≤b−a2Nxmaxx~,xi∈−1,1∏i=0Nxx~−x~i≤b−aNx+12Nx+1maxx~,xi∈−1,1∏i=0Nxx~−x~i.Using equation (36),
(50)Lix~=−1i+11−x~2TNxx~′σ¯iNx2x~−x~i,i=0,1,2,⋯,Nx,so that the Nx+1-th degree Lagrange polynomial is
(51)LNx+1x~=−1Nx+21−x~2TNxx~′σ¯iNx2x~−x~i.Taking the modulus of equation (51) and rearranging the terms yield
(52)x~−x~i=1−x~2TNxx~′Nx2LNx+1x~.Since the Chebyshev interpolating polynomial, equation (33) is the solution to the Chebyshev differential equation [85]
(53)1−x2d2ydx2−xdydx+n2y=0,for nonnegative integers n,n∈ℤ0 and x<1, it follows that
(54)1−x2TNxx′=−xTNxx−n2.Using equation (54) and applying the triangular inequality, equation (52) can then be expressed as
(55)x~−x~i≤−x~TNxx~+−n2Nx2LNx+1x~.Employing the inequality relation stated in [84] and noting that TNxx~≤1,∀xi∈−1,1, then
(56)1−x~2TNx′x~=−x~TNxx~+−n2≤2Nx,so that (55) becomes
(57)x~−x~i≤2Nx.Therefore, from (57), maxx~,xi∈−1,1∏i=0Nxx~−x~i≤2Nx+2/Nx. It then also follows from (49) that
(58)∏i=0Nxx−xi∞≤b−aNx+12Nx+1maxx~,xi∈−1,1∏i=0Nxx~−x~i≤2Nx+2b−aNx+1Nx2Nx+1.In a similar manner,(59)∏j=0Nzz−zj∞≤d−cNz+12Nz+1maxz~,zj∈−1,1∏j=0Nzz~−z~i≤2Nz+2d−cNz+1Nz2Nz+1,(60)∏k=0Ntt−tk∞≤12Nt+1maxt~,tk∈−1,1∏k=0Ntt~−t~k.Substituting (58), (59), and (60) into (47) yields
(61)ux,z,t−VNx,z,t∞≤K12Nx+2b−aNx+1Nx2Nx+1Nx+1!+K22Nz+2d−cNz+1Nz2Nz+1Nz+1!+K32Nt+2Nt2Nt+1Nt+1!−K4×2Nx+2b−aNx+1×2Nz+2d−cNz+1×2Nt+2NxNzNt2Nx+Nz+Nt+3Nt+1!Nz+1!Nx+1!.Hence, by combining (61) and (43), an upper bound of the error for the approximate solutions is obtained. Since the model, equations (10)–(13) and (19), does not have the exact solution, it is inferred that the convergence of the proposed method depends basically on the above error bound.


## 3.2.2. Stability

According to Noborio et al. [30], the display of dispersions or oscillatory errors in the numerical solutions is completely unavoidable regardless of whether one uses FD or FE methods to simulate transport due to the characteristics of the solute transport equation. Hence, to ensure a generally stable MV-SLQLM, two complementary measures are employed on the basis of solute transport equation, equation (13). First, in reference [86], we investigated the existence and uniqueness of the classical solution to equation (13) using the proposed “VF-EM” technique. Defining an error e~≔ct−c~t and making inference of Ilati and Dehghan [82], the result (see equation 3.10 in [86]) means e~ is strongly stable. Moreover, the well-posed analysis of equation (13) in time integral also demonstrates that the numerical technique is stable in L2-norm. Second, we examine numerically the stability of a discretized scheme for convection-diffusion terms in equation (13). Among the existing numerical methods for analyzing scheme’s stability on the basis of convection-diffusion terms, Yu et al. [87] observed that majority are usually based on five assumptions, i.e., one-dimensional, linear, first kind boundary condition, free source term, and uniform grid system.However, from [87], any deviation from the aforementioned five assumptions may attribute to enhanced convective stability. In this study, except for the linear nature, equation (13) is of (2+1) dimensions, together with a nonfree source term, and boundary conditions of first and second types. It is worth also mentioning that the convective terms of equations (12) and (13) use flow velocity functions derived from the pair of highly nonlinear momentum equations, equations (10) and (11). Moreover, the MV-SLQLM is built on Lagrange interpolation bases that do not depend on the order of nodes (nonuniform) and are optimally stable [88, 89]. In spite of the fact that these properties naturally contribute to a more stable numerical scheme, we further check the stability of the convective-diffusion solution profiles following the study of Abdullah [90]. Using a one-dimension steady convective-diffusion equation, Abdullah [90] described a linear relationship between grid size and local Pe number;
(62)Nmin=mPe+b,3.125<Pe≤100,3,0≤Pe≤3.125,where Nmin is the minimum grid size, sufficient for the flow profile to behave physically correctly; Pe is the grid or local Peclet number; and m and b are the slope and intercept constants, respectively. Since we have a low Pe, 0≤Pe≤100, advection-diffusion-reaction equation (13), we investigate the effect of varying Pe on the solute concentration profiles at the minimum grid number, Nmin.


## 3.2.3. Convergence

Making inference of (61) and (43), convergence to the expected solution is affirmed by considering the L∞-norm of the difference between successive iterations, using the residual error analysis
(63)EU~=max0≤i≤NtU~s+1,i−U~s,i∞,EW~=max0≤i≤NtW~s+1,i−W~s,i∞,ET~=max0≤i≤NtT~s+1,i−T~s,i∞,EΘ~=max0≤i≤NtΘ~s+1,i−Θ~s,i∞,EC~=max0≤i≤NtC~s+1,i−C~s,i∞.,defined for the difference between values of successive iterations.


## 3.3. Numerical Application

Now, we apply the above described procedure, the MV-SLQLM, to the model under consideration, equations (10)–(13) and (19), together with (20a), (20b), (20c), (20d), (20e), and (20f). Beginning with the nonlinear Darcy-Brinkman-Forchheimer equation (equation (10)) in the x−direction
(64)∂u∂t+ϕuux+ϕwuz=a1uxx+uzz−a2u−a3u2,where a1=1/Re, a2=1/Re·Λ, and a3=F∗/Re, we employ the QLM to look for the linear expression
(65)α2,sx,z,tus"+1+α1,sx,z,tus′+1+α0,sx,z,tus+1−uts+1=Rsx,z,t,where
(66a)Rsx,z,t=∑m=02αm,sx,z,tusm−Fus,us′,us",(66b)Fus,us′,us"=a1uxxs+uzzs−ϕuux−ϕwuz−a2us−a3us2,(66c)αm,sx,z,t=∂Fus,us′,us"∂usm,m=0,1,2.Using equations (66b) and (66c), we have
(67)α0,sx,z,t=−ϕuxs−a2−2a3us;α1,sx,z,t=−ϕus−ϕws;α2,sx,z,t=a1,so that(68)Rsx,z,t=−ϕuxsus+ϕwsuzs+ϕws+1uzs−a3us2.Hence, equation (65) becomes
(69)a1uxxs+1+uzzs+1−ϕusuxs+1−ϕwsuzs+1−ϕuxs+a2+2a3usus+1−uts+1−ϕuzsws+1=−ϕuxsus+wsuzs−a3us2.Applying the Chebyshev spectral collocation method to equation (69) yields
(70)a1I⟶Nt+1⊗I⟶Nz+1⊗1hx2D⟶x2+I⟶Nt+1⊗1hz2D⟶z2⊗I⟶Nx+1−ϕdiagUsI⟶Nt+1⊗I⟶Nz+1⊗1hxD⟶x−ϕdiagWsI⟶Nt+1⊗1hzD⟶z⊗I⟶Nx+1−diagϕUxs+a2+2a3Us−1htD⟶t⊗I⟶Nz+1⊗I⟶Nx+1Us+1−ϕdiagUzsWs+1=−ϕUxsUs+UzsWs−a3UsUs,where hx=xl−x0/2, hz=zl−z0/2, and ht=T/2; and I⟶Nt+1, I⟶Nx+1, and I⟶Nz+1 are, respectively, identity matrices of sizes Nt+1×Nt+1, Nx+1×Nx+1, and Nz+1×Nz+1.SinceU in equation (70) is known at the boundaries of the domain, rewriting equation (70) for the interior points of the domain, yields the new expression
(71)a1I⟶Nt⊗I⟶Nz−1⊗1hx2D⟶x22:Nx,:+I⟶Nt⊗1hz2D⟶z22:Nz,:⊗I⟶Nx−1−diagϕUxs+a2+2a3Us−ϕdiagUsI⟶Nt⊗I⟶Nz−1⊗1hxD⟶x2:Nx,:−ϕdiagWsI⟶Nt⊗1hzD⟶z2:Nz,:⊗I⟶Nx−1−1htD⟶t1:Nt,:⊗I⟶Nz−1⊗I⟶Nx−1U~s+1−ϕdiagUzsWs+1=−ϕUxsUs+UzsWs−a3U~sU~s−ψijk1+ψ~ijk1−σijk1+σ~ijk1+δijk1,where U~=U2:Nx,2:Nz,1:Nt, and
(72a)ψ1=a1D⟶x22:Nx,1⊗U1,2:Nz,1:Nt−UsD⟶x2:Nx,1⊗U1,2:Nz,1:Nt,(72b)ψ~ijk1=a1D⟶x22:Nx,Nx+1⊗UNx+1,2:Nz,1:Nt−UsD⟶x2:Nx,Nx+1⊗UNx+1,2:Nz,1:Nt,(72c)σijk1=a1D⟶z22:Nz,1⊗U2:Nx,1,1:Nt−WsD⟶z2:Nz,1⊗U2:Nx,1,1:Nt,(72d)σ~ijk1=a1D⟶z22:Nz,Nz+1⊗U2:Nx,Nz+1,1:Nt−WsD⟶z2:Nz,1⊗U2:Nx,1,1:Nt,(72e)δijk1=−D⟶t1:Nt,Nt+1⊗U2:Nx,2:Nz,Nt+1.Thus, at step (a) of the LLM, we need to solve forU~s+1 in the matrix equation (equation (71)) defined as
(73)A⟶1U~s+1=B⟶1+K⟶1,where
(74)A⟶1=a1I⟶Nt⊗I⟶Nz−1⊗1hx2D⟶x22:Nx,:+I⟶Nt⊗1hz2D⟶z22:Nz,:⊗I⟶Nx−1−diagϕUxs+a2+2a3Us−ϕdiagUsI⟶Nt⊗I⟶Nz−1⊗1hxD⟶x2:Nx,:−ϕdiagWsI⟶Nt⊗1hzD⟶z2:Nz,:⊗I⟶Nx−1−1htD⟶t1:Nt,:⊗I⟶Nz−1⊗I⟶Nx−1,and
(75a)B⟶1=−ϕUxsUs+UzsWs−a3U~sU~s+ϕdiagUzsWs+1,(75b)K⟶1=−ψijk1+ψ~ijk1−σijk1+σ~ijk1+δijk1.In a similar manner, the Darcy-Brinkman-Forchheimer equation (equation (11)) in z−direction is expressed as
(76)∂w∂t=−b0Pz−ϕuwx−ϕwwz+b1wxx+wzz−b2w−b3w2+b4,where b0=1/ϕ, b1=1/Re, b2=1/Re·Λ, b3=F∗/Re, and b4=1/ϕGrT+GrC/Re2.Applying the QLM, we obtain(77)b1wxxs+1+wzzs+1−ϕuswxs+1−ϕwswzs+1−ϕwzs+b2+2b3wsws+1−b0Pzs+1−wts+1−ϕwxsus+1=−ϕwxus+wzsws−b3wsws−b4.Then, by applying the Chebyshev spectral collocation method to equation (77), we get
(78)b1I⟶Nt+1⊗I⟶Nz+1⊗1hx2D⟶x2+I⟶Nt+1⊗1hz2D⟶z2⊗I⟶Nx+1−ϕdiagUsI⟶Nt+1⊗I⟶Nz+1⊗1hxD⟶x−ϕdiagWsI⟶Nt+1⊗1hzD⟶z⊗I⟶Nx+1−diagϕWzs+b2+2b3Ws−1htD⟶t⊗I⟶Nz+1⊗I⟶Nx+1Ws+1−b0Pzs+1−ϕdiagWxsUs+1=−ϕWxsUs+WzsWs−b3WsWs−b4,where hx, hz, ht, I⟶Nt+1, I⟶Nx+1, and I⟶Nz+1 are all defined as before.Considering the known values ofW at the initial and boundary conditions, equation (78) can be rearranged as
(79)b1I⟶Nt⊗I⟶Nz−1⊗1hx2D⟶x22:Nx,:+I⟶Nt⊗1hz2D⟶z22:Nz,:⊗I⟶Nx−1−ϕdiagUsI⟶Nt⊗I⟶Nz−1⊗1hxD⟶x2:Nx,:−ϕdiagWsI⟶Nt⊗1hzD⟶z2:Nz,:⊗I⟶Nx−1−diagϕWzs+b2+2b3Ws−1htD⟶t1:Nt,:⊗I⟶Nz−1⊗I⟶Nx−1W~s+1=b0P~zs+1+ϕdiagWxsU~s+1−ϕW~xsU~s+W~zsW~s−b3W~sW~s−ψijk2+ψ~ijk2−σijk2+σ~ijk2+δijk2−b4,where W~=W2:Nx,2:Nz,1:Nt, and
(80a)ψ2=b1D⟶x22:Nx,1⊗W1,2:Nz,1:Nt−UsD⟶x2:Nx,1⊗W1,2:Nz,1:Nt,(80b)ψ~ijk2=b1D⟶x22:Nx,Nx+1⊗WNx+1,2:Nz,1:Nt−UsD⟶x2:Nx,Nx+1⊗WNx+1,2:Nz,1:Nt,(80c)σijk2=b1D⟶z22:Nz,Nz+1⊗U2:Nx,Nz+1,1:Nt−WsD⟶z2:Nz,1⊗W2:Nx,1,1:Nt,(80d)σ~ijk2=b1D⟶z22:Nz,Nz+1⊗U2:Nx,Nz+1,1:Nt−WsD⟶z2:Nz,Nz+1⊗W2:Nx,Nz+1,1:Nt,(80e)δijk2=−D⟶t1:Nt,Nt+1⊗W2:Nx,2:Nz,Nt+1.Thus, in step (b) of the LLM, we solve forW~s+1 in the matrix equation
(81)A⟶2W~s+1=B⟶2+K⟶2,where
(82)A⟶2=b1I⟶Nt⊗I⟶Nz−1⊗1hx2D⟶x22:Nx,:+I⟶Nt⊗1hz2D⟶z22:Nz,:⊗I⟶Nx−1−ϕdiagUsI⟶Nt⊗I⟶Nz−1⊗1hxD⟶x2:Nx,:−ϕdiagWsI⟶Nt⊗1hzD⟶z2:Nz,:⊗I⟶Nx−1−diagϕWzs+b2+2b3Ws−1htD⟶t1:Nt,:⊗I⟶Nz−1⊗I⟶Nx−1,and
(83a)B⟶2=b0P~zs+1+ϕdiagWxsU~s+1−ϕW~xsU~s+W~zsW~s−b3W~sW~s,(83b)K⟶2=−ψijk2+ψ~ijk2−σijk2+σ~ijk2+δijk2−b4.The pressure termP~zs+1 is iteratively estimated from an initial guess
(84)Px,z=x4−4π2sinh2πsinhπxcosπz,deduced from the analytical solution [91]
(85)Px,z=x4−4∑n=1,odd∞1nπ2sinh2nπsinhnπxcosnπz,of the pressure equation ∇2P=0, assuming an initial state of P=0 everywhere and the zero-flux top and bottom boundary conditions.Next, we express the heat equation as(86)∂T∂t=c¯1∂2T∂x2+∂2T∂z2−c¯2u∂T∂x+w∂T∂z+c¯3∂2θ∂x2+∂2θ∂z2+c¯4T+ϕ2EcReϕ∂u∂x2,where c¯1=1/PrRe, c¯2=ϱ1ϕ2, c¯3=ϕEcϱ2L0Dθ/Reϕe, c¯4=δ1/PrRe, and c¯5=ϕ2Ec/Re.In terms of the QLM, (86) can be expressed as
(87)c¯1Txxs+1+Tzzs+1−c¯2usTxs+1−c¯2wsTzs+1+c¯3θxxs+1+θzzs+1+c¯4Ts+1−Tts+1=−ϕ2c¯5uxs2.Applying the Chebyshev spectral collocation method to equation (87), we get
(88)c¯1I⟶Nt+1⊗I⟶Nz+1⊗1hx2D⟶x2+I⟶Nt+1⊗1hz2D⟶z2⊗I⟶Nx+1−diagc¯2UsI⟶Nt+1⊗I⟶Nz+1⊗1hxD⟶x−diagc¯2WsI⟶Nt+1⊗1hzD⟶z⊗I⟶Nx+1+diagc4−1htD⟶t⊗I⟶Nz+1⊗I⟶Nx+1Ts+1+c¯3I⟶Nt+1⊗I⟶Nz+1⊗1hx2D⟶x2+I⟶Nz+1⊗1hz2D⟶z2⊗I⟶Nx+1Θs+1=−ϕ2c¯5Uxs2,where hx, hz, ht, I⟶Nt+1, I⟶Nx+1, and I⟶Nz+1 are all as previously described.In a similar manner, with the known values at the boundaries and initial points, we solve(89)c¯1I⟶Nt⊗I⟶Nz−1⊗1hx2D⟶x22:Nx,:+I⟶Nt⊗1hz2D⟶z22:Nz,:⊗I⟶Nx−1−diagc¯2UsI⟶Nt⊗I⟶Nz−1⊗1hxD⟶x2:Nx,:−diagc¯2WsI⟶Nt⊗1hzD⟶z2:Nz,:⊗I⟶Nx−1+diagc4−1htD⟶t1:Nt,:⊗I⟶Nz−1⊗I⟶Nx−1T~s+1+c¯3Θ~xxs+1+Θ~zzs+1=−ϕ2c¯5Uxs2−ψijk3+ψ~ijk3−σijk3+σ~ijk3+δijk3,where T~=T2:Nx,2:Nz,1:Nt, U~s≈U~s+1, W~s≈W~s+1, Θ~=Θ2:Nx,2:Nz,1:Nt, and
(90a)ψ3=c¯1D⟶x22:Nx,1⊗T1,2:Nz,1:Nt−c¯2U~sD⟶x2:Nx,1⊗T1,2:Nz,1:Nt,(90b)ψ~ijk3=c¯1D⟶x22:Nx,Nx+1⊗TNx+1,2:Nz,1:Nt−c¯2U~sD⟶x2:Nx,Nx+1⊗TNx+1,2:Nz,1:Nt,(90c)σijk3=c¯1D⟶z22:Nz,1⊗T2:Nx,1,1:Nt−c¯2W~sD⟶z2:Nz,1⊗T2:Nx,1,1:Nt,(90d)σ~ijk3=c¯1D⟶z22:Nz,Nz+1⊗T2:Nz,Nz+1,1:Nt−c¯2W~sD⟶z2:Nz,Nz+1⊗T2:Nx,Nz+1,1:Nt,(90e)δijk3=−D⟶t1:Nt,Nt+1⊗T2:Nx,2:Nz,Nt+1.Henceforth, at step (c) of the LLM, we solve the matrix equation,(91)A⟶3T~s+1=B⟶3+K⟶3,where
(92)A⟶3=c¯1I⟶Nt⊗I⟶Nz−1⊗1hx2D⟶x22:Nx,:+I⟶Nt⊗1hz2D⟶z22:Nz,:⊗I⟶Nx−1−diagc¯2UsI⟶Nt⊗I⟶Nz−1⊗1hxD⟶x2:Nx,:−diagc¯2WsI⟶Nt⊗1hzD⟶z2:Nz,:⊗I⟶Nx−1+diagc41htD⟶t1:Nt,:⊗I⟶Nz−1⊗I⟶Nx−1,and,
(93a)K⟶3=−ψijk3+ψ~ijk3−σijk3+σ~ijk3+δijk3,(93b)B⟶3=−c¯3Θ~xxs+Θ~zzs−ϕ2c¯5U~xs2.Since solving the system (91) requires knowledge of θ~s+1, we then consider the moisture equation, which is also coupled with temperature T,
(94)Rn∂θ∂t=θxx+θzz−δ2θz−b¯Txx+Tzz,where b¯=LeT1−T0. Applying the QLM, we need to solve
(95)θxxs+1+θzzs+1−δ2θzs+1−b¯Txxs+Tzzs=Rnθts+1.Then, applying the Chebyshev spectral collocation method to equation (95), we obtain
(96)I⟶Nt+1⊗I⟶Nz+1⊗1hx2D⟶x2+I⟶Nt+1⊗1hz2D⟶z2⊗I⟶Nx+1−diagδ2I⟶Nt+1⊗1hzD⟶z⊗I⟶Nx+1−Rn1htD⟶t⊗I⟶Nz+1⊗I⟶Nx+1Θs+1−b¯I⟶Nt+1⊗I⟶Nz+1⊗1hx2D⟶x2+I⟶Nt+1⊗1hz2D⟶z2⊗I⟶Nx+1Ts+1=0,where hx, hz, ht, I⟶Nt+1, I⟶Nx+1, and I⟶Nz+1 are all as previously described.Taking all the known initial and boundary values ofθs+1 and Ts+1 to the right-hand side of equation (96), we now solve the equation
(97)I⟶Nt⊗I⟶Nz−1⊗1hx2D⟶x22:Nx,:+I⟶Nt⊗1hz2D⟶z22:Nz,:⊗I⟶Nx+1−diagδ2I⟶Nt⊗1hzD⟶z2:Nz,:⊗I⟶Nx−1−Rn1htD⟶t1:Nt,:⊗I⟶Nz−1⊗I⟶Nx−1Θ~s+1−b¯Txxs+1+Tzzs+1=−ψijk4+ψ~ijk4−σijk4+σ~ijk4+δijk4,where
(98a)ψ4=D⟶x22:Nx,1⊗Θ1,2:Nz,1:Nt,(98b)ψ~ijk4=D⟶x22:Nx,Nx+1⊗ΘNx+1,2:Nz,1:Nt,(98c)σijk4=D⟶z22:Nz,1⊗Θ2:Nx,1,1:Nt−δ2D⟶z2:Nz,1⊗θ2:Nx,1,1:Nt,(98d)σ~ijk4=D⟶z22:Nz,Nz+1⊗Θ2:Nz,Nz+1,1:Nt−δ2D⟶z2:Nz,Nz+1⊗Θ2:Nx,Nz+1,1:Nt,(98e)δijk4=−D⟶t1:Nt,Nt+1⊗θ2:Nx,2:Nz,Nt+1,so that, at step (d) of the LLM, we solve the matrix equation
(99)A⟶4Θ~s+1=B⟶4+K⟶4,where
(100)A⟶4=I⟶Nt⊗I⟶Nz−1⊗1hx2D⟶x22:Nx,:+I⟶Nt⊗1hz2D⟶z22:Nz,:⊗I⟶Nx−1−diagδ2I⟶Nt⊗1hzD⟶z2:Nz,:⊗I⟶Nx−1−Rn1htD⟶t1:Nt,:⊗I⟶Nz−1⊗I⟶Nx−1,and
(101a)B⟶4=b¯Txxs+Tzzs,(101b)K⟶4=−ψijk4+ψ~ijk4−σijk4+σ~ijk4+δijk4.Finally, applying the QLM, the solute equation(102)Rm∂C∂t+e¯1u∂C∂x+w∂C∂z=e¯1∂2C∂x2+∂2C∂z2+ζ¯∂2T∂x2+∂2T∂z2+ϱ3C+ϱ4Ca+ϱ5,where e¯1=ϕ2 and e¯2=θw/Pe become
(103)e¯2Cxxs+1+Czzs+1−e¯1usCxs+1−e¯1wsCzs+1+ϱ3+ϱ4Cs+1−Cts+1=−ζ¯Txxs+Tzzs−ϱ5.Now, applying the Chebyshev spectral collocation to equation (103), we obtain
(104)e¯2I⟶Nt+1⊗I⟶Nz+1⊗1hx2D⟶x2+I⟶Nt+1⊗1hz2D⟶z2⊗I⟶Nx+1−diage¯1UsI⟶Nt+1⊗I⟶Nz+1⊗1hxD⟶x−diage¯1WsI⟶Nt+1⊗1hzD⟶z⊗I⟶Nx−1+diagϱ3+ϱ4−1htD⟶t⊗I⟶Nz+1⊗I⟶Nx+1Cs+1=ζ¯I⟶Nt+1⊗I⟶Nz+1⊗1hx2D⟶x2+I⟶Nt+1⊗1hz2D⟶z2⊗I⟶Nx+1Ts−diagϱ5,where hx, hz, ht, I⟶Nt+1, I⟶Nx+1, and I⟶Nz+1 are all as previously described. Taking out known values of initial and boundary conditions, we obtain
(105)e¯2I⟶Nt⊗I⟶Nz−1⊗1hx2D⟶x22:Nx,:+I⟶Nt⊗1hz2D⟶z22:Nz,:⊗I⟶Nx−1−diage¯1UsI⟶Nt⊗I⟶Nz−1⊗1hxD⟶x2:Nx,:−diage¯1WsI⟶Nt⊗1hzD⟶z2:Nz,:⊗I⟶Nx−1+diagϱ3+ϱ4−1htD⟶t1:Nt,:⊗I⟶Nz−1⊗I⟶Nx−1C~s+1=−ζ¯T~xxs+T~zzs−ψijk5+ψ~ijk5−σijk5+σ~ijk5+δijk5−diagϱ5,where C~=C2:Nx,2:Nz,1:Nt, U~s≈U~s+1, W~s≈W~s+1, and
(106a)ψijk5=e¯2D⟶x22:Nx,1⊗C1,2:Nz,1:Nt−e¯1U~sD⟶x2:Nx,1⊗C1,2:Nz,1:Nt,(106b)ψ~ijk5=e¯2D⟶x22:Nx,Nx+1⊗CNx+1,2:Nz,1:Nt−e¯1U~sD⟶x2:Nx,Nx+1⊗CNx+1,2:Nz,1:Nt,(106c)σijk5=e¯2D⟶z22:Nz,1⊗C2:Nx,1,1:Nt−e¯1W~sD⟶z2:Nz,1⊗C2:Nx,1,1:Nt,(106d)σ~ijk5=e¯2D⟶z22:Nz,Nz+1⊗C2:Nz,Nz+1,1:Nt−e¯1W~sD⟶z2:Nz,Nz+1⊗C2:Nx,Nz+1,1:Nt,(106e)δijk5=−D⟶t1:Nt,Nt+1⊗C2:Nx,2:Nz,Nt+1.Thus, at step (e), we solve forC~s+1 in the matrix equation
(107)A⟶5C~s+1=B⟶5+K⟶5,where
(108)A⟶5=e¯2I⟶Nt⊗I⟶Nz−1⊗1hx2D⟶x22:Nx,:+I⟶Nt⊗1hz2D⟶z22:Nz,:⊗I⟶Nx−1−diage¯1UsI⟶Nt⊗I⟶Nz−1⊗1hxD⟶x2:Nx,:−diage¯1WsI⟶Nt⊗1hzD⟶z2:Nz,:⊗I⟶Nx−1+diagϱ3+ϱ4−1htD⟶t1:Nt,:⊗I⟶Nz−1⊗I⟶Nx−1,and
(109a)K⟶5=−ψijk5+ψ~ijk5−σijk7+σ~ijk7+δijk5−diagρ5,(109b)B⟶5=−ζ¯T~xxs+T~zzs.The spectral-collocated field variables at the initial and boundary points are then given by(110a)U~⟶s+1:,:,1=0,T~s+1:,:,1=Ti,Θ~s+1:,:,1=θ0,C~s+1:,:,1=Ci,(110b)U~s+1Nx+1,:,:=UwmU0ϕ,W~s+11,:,:=0,W~s+1Nx+1,:,:=0,U~s+11,:,:=UwmU0ϕ,(110c)U~s+1:,Nz+1,:=0,U~s+1:,1,:=0,W~s+1:,Nz+1,:=0,W~s+1:,1,:=VwmU02ϕ,(110d)I⟶Nt+1⊗I⟶Nz+1⊗D⟶Nx+1,:,:1T~s+1=q′dT0−Ti,I⟶Nt+1⊗I⟶Nz+1⊗D⟶1,:,:1T~s+1=q′dT0−Ti,(110e)T~s+1:,Nz+1,:=−TiT0−Ti,T~s+1:,1,:=T0−2TiT0−Ti,C~s+1:,Nz+1,:=0,(110f)I⟶Nt+1⊗I⟶Nz+1⊗D⟶Nx+1,:,:1C~s+1=0,I⟶Nt+1⊗I⟶Nz+1⊗D⟶1,:,:1C~s+1=0,C~s+1:,1,:=1,(110g)I⟶Nt+1⊗I⟶Nz+1⊗D⟶Nx+1,:,:1Θ~s+1=Exdθs−θr,I⟶Nt+1⊗I⟶Nz+1⊗D⟶1,:,:1Θ~s+1=Exdθs−θr,(110h)Θ~s+1:,Nz+1,:=0,Θ~s+1:,1,:=θ0−θrθs−θr.Then, starting from initial guesses, we obtain approximate solutions forux,z,t, wx,z,t, Tx,z,t, θx,z,t, and Cx,z,t by iteratively solving the matrix equations (73), (81), (91), (99), and (107), in turn, for s=0,1,2,⋯. We only confirm convergence to the expected solutions by measuring the norm residual errors of the difference between successive iterations as defined in (64).


## 4. Numerical Simulations and Discussion

In this section, numerical results of the stability and convergence analyses of the proposed method over the given model are presented and discussed. A convergence test using residual errors is conducted to examine the accuracy of the scheme. To further validate its performance, numerical computations of the various flow profiles of the model are presented to demonstrate applicability of the proposed MV-SLQLM. In computing the numerical results, we use the model parameter values presented in Tables1 and 2. All computations were performed on an Intel core i7 processor of 8th Generation using MATLAB R2015a.Table 1
Dimensionless model numbers and parameter values.
SymbolDescriptionValue rangeReference(s)NumbersRepPores’ Reynolds number, laminar Darcy flows10−2−101[40]D′aDarcy number3.33×10−6−3.33×10−3[92]F′Forchheimer numbere∗[40]GrCConcentration Grashof number1.12−4.2[93]GrTTemperature Grashof number0.67−3.31[93]PrPrandtl number1.0—15[94]EcEckert number1.8[74, 95]PePeclet number0.05−5.0[23, 96]DaDamko¨hler number1.0−11.5[23]Le′Lewis numbere∗Eq. (19)Parametersϱ1Heat capacity ratio0.01[95]ϱ2Characteristic timee∗L/U0δ1Heat generation/absorption2.0[95]ϱ′3Apparent Danko¨hler numbere∗equation (27)ϱ′4Weighted Danko¨hler numbere∗equation (29)ϱ′5Apparent chemical sinke∗equation (29)RmMatrix retardation factore∗equation (13)δ2Apparent diffusion ratioe∗equation (19)RnMoisture retardation factore∗equation (19)Note:e∗ means that the value is empirically estimated. The empirical values are computed at θr=0.124m3m−3, θs=0.495m3m−3, C0=0.125Mgm−3, Ci=0.1Mgm−3, and Cs=0.211Mgm−3.Table 2
Parameter values for loamy sand soil considered for laminar flow.
Thermodynamic, transport, and soil propertiesParameterDescriptionSI unitValue (s)SourceKhSoil thermal conductivityJs−1∘C−145.318×θ∗0.6619[39]ϕPorosity−0−1[36]kPermeabilitym210−20−10−5[40]KθHydraulic conductivityms−110−7−10−3[41]ζsWater uptake/sinkday−10.013−0.386[97][22]DoMolecular diffusion of NaCl in waterm2s−11×10−9−2×10−9[18][31, 39]KsSaturated hydraulic conductivityms−11.23×10−7−6.67×10−6[1][97][98]DxSoil thermal diffusivitym2s−15.78×10−9−4.44×10−7[1]θrResidual moisture content−0.005−0.027[31, 39]θsSaturated moisture content−0.350−0.396[31, 39][97, 98]θwMedium water content−0.029−0.454[30][99]FForchheimer coefficient−0.0−0.5[40]ρbAverage bulk densityMgm−31.470[31, 39]LoLatent heat of vaporizationJkg−12.5×106−2369.2T∗[31, 100]DTMoisture diffusivity under temperature gradientm2s−1∘C−13.6×10−14.012×T∗6.9146[39]κdDistribution coefficientL3Mg−10.013−0.086[22, 101]Note:∗-empirically estimated at T∗=25∘C; θ∗=θw.


### 4.1. Stability and Rate of Convergence

At varying localPe, changes in the solute concentration profile for the given grid numbers are shown in Figure 2.Figure 2
Effect of varying localPe on the distribution of solute concentration profile at different collocation points.
(a)
0<Pe≤3, Nx=Nz=20, and Nt=10(b)
Pe>2, Nx=Nz=20, and Nt=5(c)
0<Pe≤3, Nx=Nz=25, and Nt=10(d)
Pe>2, Nx=Nz=25, and Nt=5From Figure2, stability of the numerical method is more pronounced when solving equation (13) with local Pe>2, particularly for the collocation points Nx=Nz=20 and Nt=5 (see Figure 2(b)). The reason is that, unlike in Figure 2(b), the convection-diffusion-reaction numerical solution is more sensitive to the data perturbations in Figures 2(a), 2(c), and 2(d). Particularly, we notice from Figure 2 that stability deteriorates as the number of spatial collocation points increases (see Figures 2(c) and 2(d)). Further, for successively increased number of spatial collocation points, Figures 3–7 show the infinity error norms at different iteration points for each equation in the model.Figure 3
Error norms at different iteration points of the model, forNx=Nz=5 and Nt=5.
(a)
u-velocity equation(b)
w-velocity equation(c)
Heat equation(d)
Moisture equation(e)
Solute equationFigure 4
Error norms at different iteration points of the model, forNx=Nz=10 and Nt=5.
(a)
u-velocity equation(b)
w-velocity equation(c)
Heat equation(d)
Moisture equation(e)
Solute equationFigure 5
Error norms at different iteration points of the model, forNx=Nz=15 and Nt=5.
(a)
u-velocity equation(b)
w-velocity equation(c)
Heat equation(d)
Moisture equation(e)
Solute equationFigure 6
Error norms at different iteration points of the model, forNx=Nz=20 and Nt=5.
(a)
u-velocity equation(b)
w-velocity equation(c)
Heat equation(d)
Moisture equation(e)
Solute equationFigure 7
Error norms at different iteration points of the model, forNx=Nz=25 and Nt=5.
(a)
u-velocity equation(b)
w-velocity equation(c)
Heat equation(d)
Moisture equation(e)
Solute equationSince the present scheme is limited to short temporal domains, the number of collocation points in time was fixed toNt=5. Nonetheless, overall, the results show that solution errors introduced in one iteration step do not accumulate in the successive iteration points, affirming general stability of the MV-SLQLM, perhaps as attributed to the optimally stable property of Lagrange interpolation basis function also [88, 89]. Moreover, despite using small values of Nx, Nz, and Nt (see Figures 3 and 4), we obtain highly accurate solutions to the model using the present method. However, from the inspection of Figure 6, the Nx=Nz=20 and Nt=5 collocation points are sufficient to give accurate and consistent results. The reason is that, unlike at 20×20 spatial collocation points, jumps and oscillating wiggles appear at the right end of iteration domain in the solute equation for 15×15 (see Figure 5) and 25×25 (see Figure 7) collocation points, respectively. The excellent superexponential convergence of residual errors to almost machine precision level (Figure 6) is an indication that the proposed MV-SLQLM is reliable to model the distribution and interaction of heat and mass transfer in the porous media.In Table3, the numerical results based on running/CPU time and condition number of the coefficient matrices, obtained using the proposed method, are compared with those presented by the TV-SCM [84] and MV-SQLM [72] at two subinterval lengths 0,T in t variable and using the same number of collocation points in the spatial and time domains.Table 3
Comparing running time and matrix coefficient condition number of the proposed method and the TV-SCM [84] and MV-SQLM [72] for problem, equations (10)–(13) and (19); Nx=Nz=20, Nt=5, and It=20.
TMV-SLQLMTV-SCM [84]MV-SQLM [72]0.51.00.51.00.51.0CPU time (sec)0.210530.237411.324471.4632820.484650.45792Condition number1.1804×1041.3548×1043.34237×1063.43821×1062.4567×1042.5045×104From the inspection of Table3, we observe that for each method, both the computational time and condition number increase with time, T but with small values recorded for the MV-SLQLM. The present results suggest that employing the LLM to decouple the quasilinearized equations saves more on computational time than solving it as a compact matrix system of quasilinearized equations (see [72, 84]). Furthermore, in agreement with Hidayat et al. [49], Table 3 shows that low computational effort of collocation is indeed guaranteed for schemes with Kronecker delta properties. This is evidenced from a comparison of MV-SLQLM and MV-SQLM in terms of CPU times, which both use Kronecker delta properties, suggesting that the current proposed method (MV-SLQLM) is more efficient than the TV-SCM [84] and MV-SQLM [72] at solving large coupled systems of PDEs.


### 4.2. Validation

The local linearization method (LLM) described by equations (41a)–(41e) and the subsequent steps (a)–(e) are then used to generate approximate numerical solutions for the system of governing equations (10)–(13). The linearized equations (73), (81), (91), (99), and (107) are then subsequently solved using the MV-SLQLM. The numerical values of pertinent parameters are kept fixed as ϕ=0.62, F=0.5, Λ=1.3, GrC=2.02, GrT=1.07, Rn=1.0, Le=0.31, Ec=1.8, Rep=1.05, Da=2.5, Pr=0.2, Pe=0.25, κd=0.013, and ζs=0.257, otherwise stated. The numerical simulations are discussed in comparison with published results in the literature. In addition, Nx=Nz=20 and Nt=5 collocation points in the space and time domain, respectively, sufficient to give accurate and consistent results are used.


#### 4.2.1. Velocity Distribution

Figure8 illustrates the variation of velocity profile against depth z for the vertical preferential flow at changes in ϕ, Rep, GrT, GrC, F, and Λ, while the other parameters keep constant and read as described above.Figure 8
Effect ofϕ, Rep, GrT, GrC, F, and Λ on vertical velocity flow profile.
(a)
Varyingϕ(b)
VaryingRep(c)
VaryingGrT(d)
VaryingGrC(e)
VaryingF(f)
VaryingΛAs expected, velocity reduces for greater values of porosity and Forchheimer number due to the increase in viscosity. The effect of varying porosityϕ and Forchheimer number on velocity profile is shown in Figures 8(a) and 8(e), respectively. As it is found, for higher values of ϕ or Forchheimer number, indeed velocity reduces. Physically, the porous medium characterizes one of the types of external resistance forces that work to resist the movement of the fluid through it, and thus, the velocity of the fluid decreases. Figure 8(f) shows the effect of varying the Darcy number Λ on the velocity distribution W,ms−1. Literature retains that increasing the values of Darcy number changes the rate of fluid-permeability through the porous medium, thereby leading to a reduced velocity distribution. However, through the present study, the velocity distribution of the fluid continually increases with the Darcy number Λ enhancement, consistent with [35]. Perhaps, the reason is that at low values Λ, as caused by a reduction in the medium’s permeability k, there is an increase in flow resistance. It is therefore an increase in capillary pressure, due to flow resistance, which causes a reduction in velocity of flow at low Λ numbers. For higher values of Grashof number Gr, however, velocity distribution enhances because of the decrease in viscous forces. In agreement, Figures 8(c) and 8(d) show that the velocity distribution increases with increase in GrT and GrC.Unlike in porous channel (single-phase flow) studies, Figure8(b) shows that the velocity profile decreases with an increase in Rep. The reason is that the velocity distribution decreases to satisfy the conservation of pore momentum as Rep increases because of the condensation process (Figure 1(b)) that occurs at a low Rep number, consistent with [102]. In addition, in agreement with [102], Figure 8(b) shows that the velocity decay rate is slower at high Rep than at low Rep. Moreover, for liquid flows in porous media such as soil, the Reynolds number is generally low, leading to more occurrence of the condensation process.


#### 4.2.2. Heat Distribution

In [92], it was claimed that the local thermal equilibrium in a porous medium becomes stronger as anyone of Rep, Pr, or the Λ decreases, or the Nusselt number Nu increases. Therefore, we also investigated the effect of Pr and Λ in addition to Ec on the soil temperature distribution profile. Figure 9 shows the variation in temperature profile for varied values of the Prandtl number, Reynolds number, and heat generation/absorption parameter.Figure 9
Effect ofRep, Pr, and δ1 on soil temperature profile.
(a)
VaryingRep(b)
VaryingPr(c)
Varyingδ1Figure9(a) shows that the heat distribution profile advances at greater values of Rep, consistent with [102]. When [94] studied the effect of porosity on natural convection of moderate Prandtl number fluids (Pr≥1) in a porous medium, heat transfer rate increased with porosity. The reason is that higher porosity leads to a larger contact surface area, thereby intensifying the heat transfer process. However, even at the moderate Prandtl number of 5.2, the heat transfer profiles insignificantly changes at increasing values of ϕ. We therefore do not include its result here. Despite the case, the constant increased values of heat transfer at higher porosity values could be that, for moderately high Prandtl number fluids in porous media, Nu increases nearly linearly with ϕ [94]. The linear relationship between Nu and ϕ, as reported in [94], could explain the observed insignificant changes but at high-temperature values.The Prandtl numberPr is the ratio of momentum to thermal diffusivity. Figure 9(b) shows the effect of Pr on the temperature distribution profile. Since Pr is inversely proportional to thermal conductivity, several studies indicate that enhancing the Prandtl number causes a decline in the thermal profile because high estimates of the Prandtl number mean weaker thermal diffusivity. Hence, heat diffuses gradually and thus affecting the temperature. However, in the current numerical simulations, soil temperature profiles increase with Prandtl number because the heat transfer characteristics are now influenced by several interactive variables inside the natural convection-dominated porous zone (Figure 1(b)). For heat and mass transfer characteristics of natural convection of large Prandtl number 1≤Pr≤200 fluids in porous media such as soil, a recent study found that there exist three characteristic heat transfer zones where the Nusselt number Nu (ratio of heat convection to conduction) depends on Pr monotonically [94]. In particular, [94] found that Nu increases with Pr monotonically for moderate 1≤Pr≤30 Prandtl number fluids but decreases in the 30<Pr≤140 zone and slightly increases again for Pr>140. Therefore, with reference to the zone of Pr values (Figure 9(b)), soil heat increases with Pr due to the enhanced Nu for natural convection flow [92].Eckert numberEc characterizes the influence of self-heating of fluid as a consequence of dissipation effects due to internal friction of the fluid. However, in the current study, varying values of Ec did not show any significant effect on the temperature profiles. Therefore, in the current investigated porous domain (Figure 1(b)), the effects of self-heating due to dissipation are negligible. However, under effect of the heat generation/absorption parameter δ1, the temperature profile increases (Figure 9(c)), where δ>0 means a heat generation and δ1<0 means a heat absorption, because heat generation escalates the spread and distribution of the temperature. With heat absorption, on the other hand, the temperature of the fluid begins to decline. Hence, the response to the negative temperature gradient from the heat absorption to the heat generation mechanism results in a high thermal diffusion that makes the temperature of the boundary layer rise, leading to a rise in the temperature distribution of the fluid. Finally, we also investigated the effect of varying values of Darcy numbers Λ on temperature profiles. Unfortunately, like a case of porosity ϕ that monotonically depends on permeability k, temperature distribution profiles do not significantly change with Darcy numbers.


#### 4.2.3. Moisture Distribution

As demonstrated by several studies, changes in retardation, temperature gradient, diffusivity ratio, and sink can alter moisture transfer characteristics through the porous medium. Figure10 shows the effect of varying ϕ, Rn, δ2, and Le′ on changes to moisture distribution. We chose these parameters because, elsewhere, they showed to be sensitive to calculating moisture content.Figure 10
Effect ofϕ, Rn, δ2, and Le′ on distribution of moisture.
(a)
Varyingϕ(b)
VaryingRn(c)
Varyingδ2(d)
VaryingLe′At varying values ofϕ, moisture distribution profile does not significantly change with porosity (Figure 10(a)). The result correlates with nonsignificant changes to temperature profiles at varying ϕ since increasing porosity enhances evaporation rate, which may lead to reduced moisture content. But, at low Prandtl and Eckert numbers (minimal heat effect), Figure 10(b) shows that soil moisture profiles decrease with an increase in retardation Rn, attributed to the sink effect when the plant-root uptake mechanism has a significant impact on retardation factor. Consistent with the previous studies, Figure 10(c) shows that moisture profiles decrease with the increasing values of δ2 due to the redistribution process [25, 103] and low values of isothermal moisture diffusivity Dθ. Using Le′=DT/Dθ, Figure 10(d) shows the effect of thermal moisture diffusivity DT on moisture profile in the investigated domain under combined moisture and temperature gradients. Numerical simulations show that for high Le′, the evaporation rate is very marginal, leading to a significant impact on moisture content.


#### 4.2.4. Solute Concentration Distribution

Using analytical or semianalytical solutions, studies have shown that solute concentration profiles get affected by the changes in several model parameters, includingPe, ζ¯, and θw. In addition, [23] observed the changes in the solute transfer characteristics with the chemical reaction parameter by varying the Damko¨hler number Da. Similarly, the effect of reaction and absorption coefficients on solute concentration was investigated, respectively, using the deduced “apparent” Damko¨hler ρ′3 and “weighted” Damko¨hler ρ′4 numbers. However, we specifically observed the effects of varying Da and κd present in ρ′3 and ρ′4, respectively. To investigate the effect of plant-root uptake activity, the effect of varying the sink coefficient ζs was observed. The effect of soil temperatures on solute distribution was investigated by varying the thermophoresis parameter ζ¯. Figure 11 shows the results.Figure 11
Effect ofPe, ζs, ζ¯, Da, κd, and θw on solute concentration profile.
(a)
VaryingPe(b)
Varyingζs(c)
Varyingζ¯(d)
VaryingDa(e)
Varyingκd(f)
VaryingθwFigure11(a) shows the effect of varying low Peclet numbers 0≤Pe≤5 on the solute concentration at low Rep and moderate Da. As mentioned earlier on, we chose the low range of Pe numbers because for the investigated domain that is characterized by laminar flow, diffusion often dominates the transport in the direction perpendicular to the flow streamlines. According to Figure 11(a), solute concentration is generally low at low Peclet numbers with more impact from small than higher Pe values. The reason is that high Pe values, as characterized by the increase in flow velocity or diameter of pore size, tend to enhance the distribution. Since the “apparent” chemical sink parameter ϱ´5 largely depends on ζs (see, equation (29)), the effect of ζs on the solute concentration is investigated, and Figure 11(b) shows the results. As expected, solute concentration decreases for an increase in sink parameters, attributed to the plant-root uptake process.Regarding the change that occurs in the concentration of solute distribution, in the case of the effect of the thermophoresis parameterζ¯, Figure 11(c) shows that escalation in the values of ζ¯ leads to a higher concentration of the solute distribution in the natural convection zone. The reason is that with an increase in the thermophoresis coefficient, transport force due to the presence of temperature gradient increases, leading to the transfer of more solute particles from the boundary layer region to the cooler region. Hence, by increasing ζ¯, the concentration in the boundary layer decreases, thereby enhancing that of the natural convection zone (cool region). Dejam [23] reported that an increase in Damko¨hler numbers Da decreases solute concentration because the increased reaction rate leads to a generation of more species. To the contrary, the results of Figure 11(d) show that the distribution of solute concentration increases with Da. The inconsistency arises due to differences in the investigated domain where [23] studied solute transfer in a fractured porous medium.Last but not least, the role of distribution coefficientκd in the effect of “weighted” Damko¨hler number ϱ´4 (see equation (29)) on the distribution profiles of solute concentration is shown in Figure 11(e). Other studies found that an increase in κd decreases solute concentration because increasing κd influences the decaying process. In constrast, Figure 11(e) suggests that when both heat and water contents are present in the media, κd insignificantly affects changes of the solute distribution. The reason is that an increase in water content or heat somewhat promotes contact areas, leading to a high rate of solute mixing [104, 105]. Finally, Figure 11(f) shows the effect of varying volumetric water content θw on the distribution of solute concentration in which the distribution is low at higher average values of θw. The result suggests that the distribution mechanism is now largely facilitated by the increase in dispersion coefficient as the average volumetric water content θw reduces [105]. However, it is important to point out that these results are in contrast with other studies. This difference can be attributed to the behavior of molecular diffusion Dw under saturation. Unlike in the unsaturated flow, Dw contributes greatly to the dispersion coefficient because high water contents enhance solute mixing so fast so that the media reaches the Fickian region (equilibrium condition) faster.


## 4.1. Stability and Rate of Convergence

At varying localPe, changes in the solute concentration profile for the given grid numbers are shown in Figure 2.Figure 2
Effect of varying localPe on the distribution of solute concentration profile at different collocation points.
(a)
0<Pe≤3, Nx=Nz=20, and Nt=10(b)
Pe>2, Nx=Nz=20, and Nt=5(c)
0<Pe≤3, Nx=Nz=25, and Nt=10(d)
Pe>2, Nx=Nz=25, and Nt=5From Figure2, stability of the numerical method is more pronounced when solving equation (13) with local Pe>2, particularly for the collocation points Nx=Nz=20 and Nt=5 (see Figure 2(b)). The reason is that, unlike in Figure 2(b), the convection-diffusion-reaction numerical solution is more sensitive to the data perturbations in Figures 2(a), 2(c), and 2(d). Particularly, we notice from Figure 2 that stability deteriorates as the number of spatial collocation points increases (see Figures 2(c) and 2(d)). Further, for successively increased number of spatial collocation points, Figures 3–7 show the infinity error norms at different iteration points for each equation in the model.Figure 3
Error norms at different iteration points of the model, forNx=Nz=5 and Nt=5.
(a)
u-velocity equation(b)
w-velocity equation(c)
Heat equation(d)
Moisture equation(e)
Solute equationFigure 4
Error norms at different iteration points of the model, forNx=Nz=10 and Nt=5.
(a)
u-velocity equation(b)
w-velocity equation(c)
Heat equation(d)
Moisture equation(e)
Solute equationFigure 5
Error norms at different iteration points of the model, forNx=Nz=15 and Nt=5.
(a)
u-velocity equation(b)
w-velocity equation(c)
Heat equation(d)
Moisture equation(e)
Solute equationFigure 6
Error norms at different iteration points of the model, forNx=Nz=20 and Nt=5.
(a)
u-velocity equation(b)
w-velocity equation(c)
Heat equation(d)
Moisture equation(e)
Solute equationFigure 7
Error norms at different iteration points of the model, forNx=Nz=25 and Nt=5.
(a)
u-velocity equation(b)
w-velocity equation(c)
Heat equation(d)
Moisture equation(e)
Solute equationSince the present scheme is limited to short temporal domains, the number of collocation points in time was fixed toNt=5. Nonetheless, overall, the results show that solution errors introduced in one iteration step do not accumulate in the successive iteration points, affirming general stability of the MV-SLQLM, perhaps as attributed to the optimally stable property of Lagrange interpolation basis function also [88, 89]. Moreover, despite using small values of Nx, Nz, and Nt (see Figures 3 and 4), we obtain highly accurate solutions to the model using the present method. However, from the inspection of Figure 6, the Nx=Nz=20 and Nt=5 collocation points are sufficient to give accurate and consistent results. The reason is that, unlike at 20×20 spatial collocation points, jumps and oscillating wiggles appear at the right end of iteration domain in the solute equation for 15×15 (see Figure 5) and 25×25 (see Figure 7) collocation points, respectively. The excellent superexponential convergence of residual errors to almost machine precision level (Figure 6) is an indication that the proposed MV-SLQLM is reliable to model the distribution and interaction of heat and mass transfer in the porous media.In Table3, the numerical results based on running/CPU time and condition number of the coefficient matrices, obtained using the proposed method, are compared with those presented by the TV-SCM [84] and MV-SQLM [72] at two subinterval lengths 0,T in t variable and using the same number of collocation points in the spatial and time domains.Table 3
Comparing running time and matrix coefficient condition number of the proposed method and the TV-SCM [84] and MV-SQLM [72] for problem, equations (10)–(13) and (19); Nx=Nz=20, Nt=5, and It=20.
TMV-SLQLMTV-SCM [84]MV-SQLM [72]0.51.00.51.00.51.0CPU time (sec)0.210530.237411.324471.4632820.484650.45792Condition number1.1804×1041.3548×1043.34237×1063.43821×1062.4567×1042.5045×104From the inspection of Table3, we observe that for each method, both the computational time and condition number increase with time, T but with small values recorded for the MV-SLQLM. The present results suggest that employing the LLM to decouple the quasilinearized equations saves more on computational time than solving it as a compact matrix system of quasilinearized equations (see [72, 84]). Furthermore, in agreement with Hidayat et al. [49], Table 3 shows that low computational effort of collocation is indeed guaranteed for schemes with Kronecker delta properties. This is evidenced from a comparison of MV-SLQLM and MV-SQLM in terms of CPU times, which both use Kronecker delta properties, suggesting that the current proposed method (MV-SLQLM) is more efficient than the TV-SCM [84] and MV-SQLM [72] at solving large coupled systems of PDEs.


## 4.2. Validation

The local linearization method (LLM) described by equations (41a)–(41e) and the subsequent steps (a)–(e) are then used to generate approximate numerical solutions for the system of governing equations (10)–(13). The linearized equations (73), (81), (91), (99), and (107) are then subsequently solved using the MV-SLQLM. The numerical values of pertinent parameters are kept fixed as ϕ=0.62, F=0.5, Λ=1.3, GrC=2.02, GrT=1.07, Rn=1.0, Le=0.31, Ec=1.8, Rep=1.05, Da=2.5, Pr=0.2, Pe=0.25, κd=0.013, and ζs=0.257, otherwise stated. The numerical simulations are discussed in comparison with published results in the literature. In addition, Nx=Nz=20 and Nt=5 collocation points in the space and time domain, respectively, sufficient to give accurate and consistent results are used.


### 4.2.1. Velocity Distribution

Figure8 illustrates the variation of velocity profile against depth z for the vertical preferential flow at changes in ϕ, Rep, GrT, GrC, F, and Λ, while the other parameters keep constant and read as described above.Figure 8
Effect ofϕ, Rep, GrT, GrC, F, and Λ on vertical velocity flow profile.
(a)
Varyingϕ(b)
VaryingRep(c)
VaryingGrT(d)
VaryingGrC(e)
VaryingF(f)
VaryingΛAs expected, velocity reduces for greater values of porosity and Forchheimer number due to the increase in viscosity. The effect of varying porosityϕ and Forchheimer number on velocity profile is shown in Figures 8(a) and 8(e), respectively. As it is found, for higher values of ϕ or Forchheimer number, indeed velocity reduces. Physically, the porous medium characterizes one of the types of external resistance forces that work to resist the movement of the fluid through it, and thus, the velocity of the fluid decreases. Figure 8(f) shows the effect of varying the Darcy number Λ on the velocity distribution W,ms−1. Literature retains that increasing the values of Darcy number changes the rate of fluid-permeability through the porous medium, thereby leading to a reduced velocity distribution. However, through the present study, the velocity distribution of the fluid continually increases with the Darcy number Λ enhancement, consistent with [35]. Perhaps, the reason is that at low values Λ, as caused by a reduction in the medium’s permeability k, there is an increase in flow resistance. It is therefore an increase in capillary pressure, due to flow resistance, which causes a reduction in velocity of flow at low Λ numbers. For higher values of Grashof number Gr, however, velocity distribution enhances because of the decrease in viscous forces. In agreement, Figures 8(c) and 8(d) show that the velocity distribution increases with increase in GrT and GrC.Unlike in porous channel (single-phase flow) studies, Figure8(b) shows that the velocity profile decreases with an increase in Rep. The reason is that the velocity distribution decreases to satisfy the conservation of pore momentum as Rep increases because of the condensation process (Figure 1(b)) that occurs at a low Rep number, consistent with [102]. In addition, in agreement with [102], Figure 8(b) shows that the velocity decay rate is slower at high Rep than at low Rep. Moreover, for liquid flows in porous media such as soil, the Reynolds number is generally low, leading to more occurrence of the condensation process.


### 4.2.2. Heat Distribution

In [92], it was claimed that the local thermal equilibrium in a porous medium becomes stronger as anyone of Rep, Pr, or the Λ decreases, or the Nusselt number Nu increases. Therefore, we also investigated the effect of Pr and Λ in addition to Ec on the soil temperature distribution profile. Figure 9 shows the variation in temperature profile for varied values of the Prandtl number, Reynolds number, and heat generation/absorption parameter.Figure 9
Effect ofRep, Pr, and δ1 on soil temperature profile.
(a)
VaryingRep(b)
VaryingPr(c)
Varyingδ1Figure9(a) shows that the heat distribution profile advances at greater values of Rep, consistent with [102]. When [94] studied the effect of porosity on natural convection of moderate Prandtl number fluids (Pr≥1) in a porous medium, heat transfer rate increased with porosity. The reason is that higher porosity leads to a larger contact surface area, thereby intensifying the heat transfer process. However, even at the moderate Prandtl number of 5.2, the heat transfer profiles insignificantly changes at increasing values of ϕ. We therefore do not include its result here. Despite the case, the constant increased values of heat transfer at higher porosity values could be that, for moderately high Prandtl number fluids in porous media, Nu increases nearly linearly with ϕ [94]. The linear relationship between Nu and ϕ, as reported in [94], could explain the observed insignificant changes but at high-temperature values.The Prandtl numberPr is the ratio of momentum to thermal diffusivity. Figure 9(b) shows the effect of Pr on the temperature distribution profile. Since Pr is inversely proportional to thermal conductivity, several studies indicate that enhancing the Prandtl number causes a decline in the thermal profile because high estimates of the Prandtl number mean weaker thermal diffusivity. Hence, heat diffuses gradually and thus affecting the temperature. However, in the current numerical simulations, soil temperature profiles increase with Prandtl number because the heat transfer characteristics are now influenced by several interactive variables inside the natural convection-dominated porous zone (Figure 1(b)). For heat and mass transfer characteristics of natural convection of large Prandtl number 1≤Pr≤200 fluids in porous media such as soil, a recent study found that there exist three characteristic heat transfer zones where the Nusselt number Nu (ratio of heat convection to conduction) depends on Pr monotonically [94]. In particular, [94] found that Nu increases with Pr monotonically for moderate 1≤Pr≤30 Prandtl number fluids but decreases in the 30<Pr≤140 zone and slightly increases again for Pr>140. Therefore, with reference to the zone of Pr values (Figure 9(b)), soil heat increases with Pr due to the enhanced Nu for natural convection flow [92].Eckert numberEc characterizes the influence of self-heating of fluid as a consequence of dissipation effects due to internal friction of the fluid. However, in the current study, varying values of Ec did not show any significant effect on the temperature profiles. Therefore, in the current investigated porous domain (Figure 1(b)), the effects of self-heating due to dissipation are negligible. However, under effect of the heat generation/absorption parameter δ1, the temperature profile increases (Figure 9(c)), where δ>0 means a heat generation and δ1<0 means a heat absorption, because heat generation escalates the spread and distribution of the temperature. With heat absorption, on the other hand, the temperature of the fluid begins to decline. Hence, the response to the negative temperature gradient from the heat absorption to the heat generation mechanism results in a high thermal diffusion that makes the temperature of the boundary layer rise, leading to a rise in the temperature distribution of the fluid. Finally, we also investigated the effect of varying values of Darcy numbers Λ on temperature profiles. Unfortunately, like a case of porosity ϕ that monotonically depends on permeability k, temperature distribution profiles do not significantly change with Darcy numbers.


### 4.2.3. Moisture Distribution

As demonstrated by several studies, changes in retardation, temperature gradient, diffusivity ratio, and sink can alter moisture transfer characteristics through the porous medium. Figure10 shows the effect of varying ϕ, Rn, δ2, and Le′ on changes to moisture distribution. We chose these parameters because, elsewhere, they showed to be sensitive to calculating moisture content.Figure 10
Effect ofϕ, Rn, δ2, and Le′ on distribution of moisture.
(a)
Varyingϕ(b)
VaryingRn(c)
Varyingδ2(d)
VaryingLe′At varying values ofϕ, moisture distribution profile does not significantly change with porosity (Figure 10(a)). The result correlates with nonsignificant changes to temperature profiles at varying ϕ since increasing porosity enhances evaporation rate, which may lead to reduced moisture content. But, at low Prandtl and Eckert numbers (minimal heat effect), Figure 10(b) shows that soil moisture profiles decrease with an increase in retardation Rn, attributed to the sink effect when the plant-root uptake mechanism has a significant impact on retardation factor. Consistent with the previous studies, Figure 10(c) shows that moisture profiles decrease with the increasing values of δ2 due to the redistribution process [25, 103] and low values of isothermal moisture diffusivity Dθ. Using Le′=DT/Dθ, Figure 10(d) shows the effect of thermal moisture diffusivity DT on moisture profile in the investigated domain under combined moisture and temperature gradients. Numerical simulations show that for high Le′, the evaporation rate is very marginal, leading to a significant impact on moisture content.


### 4.2.4. Solute Concentration Distribution

Using analytical or semianalytical solutions, studies have shown that solute concentration profiles get affected by the changes in several model parameters, includingPe, ζ¯, and θw. In addition, [23] observed the changes in the solute transfer characteristics with the chemical reaction parameter by varying the Damko¨hler number Da. Similarly, the effect of reaction and absorption coefficients on solute concentration was investigated, respectively, using the deduced “apparent” Damko¨hler ρ′3 and “weighted” Damko¨hler ρ′4 numbers. However, we specifically observed the effects of varying Da and κd present in ρ′3 and ρ′4, respectively. To investigate the effect of plant-root uptake activity, the effect of varying the sink coefficient ζs was observed. The effect of soil temperatures on solute distribution was investigated by varying the thermophoresis parameter ζ¯. Figure 11 shows the results.Figure 11
Effect ofPe, ζs, ζ¯, Da, κd, and θw on solute concentration profile.
(a)
VaryingPe(b)
Varyingζs(c)
Varyingζ¯(d)
VaryingDa(e)
Varyingκd(f)
VaryingθwFigure11(a) shows the effect of varying low Peclet numbers 0≤Pe≤5 on the solute concentration at low Rep and moderate Da. As mentioned earlier on, we chose the low range of Pe numbers because for the investigated domain that is characterized by laminar flow, diffusion often dominates the transport in the direction perpendicular to the flow streamlines. According to Figure 11(a), solute concentration is generally low at low Peclet numbers with more impact from small than higher Pe values. The reason is that high Pe values, as characterized by the increase in flow velocity or diameter of pore size, tend to enhance the distribution. Since the “apparent” chemical sink parameter ϱ´5 largely depends on ζs (see, equation (29)), the effect of ζs on the solute concentration is investigated, and Figure 11(b) shows the results. As expected, solute concentration decreases for an increase in sink parameters, attributed to the plant-root uptake process.Regarding the change that occurs in the concentration of solute distribution, in the case of the effect of the thermophoresis parameterζ¯, Figure 11(c) shows that escalation in the values of ζ¯ leads to a higher concentration of the solute distribution in the natural convection zone. The reason is that with an increase in the thermophoresis coefficient, transport force due to the presence of temperature gradient increases, leading to the transfer of more solute particles from the boundary layer region to the cooler region. Hence, by increasing ζ¯, the concentration in the boundary layer decreases, thereby enhancing that of the natural convection zone (cool region). Dejam [23] reported that an increase in Damko¨hler numbers Da decreases solute concentration because the increased reaction rate leads to a generation of more species. To the contrary, the results of Figure 11(d) show that the distribution of solute concentration increases with Da. The inconsistency arises due to differences in the investigated domain where [23] studied solute transfer in a fractured porous medium.Last but not least, the role of distribution coefficientκd in the effect of “weighted” Damko¨hler number ϱ´4 (see equation (29)) on the distribution profiles of solute concentration is shown in Figure 11(e). Other studies found that an increase in κd decreases solute concentration because increasing κd influences the decaying process. In constrast, Figure 11(e) suggests that when both heat and water contents are present in the media, κd insignificantly affects changes of the solute distribution. The reason is that an increase in water content or heat somewhat promotes contact areas, leading to a high rate of solute mixing [104, 105]. Finally, Figure 11(f) shows the effect of varying volumetric water content θw on the distribution of solute concentration in which the distribution is low at higher average values of θw. The result suggests that the distribution mechanism is now largely facilitated by the increase in dispersion coefficient as the average volumetric water content θw reduces [105]. However, it is important to point out that these results are in contrast with other studies. This difference can be attributed to the behavior of molecular diffusion Dw under saturation. Unlike in the unsaturated flow, Dw contributes greatly to the dispersion coefficient because high water contents enhance solute mixing so fast so that the media reaches the Fickian region (equilibrium condition) faster.


## 4.2.1. Velocity Distribution

Figure8 illustrates the variation of velocity profile against depth z for the vertical preferential flow at changes in ϕ, Rep, GrT, GrC, F, and Λ, while the other parameters keep constant and read as described above.Figure 8
Effect ofϕ, Rep, GrT, GrC, F, and Λ on vertical velocity flow profile.
(a)
Varyingϕ(b)
VaryingRep(c)
VaryingGrT(d)
VaryingGrC(e)
VaryingF(f)
VaryingΛAs expected, velocity reduces for greater values of porosity and Forchheimer number due to the increase in viscosity. The effect of varying porosityϕ and Forchheimer number on velocity profile is shown in Figures 8(a) and 8(e), respectively. As it is found, for higher values of ϕ or Forchheimer number, indeed velocity reduces. Physically, the porous medium characterizes one of the types of external resistance forces that work to resist the movement of the fluid through it, and thus, the velocity of the fluid decreases. Figure 8(f) shows the effect of varying the Darcy number Λ on the velocity distribution W,ms−1. Literature retains that increasing the values of Darcy number changes the rate of fluid-permeability through the porous medium, thereby leading to a reduced velocity distribution. However, through the present study, the velocity distribution of the fluid continually increases with the Darcy number Λ enhancement, consistent with [35]. Perhaps, the reason is that at low values Λ, as caused by a reduction in the medium’s permeability k, there is an increase in flow resistance. It is therefore an increase in capillary pressure, due to flow resistance, which causes a reduction in velocity of flow at low Λ numbers. For higher values of Grashof number Gr, however, velocity distribution enhances because of the decrease in viscous forces. In agreement, Figures 8(c) and 8(d) show that the velocity distribution increases with increase in GrT and GrC.Unlike in porous channel (single-phase flow) studies, Figure8(b) shows that the velocity profile decreases with an increase in Rep. The reason is that the velocity distribution decreases to satisfy the conservation of pore momentum as Rep increases because of the condensation process (Figure 1(b)) that occurs at a low Rep number, consistent with [102]. In addition, in agreement with [102], Figure 8(b) shows that the velocity decay rate is slower at high Rep than at low Rep. Moreover, for liquid flows in porous media such as soil, the Reynolds number is generally low, leading to more occurrence of the condensation process.


## 4.2.2. Heat Distribution

In [92], it was claimed that the local thermal equilibrium in a porous medium becomes stronger as anyone of Rep, Pr, or the Λ decreases, or the Nusselt number Nu increases. Therefore, we also investigated the effect of Pr and Λ in addition to Ec on the soil temperature distribution profile. Figure 9 shows the variation in temperature profile for varied values of the Prandtl number, Reynolds number, and heat generation/absorption parameter.Figure 9
Effect ofRep, Pr, and δ1 on soil temperature profile.
(a)
VaryingRep(b)
VaryingPr(c)
Varyingδ1Figure9(a) shows that the heat distribution profile advances at greater values of Rep, consistent with [102]. When [94] studied the effect of porosity on natural convection of moderate Prandtl number fluids (Pr≥1) in a porous medium, heat transfer rate increased with porosity. The reason is that higher porosity leads to a larger contact surface area, thereby intensifying the heat transfer process. However, even at the moderate Prandtl number of 5.2, the heat transfer profiles insignificantly changes at increasing values of ϕ. We therefore do not include its result here. Despite the case, the constant increased values of heat transfer at higher porosity values could be that, for moderately high Prandtl number fluids in porous media, Nu increases nearly linearly with ϕ [94]. The linear relationship between Nu and ϕ, as reported in [94], could explain the observed insignificant changes but at high-temperature values.The Prandtl numberPr is the ratio of momentum to thermal diffusivity. Figure 9(b) shows the effect of Pr on the temperature distribution profile. Since Pr is inversely proportional to thermal conductivity, several studies indicate that enhancing the Prandtl number causes a decline in the thermal profile because high estimates of the Prandtl number mean weaker thermal diffusivity. Hence, heat diffuses gradually and thus affecting the temperature. However, in the current numerical simulations, soil temperature profiles increase with Prandtl number because the heat transfer characteristics are now influenced by several interactive variables inside the natural convection-dominated porous zone (Figure 1(b)). For heat and mass transfer characteristics of natural convection of large Prandtl number 1≤Pr≤200 fluids in porous media such as soil, a recent study found that there exist three characteristic heat transfer zones where the Nusselt number Nu (ratio of heat convection to conduction) depends on Pr monotonically [94]. In particular, [94] found that Nu increases with Pr monotonically for moderate 1≤Pr≤30 Prandtl number fluids but decreases in the 30<Pr≤140 zone and slightly increases again for Pr>140. Therefore, with reference to the zone of Pr values (Figure 9(b)), soil heat increases with Pr due to the enhanced Nu for natural convection flow [92].Eckert numberEc characterizes the influence of self-heating of fluid as a consequence of dissipation effects due to internal friction of the fluid. However, in the current study, varying values of Ec did not show any significant effect on the temperature profiles. Therefore, in the current investigated porous domain (Figure 1(b)), the effects of self-heating due to dissipation are negligible. However, under effect of the heat generation/absorption parameter δ1, the temperature profile increases (Figure 9(c)), where δ>0 means a heat generation and δ1<0 means a heat absorption, because heat generation escalates the spread and distribution of the temperature. With heat absorption, on the other hand, the temperature of the fluid begins to decline. Hence, the response to the negative temperature gradient from the heat absorption to the heat generation mechanism results in a high thermal diffusion that makes the temperature of the boundary layer rise, leading to a rise in the temperature distribution of the fluid. Finally, we also investigated the effect of varying values of Darcy numbers Λ on temperature profiles. Unfortunately, like a case of porosity ϕ that monotonically depends on permeability k, temperature distribution profiles do not significantly change with Darcy numbers.


## 4.2.3. Moisture Distribution

As demonstrated by several studies, changes in retardation, temperature gradient, diffusivity ratio, and sink can alter moisture transfer characteristics through the porous medium. Figure10 shows the effect of varying ϕ, Rn, δ2, and Le′ on changes to moisture distribution. We chose these parameters because, elsewhere, they showed to be sensitive to calculating moisture content.Figure 10
Effect ofϕ, Rn, δ2, and Le′ on distribution of moisture.
(a)
Varyingϕ(b)
VaryingRn(c)
Varyingδ2(d)
VaryingLe′At varying values ofϕ, moisture distribution profile does not significantly change with porosity (Figure 10(a)). The result correlates with nonsignificant changes to temperature profiles at varying ϕ since increasing porosity enhances evaporation rate, which may lead to reduced moisture content. But, at low Prandtl and Eckert numbers (minimal heat effect), Figure 10(b) shows that soil moisture profiles decrease with an increase in retardation Rn, attributed to the sink effect when the plant-root uptake mechanism has a significant impact on retardation factor. Consistent with the previous studies, Figure 10(c) shows that moisture profiles decrease with the increasing values of δ2 due to the redistribution process [25, 103] and low values of isothermal moisture diffusivity Dθ. Using Le′=DT/Dθ, Figure 10(d) shows the effect of thermal moisture diffusivity DT on moisture profile in the investigated domain under combined moisture and temperature gradients. Numerical simulations show that for high Le′, the evaporation rate is very marginal, leading to a significant impact on moisture content.


## 4.2.4. Solute Concentration Distribution

Using analytical or semianalytical solutions, studies have shown that solute concentration profiles get affected by the changes in several model parameters, includingPe, ζ¯, and θw. In addition, [23] observed the changes in the solute transfer characteristics with the chemical reaction parameter by varying the Damko¨hler number Da. Similarly, the effect of reaction and absorption coefficients on solute concentration was investigated, respectively, using the deduced “apparent” Damko¨hler ρ′3 and “weighted” Damko¨hler ρ′4 numbers. However, we specifically observed the effects of varying Da and κd present in ρ′3 and ρ′4, respectively. To investigate the effect of plant-root uptake activity, the effect of varying the sink coefficient ζs was observed. The effect of soil temperatures on solute distribution was investigated by varying the thermophoresis parameter ζ¯. Figure 11 shows the results.Figure 11
Effect ofPe, ζs, ζ¯, Da, κd, and θw on solute concentration profile.
(a)
VaryingPe(b)
Varyingζs(c)
Varyingζ¯(d)
VaryingDa(e)
Varyingκd(f)
VaryingθwFigure11(a) shows the effect of varying low Peclet numbers 0≤Pe≤5 on the solute concentration at low Rep and moderate Da. As mentioned earlier on, we chose the low range of Pe numbers because for the investigated domain that is characterized by laminar flow, diffusion often dominates the transport in the direction perpendicular to the flow streamlines. According to Figure 11(a), solute concentration is generally low at low Peclet numbers with more impact from small than higher Pe values. The reason is that high Pe values, as characterized by the increase in flow velocity or diameter of pore size, tend to enhance the distribution. Since the “apparent” chemical sink parameter ϱ´5 largely depends on ζs (see, equation (29)), the effect of ζs on the solute concentration is investigated, and Figure 11(b) shows the results. As expected, solute concentration decreases for an increase in sink parameters, attributed to the plant-root uptake process.Regarding the change that occurs in the concentration of solute distribution, in the case of the effect of the thermophoresis parameterζ¯, Figure 11(c) shows that escalation in the values of ζ¯ leads to a higher concentration of the solute distribution in the natural convection zone. The reason is that with an increase in the thermophoresis coefficient, transport force due to the presence of temperature gradient increases, leading to the transfer of more solute particles from the boundary layer region to the cooler region. Hence, by increasing ζ¯, the concentration in the boundary layer decreases, thereby enhancing that of the natural convection zone (cool region). Dejam [23] reported that an increase in Damko¨hler numbers Da decreases solute concentration because the increased reaction rate leads to a generation of more species. To the contrary, the results of Figure 11(d) show that the distribution of solute concentration increases with Da. The inconsistency arises due to differences in the investigated domain where [23] studied solute transfer in a fractured porous medium.Last but not least, the role of distribution coefficientκd in the effect of “weighted” Damko¨hler number ϱ´4 (see equation (29)) on the distribution profiles of solute concentration is shown in Figure 11(e). Other studies found that an increase in κd decreases solute concentration because increasing κd influences the decaying process. In constrast, Figure 11(e) suggests that when both heat and water contents are present in the media, κd insignificantly affects changes of the solute distribution. The reason is that an increase in water content or heat somewhat promotes contact areas, leading to a high rate of solute mixing [104, 105]. Finally, Figure 11(f) shows the effect of varying volumetric water content θw on the distribution of solute concentration in which the distribution is low at higher average values of θw. The result suggests that the distribution mechanism is now largely facilitated by the increase in dispersion coefficient as the average volumetric water content θw reduces [105]. However, it is important to point out that these results are in contrast with other studies. This difference can be attributed to the behavior of molecular diffusion Dw under saturation. Unlike in the unsaturated flow, Dw contributes greatly to the dispersion coefficient because high water contents enhance solute mixing so fast so that the media reaches the Fickian region (equilibrium condition) faster.


## 5. Conclusions

The purpose of the current study was to propose and apply a multivariate Lagrange polynomial based spectral collocation method to solve a system of coupled nonlinear partial differential equations for modelling transport and interaction of soil moisture, heat, and solute concentration inside a soil ridge when water is present in the furrows. The proposed solution method first uses a quasilinearization method to linearize any nonlinear equations and then employs a local linearization technique to decouple the linearized system of PDEs to form a sequence of equations that are solved in a computationally efficient manner. The applicability of the proposed method, termed the multivariate spectral local quasilinearization method (MV-SLQLM), was demonstrated by the systematic analysis of the effects of varying model parameters on the distribution profiles of velocity, moisture, heat, and solute concentration. The accuracy of the MV-SLQLM solutions on the model was examined through the analysis of residual errors that excellently converge to almost machine precision level. The validity of the approximate numerical results was confirmed against theory and reported findings from the literature. Despite some inconsistencies with the literature, attributed to the physical characteristics of the investigated study domain, the good agreement between the current study findings with those published previously reveals that the proposed numerical scheme is suitable for solving coupled nonlinear PDEs with environmental fluid dynamics applications. Moreover, in comparison, the proposed MV-SLQLM is more efficient than the previous TV-SCM and MV-SQLM at solving large coupled systems of PDEs.



---

*Source: 7104852-2023-03-15.xml*