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.

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 [13], 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 [918] provided analytical solutions to modelling the transport of contaminants in soil and groundwater systems. Other studies [1921] proposed analytical solutions to solute and water transport equations for verifying numerical solvers. Gao et al. [22] and others [2325] 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 [5255]. As demonstrated on various theoretical diffusion-reaction or advection-diffusion problems [5658], 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 Section 2, 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 (Figure 1). 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.

Naturally, 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 (Figure 1(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] where and are velocity components, and is the pressure of the fluid; and are density and viscosity of the fluid, respectively; and and are, respectively, the inertia resistance coefficient and the permeability coefficient. Following the Darcy-Boussinesq approximation assumption [79] on the vertical flow, is the force due to gravity, and and are, respectively, volumetric coefficients of expansion of heat and contaminants. The terms and in equation (3) represent medium surface temperature and concentration, respectively. The presence of sink terms, in equation (5) and in equation (6), mimics the plant-root uptake activity due to the nonbare ridge surface assumption.

In equation (4), is the apparent heat capacity or the heat flow retardation factor [46], where , , and are, respectively, the volumetric heat capacities of the porous medium, the liquid phase, and the vapor phase; and are the liquid flux and vapor flux densities, respectively; is the apparent thermal conductivity of soil; and is the volumetric latent heat of vaporization of liquid water. Otherwise, in equation (5), , , and are, respectively, unsaturated moisture diffusivity, coefficient of Soret effect of temperature in moisture flux, and unsaturated hydraulic conductivity.

In equation (6), is volumetric water content ; is the hydrodynamic dispersion coefficient ; is the Darcian fluid flux density ; is the soil bulk density ; is the concentration of contaminants adsorbed onto the solid surface ; is the first-order reaction rate coefficient ; and is the contaminant sink term . Additionally, is the distribution coefficient for linear adsorption process , is the volumetric flux of water per unit volume of soil matrix ; is the concentration of contaminants at sink; is the soil matrix retardation factor ; is the thermal gradient diffusion coefficient ; and and are, respectively, seepage velocity components in the horizontal and vertical directions.

The initial and boundary conditions to (2)–(6) are [78]

Applying the following nondimensional variables we obtain (upon dropping the asterisk) the dimensionless equations [78] where is the Prandtl number; is the Reynolds number; is the Eckert number; is the Forchheimer number, where is the Darcy number; is the Grashof number due to soil temperature; is the Grashof number due to soil concentration; and is the Peclet number.

In addition, is the heat capacity ratio parameter, is the heat generation/absorption parameter, is the thermophoresis parameter, is the characteristic time, is the chemical reaction parameter, is the sorption parameter, is the chemical sink parameter, and is the soil matrix retardation factor.

Equation (5) is highly nonlinear and degenerate because of and that also depend on the field variable . For simplicity, we first rewrite equation (5) as where replaces the sink term to describe the “rate of moisture uptake by the plant root in the -direction” only. Then, using the relation [20]

equation (14) is rewritten as

Using the expression [34] where is the “relative water/moisture content,” and the following nondimensional variables lead to a simplified and nondimensional moisture (normalized) equation [78] where is the moisture retardation factor; is the Lewis number or “reciprocal” Luikov number; is the ratio of “apparent” diffusion coefficients; and .

The dimensionless initial and boundary conditions are [78]

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 [80] where is the hydraulic head , is the effective porosity, and is the single direction distance over which the measured change in the hydraulic head is measured . Using the definitions of moisture specific capacity and unsaturated moisture diffusivity yields a new expression

Combining equations (21) and (22) yields such that

Using the relationship between the Darcy flux density and the average seepage velocity [80], together with equation (24), we deduce the following expressions;

Then, using the first two terms of (26), we replace the Darcy flux density components and in equations (12) and (13) with the “average” seepage velocity components and , 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 in equation (12) to be eloquently coupled with equation (19).

Lastly, the parameters , , and in equation (13) are further rearranged. Introducing a “one,” to the RHS of , results into where is the apparent diffusivity of the porous matrix. Hence, equation (27) is expressed in the form of the Damkhler (Da) number [23] when , , and , where , , and are, respectively, half-fracture aperture length, reaction rate constant, and matrix diffusion coefficient. Hence, the nondimensional parameter (equation (27)) is termed “apparent” Damkhler number. In a similar manner, the rearranged forms of and are, respectively, where in (29) is now called the “weighted” Damkhler number.

3. Derivation of Solution Method

3.1. MV-SLQLM Approximation and Construction in () Dimensions

Suppose we have a general nonlinear PDE of the form for , where 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 of each linear equation using the triple Lagrange interpolation functions , , and at the Chebyshev-Gauss-Lobatto (CGL) collocation points located in the local support domain . 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 from the original domain via the linear transformation, where [55] is the collocation points described by the Chebyshev polynomial of order ,

Then, we approximate each function by the Chebyshev interpolation , where together with the corresponding approximate derivatives at the collocation points, where , , and are, respectively, , , and Chebyshev differentiation matrices for .

Further, a Lagrange polynomial of order , for example, at the Chebyshev-Gauss-Lobatto points , is defined as [72] where for , and . In addition, the first order Chebyshev differentiation matrix at the collocation points is given by [52, 55]

For simplicity and efficiency, we express equations (34), (35a), (35b), and (35c) in matrix form, using the Kronecker product, as [72] where , and

Third, we formulate the spectral collocated coupled system, defined for each equation, as

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 for in equation (41a), assuming and are known from the previous iteration(b)With now known, solve for in equation (41b)(c)Solve for in (41c), using and , assuming is known from the previous iteration(d)With now known, proceed to solve for in (41d)(e)Finally, using , , and , solve for in (41e)

Beginning with the initial guesses, steps (a) to (e) are repeated, in turn, for 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 . Let

Given a smooth function . Let , defined in (34), be the best approximation of so that ,

Then, there exists arbitrary points , , and , so that where is the multivariate chebyshev interpolating polynomial of at , , and .

Taking the infinite norm of equation (44) yields

Since is a smooth function on ; then, its , , , and -th partial derivatives are both continuous and upper-bounded, so that where are arbitrary constants in .

Substituting (46) into (45) yields the norm error bound estimate where .

However, (47) can be further simplified. From the linear relationships defined in (31), we have so that from (47),

Using equation (36), so that the -th degree Lagrange polynomial is

Taking the modulus of equation (51) and rearranging the terms yield

Since the Chebyshev interpolating polynomial, equation (33) is the solution to the Chebyshev differential equation [85] for nonnegative integers and , it follows that

Using equation (54) and applying the triangular inequality, equation (52) can then be expressed as

Employing the inequality relation stated in [84] and noting that , then so that (55) becomes

Therefore, from (57), . It then also follows from (49) that

In a similar manner,

Substituting (58), (59), and (60) into (47) yields

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 and making inference of Ilati and Dehghan [82], the result (see equation 3.10 in [86]) means 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 () 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 number; where is the minimum grid size, sufficient for the flow profile to behave physically correctly; is the grid or local Peclet number; and and are the slope and intercept constants, respectively. Since we have a low , , advection-diffusion-reaction equation (13), we investigate the effect of varying on the solute concentration profiles at the minimum grid number, .

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 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 direction where , , and , we employ the QLM to look for the linear expression where

Using equations (66b) and (66c), we have

so that

Hence, equation (65) becomes

Applying the Chebyshev spectral collocation method to equation (69) yields where , , and ; and , , and are, respectively, identity matrices of sizes , , and .

Since 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 where , and

Thus, at step (a) of the LLM, we need to solve for in the matrix equation (equation (71)) defined as where and

In a similar manner, the Darcy-Brinkman-Forchheimer equation (equation (11)) in direction is expressed as where , , , , and .

Applying the QLM, we obtain

Then, by applying the Chebyshev spectral collocation method to equation (77), we get where , , , , , and are all defined as before.

Considering the known values of at the initial and boundary conditions, equation (78) can be rearranged as where , and

Thus, in step (b) of the LLM, we solve for in the matrix equation where and

The pressure term is iteratively estimated from an initial guess deduced from the analytical solution [91] of the pressure equation , assuming an initial state of everywhere and the zero-flux top and bottom boundary conditions.

Next, we express the heat equation as where , , , , and .

In terms of the QLM, (86) can be expressed as

Applying the Chebyshev spectral collocation method to equation (87), we get where , , , , , and are all as previously described.

In a similar manner, with the known values at the boundaries and initial points, we solve where , , , , and

Henceforth, at step (c) of the LLM, we solve the matrix equation, where and,

Since solving the system (91) requires knowledge of , we then consider the moisture equation, which is also coupled with temperature , where . Applying the QLM, we need to solve

Then, applying the Chebyshev spectral collocation method to equation (95), we obtain where , , , , , and are all as previously described.

Taking all the known initial and boundary values of and to the right-hand side of equation (96), we now solve the equation where so that, at step (d) of the LLM, we solve the matrix equation where and

Finally, applying the QLM, the solute equation where and become

Now, applying the Chebyshev spectral collocation to equation (103), we obtain where , , , , , and are all as previously described. Taking out known values of initial and boundary conditions, we obtain where , , , and

Thus, at step (e), we solve for in the matrix equation where and

The spectral-collocated field variables at the initial and boundary points are then given by

Then, starting from initial guesses, we obtain approximate solutions for , , , , and by iteratively solving the matrix equations (73), (81), (91), (99), and (107), in turn, for . 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 Tables 1 and 2. All computations were performed on an Intel core i7 processor of 8th Generation using MATLAB R2015a.

4.1. Stability and Rate of Convergence

At varying local , changes in the solute concentration profile for the given grid numbers are shown in Figure 2.

From Figure 2, stability of the numerical method is more pronounced when solving equation (13) with local , particularly for the collocation points and (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 37 show the infinity error norms at different iteration points for each equation in the model.

Since the present scheme is limited to short temporal domains, the number of collocation points in time was fixed to . 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 , , and (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 and collocation points are sufficient to give accurate and consistent results. The reason is that, unlike at spatial collocation points, jumps and oscillating wiggles appear at the right end of iteration domain in the solute equation for (see Figure 5) and (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 Table 3, 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 in variable and using the same number of collocation points in the spatial and time domains.

From the inspection of Table 3, we observe that for each method, both the computational time and condition number increase with time, 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 , , , , , , , , , , , , , and , otherwise stated. The numerical simulations are discussed in comparison with published results in the literature. In addition, and collocation points in the space and time domain, respectively, sufficient to give accurate and consistent results are used.

4.2.1. Velocity Distribution

Figure 8 illustrates the variation of velocity profile against depth for the vertical preferential flow at changes in , , , , , and , while the other parameters keep constant and read as described above.

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 . 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 , 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 , 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 and .

Unlike in porous channel (single-phase flow) studies, Figure 8(b) shows that the velocity profile decreases with an increase in . The reason is that the velocity distribution decreases to satisfy the conservation of pore momentum as increases because of the condensation process (Figure 1(b)) that occurs at a low number, consistent with [102]. In addition, in agreement with [102], Figure 8(b) shows that the velocity decay rate is slower at high than at low . 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 , , or the decreases, or the Nusselt number increases. Therefore, we also investigated the effect of and in addition to 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(a) shows that the heat distribution profile advances at greater values of , consistent with [102]. When [94] studied the effect of porosity on natural convection of moderate Prandtl number fluids () 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, increases nearly linearly with [94]. The linear relationship between and , as reported in [94], could explain the observed insignificant changes but at high-temperature values.

The Prandtl number is the ratio of momentum to thermal diffusivity. Figure 9(b) shows the effect of on the temperature distribution profile. Since 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 fluids in porous media such as soil, a recent study found that there exist three characteristic heat transfer zones where the Nusselt number (ratio of heat convection to conduction) depends on monotonically [94]. In particular, [94] found that increases with monotonically for moderate Prandtl number fluids but decreases in the zone and slightly increases again for . Therefore, with reference to the zone of values (Figure 9(b)), soil heat increases with due to the enhanced for natural convection flow [92].

Eckert number 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 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 , the temperature profile increases (Figure 9(c)), where means a heat generation and 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 , 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. Figure 10 shows the effect of varying , , , and on changes to moisture distribution. We chose these parameters because, elsewhere, they showed to be sensitive to calculating moisture content.

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 , 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 due to the redistribution process [25, 103] and low values of isothermal moisture diffusivity . Using , Figure 10(d) shows the effect of thermal moisture diffusivity on moisture profile in the investigated domain under combined moisture and temperature gradients. Numerical simulations show that for high , 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, including , , and . In addition, [23] observed the changes in the solute transfer characteristics with the chemical reaction parameter by varying the Damkhler number . Similarly, the effect of reaction and absorption coefficients on solute concentration was investigated, respectively, using the deduced “apparent” Damkhler and “weighted” Damkhler numbers. However, we specifically observed the effects of varying and present in and , respectively. To investigate the effect of plant-root uptake activity, the effect of varying the sink coefficient was observed. The effect of soil temperatures on solute distribution was investigated by varying the thermophoresis parameter . Figure 11 shows the results.

Figure 11(a) shows the effect of varying low Peclet numbers on the solute concentration at low and moderate . As mentioned earlier on, we chose the low range of 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 values. The reason is that high 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 largely depends on (see, equation (29)), the effect of 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 Damkhler numbers 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 . 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 in the effect of “weighted” Damkhler number (see equation (29)) on the distribution profiles of solute concentration is shown in Figure 11(e). Other studies found that an increase in decreases solute concentration because increasing influences the decaying process. In constrast, Figure 11(e) suggests that when both heat and water contents are present in the media, 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 on the distribution of solute concentration in which the distribution is low at higher average values of . The result suggests that the distribution mechanism is now largely facilitated by the increase in dispersion coefficient as the average volumetric water content 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 under saturation. Unlike in the unsaturated flow, 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.

Data Availability

Data sharing not applicable to this article as no datasets were generated or analysed during the current study.

Conflicts of Interest

The authors declare no conflicts of interest.

Acknowledgments

This work was produced as part of PhD study of the corresponding author at Pan African University Institute for Basic Sciences, Technology & Innovation (PAUISTI), funded by the African Union (AU). We are therefore grateful to both AU and PAUISTI.