235 American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) ISSN (Print) 2313-4410, ISSN (Online) 2313-4402 © Global Society of Scientific Research and Researchers http://asrjetsjournal.org/ Comparative Analysis of Laplace Transform and Finite Difference Modeling and Simulation of Solute Transport in Soil. (Case Study: Nitrate Solute Transport in Homogeneous Soil) Ernest Kwame Ampomaha*, Akwasi Ayensub, Barnabas Amisigoc, Obed Agbenorkud aChristian Service University College P O Box 312, Kumasi –Ghana bGraduate School of Nuclear and Allied Science, P O Box AE1, Accra-Ghana cCouncil of Scientific and Industrial Research, P O Box M32, Accra – Ghana dNuclear Regulatory Authority, P O Box AE 50, Accra – Ghana aEmail: freeamps@yahoo.com, bEmail: profayensu@gmail.com cEmail: barnyy2002@yahoo.co.uk, dEmail: abidegh@yahoo.com Abstract Analytical Laplace transform and numerical finite difference methods were used to solve solute transport model (conversion dispersion equation) for a simplified homogeneous soil and simulation of the transport were done using Matlab programming language. Nitrate solute was used for the study. The study compared the simulation results that were generated by both the Laplace transform and the finite difference methods. Spatial and Temporal simulation of nitrate transport comparing both analytical and numerical solutions were presented. The errors in the spatial and temporal numerical solution were simulated. A three dimensional simulation of the nitrate concentration, depth and time for both the Laplace transform and the finite difference method were also presented. The results showed that finite difference numerical method gave a good approximation of the Laplace transform analytical method which provide exact solution. Although there were errors associated with the numerical solution, the output of the numerical solution do not sharply deviate from that of the analytical solution. The errors associated with the finite difference numerical solution were mainly as a result of truncation of the Taylor series expression. ------------------------------------------------------------------------ * Corresponding author. http://asrjetsjournal.org/ American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2017) Volume 38, No 1, pp 235-249 236 Finite difference numerical method can be used to model solute transport in heterogeneous soil which is a more complex process that cannot be accomplished with analytical method with a reasonable level of accuracy. Modeling of solute transport in soil is essential for management of nutrients supply to plants and water resource. Keyword: Laplace Transform; Finite Difference; Model; Simulation; Solute. 1. Introduction Different kinds of solutes can be found in the soil and they include nutrients, salts, pesticides, naturally occurring chemicals and other applied chemicals. A lot of these solutes are of great benefit as they give plants and soil organisms’ food and protection against diseases. However, when these solutes move out of their desired zones in the soil into ground and surface water sources, they can cause a substantial agronomic, environmental, economic consequences [1]. Huge cost are usually associated with solutes that move off site as a result of over application of agrichemicals, ineffective treatment of targeted pests (weeds, insects or diseases) or the remediation of contaminated water sources [2]. Hence, knowledge and understanding of solute transport in soil and how to minimize off-site contamination is required to be able to effectively use chemicals and protect the quality of water resources. Studying transport of solute in soils is necessary to improving crop production and assessing its impact on the environment [3]. Mathematical models are critical for any attempt to effectively understand and estimate site specific subsurface water flow and solute transport processes. Mathematical modelling helps to analyze the existing situation, allows forecasting, and to evaluate the effects of nutrients transport on surrounding water quality. Models are helpful tools for designing, testing and implementing soil, water, and crop management practices that minimize soil and water pollution [4]. There are two methods to mathematically model a system and these are analytical methods and numerical methods. The analytical methods provide exact solution to mathematical problems, however, analytical methods cannot be used to model complex situations. The numerical methods provide approximate solution to mathematical problems but they can be used to model very complex situations. Soil is a heterogeneous medium and hence transport phenomena within it is very complex, as a result only numerical approach can be used to model solute transport within the soil without oversimplifying the soil medium. Since numerical approach provides approximate solution, this paper assessed the margin of error of Finite Difference numerical method as compared with Laplace Transform analytical methods when both approaches are used to model nitrate solute transport in a simplified homogeneous soil. 1.1 Laplace Transform analytical method Laplace transform method uses algebra to solve differential equations. The Laplace transform is defined by: ℒ[𝑓𝑓(𝑡𝑡)](𝑠𝑠) ≡ ∫ 𝑓𝑓(𝑡𝑡)𝑒𝑒−𝑠𝑠𝑠𝑠∞ 0 𝑑𝑑𝑡𝑡……………………………………………………………...1.1 Where 𝑓𝑓(𝑡𝑡) is defined for 𝑡𝑡 ≥ 0 [5] Laplace transforms are used to simplify the governing equations for solute transport. The transform eliminate American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2017) Volume 38, No 1, pp 235-249 237 one independent variable, usually time, and also convert the original transport equation from a partial to an ordinary differential equation. A governing Convection Dispersion Equation (CDE) in the Laplace domain is obtained with the Laplace transform and the equation in this form is much simpler to solve analytically than the original equation. The Laplace Transform of the solute concentration with respect to time is defined as 𝐿𝐿[𝑐𝑐(𝑧𝑧, 𝑡𝑡)] = 𝑐𝑐 �(𝑧𝑧, 𝑠𝑠) = ∫ 𝑐𝑐(𝑧𝑧, 𝑡𝑡)𝑒𝑒−𝑠𝑠𝑠𝑠∞ 0 𝑑𝑑𝑡𝑡…………………………………………………..1.2 Where s is the Laplace transform variable. Analytical solution in the Laplace space must be subsequently inverted to the real space using either a Table of Laplace transforms or by applying inversion theorems and the solution is expressed as an error function or complimentary error function [6]. 1.2 Finite Difference The finite difference method consists of approximation of the differential operator by replacing the derivatives in the equation using differential quotients. The domain of interest is partitioned in space and time, and approximate solutions are computed at the space or time points. Time and space are both divided into small increments t∆ and z∆ or ( x∆ and z∆ ) known as step size as shown in Figure 1.1. Temporal and spatial derivatives are approximated by Taylor series expansion. The accuracy of the approximation is determined by the scheme selected and the mesh sizes of the spatial and temporal domains [6]. Figure 1.1: Spatial and temporal finite difference discretization [6] 1.3 Solute Transport Solute transport in soil results from convection of the dissolved substances, molecular or ionic diffusion and mechanical dispersion but the transport process is generally assumed to be convection-dominated process. Other factors including soil matrix-solute interaction and decay phenomena may affect transport of solute in soils. 1.3.1 Convection (Mass) Flow Convection is passive movement of dissolved constituents of solute with water flowing through the soil, American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2017) Volume 38, No 1, pp 235-249 238 whereby solute and the water move at the same average rate. It is also known as Darcian flow [7]. Convection describes the bulk movement of solute particles along the mean direction of fluid flow at a rate equal to the average interstitial fluid velocity. Convection does not consider microscopic processes but follows the bulk Darcian flow vectors, and is therefore described as the transport along path lines [8]. The velocity at which solutes move through soil matrix is known as pore water velocity, and is described by the ratio of Darcian velocity and moisture content. In general pore water velocity accounts for straight line of length of path traversed in the soil in a given time [9]. The solute flux density, 𝐽𝐽𝑐𝑐 for convective transport is defined; 𝐽𝐽𝑐𝑐 = 𝑞𝑞𝐶𝐶…………………………………………………………………................................... 1.3 𝐶𝐶 = dissolved solute concentration, The water flux density, q is expressed as 𝑞𝑞 = 𝑣𝑣𝑣𝑣………………………………………………………………………………………... 1.4 𝑣𝑣 = pore water velocity, 1.3.2 Diffusion Diffusion is a spontaneous process where solute ions and molecules move from locations of higher to lower concentration as a result of thermal random motion of dissolved ions and molecules. Diffusion is an active process and tends to decrease existing concentration gradients and move the process towards homogeneity rather rapidly [9]. Rates of molecular diffusion are independent of soil water velocity, and diffusion occurs even in the absence of fluid movement [8]. Diffusion flux spreads solute through a concentration gradient. Diffusion is a dominant transport mechanism when convection is insignificant, and is usually a negligible transport mechanism when convection process is very high. Fick’s law defines the diffusive transport as: 𝐽𝐽𝑑𝑑𝑖𝑖𝑓𝑓𝑓𝑓 = −𝑣𝑣𝐷𝐷𝑑𝑑𝑖𝑖𝑓𝑓𝑓𝑓 ∂C ∂z …………………………………………………………………………….. 1.5 𝐽𝐽𝑑𝑑𝑖𝑖𝑓𝑓𝑓𝑓 = solute flux density for diffusion, z = soil depth 𝐷𝐷𝑑𝑑𝑖𝑖𝑓𝑓𝑓𝑓 = diffusion coefficient in soils. 𝐷𝐷𝑑𝑑𝑖𝑖𝑓𝑓𝑓𝑓 = 𝐷𝐷𝑜𝑜θε……………………………………………………….…………………………... 1.6 𝐷𝐷𝑜𝑜 = diffusion coefficient in pure water, ε = tortuosity 1.3.3 Dispersion Dispersion is the mixing and spreading of solutes along and transverse to the direction of flow in response to American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2017) Volume 38, No 1, pp 235-249 239 local variations in interstitial fluid velocities. Dispersive transport of solute occurs due to the uneven distribution of water flow velocities within and between different soil pores. Dispersion is a passive process. Macroscopically, dispersion is similar to diffusion; however it occurs only during water movement and not driven by concentration gradients [9]. The dispersion that occurs along the direction of flow path is called longitudinal dispersion and that in the direction normal to flow is known as transverse dispersion. Dispersion process and diffusion process are considered to be additive at macroscopic level. The dispersive transport is described by an equation similar to diffusion as: 𝐽𝐽𝑑𝑑𝑖𝑖𝑠𝑠 = −𝑣𝑣𝐷𝐷𝑑𝑑𝑖𝑖𝑠𝑠 ∂C ∂z ………………………………………..….………………………………… 1.7 𝐽𝐽𝑑𝑑𝑖𝑖𝑠𝑠 = solute flux for dispersion, 𝐷𝐷𝑑𝑑𝑖𝑖𝑠𝑠 = dispersion coefficient. Dispersion is assumed to be a function of fluid velocity as: 𝐷𝐷𝑑𝑑𝑖𝑖𝑠𝑠 = λvn……………………………………………………………………………………... 1.8 λ = dispersivity, n = empirical constant (generally assumed to be 1). Dispersion process and diffusion process are considered to be additive at macroscopic level. The two are therefore combined to define a new parameter called the apparent dispersion coefficient or hydrodynamic dispersion coefficient (D) D = 𝐷𝐷𝑑𝑑𝑖𝑖𝑓𝑓𝑓𝑓 + 𝐷𝐷𝑑𝑑𝑖𝑖𝑠𝑠 ……………………………………………………..……………………….. 1.9 1.4 Nitrate Solute Nitrate (NO3 -) is negatively charged ion. It is very mobile in soils and can easily be lost from the soil with water that moves downward laterally through a soil profile. The surface of the negatively charged clay or organic matter particles repels nitrate rather than been attracted and therefore it can be lost by leaching. Movement of the NO3 - ion through soil is governed by convection, or mass-flow, with the moving soil solution and by diffusion and dispersion within the soil solution. The widespread appearance of NO3 - in ground water is a consequence of its high solubility, mobility, and easy displacement by water [10]. Nitrate is a potential pollutant if it reaches surface and ground water supplies. 2. Methodology The governing equation for the modelling process was the one dimensional convectional dispersion equation (CDE). The CDE was solve using both Laplace transform and finite difference methods and the solutions were implemented in mat- lab programming environment. 2.1 Solute Transport Equation American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2017) Volume 38, No 1, pp 235-249 240 The convection-dispersion equation (CDE) which is the accepted deterministic solute transport equation describes the time rate of change of solute concentration for a single solute [8]. The CDE is a partial differential equation of parabolic type, derived on the principle of conservation of mass using Fick’s law. The analytical/numerical solutions of the CDE along with initial and boundary conditions help to understand the solute concentration profile or distribution behavior through an open medium like air, rivers, lakes and porous medium on the basis for which remedial measures to reduce or eliminate the damages may be implemented [11]. 2.1.1 Convection-Dispersion Equation The classic one dimensional CDE for transport of conservative species without adsorption or decay in a partially saturated porous medium can be written as [7]: 𝜕𝜕(𝜃𝜃𝜃𝜃) 𝜕𝜕𝑠𝑠 = 𝜕𝜕 𝜕𝜕𝜕𝜕 �𝑣𝑣𝐷𝐷 𝜕𝜕𝜃𝜃 𝜕𝜕𝜕𝜕 − 𝑞𝑞𝐶𝐶�…………………………………………………………..… (2.1) The comprehensive CDE for one-dimensional transport of reactive solutes, subject to reaction terms of adsorption, first-order degradation, and zero-order production is given as [12]; 𝜕𝜕 𝜕𝜕𝑠𝑠 (𝑣𝑣𝐶𝐶 + 𝜌𝜌𝑏𝑏𝑆𝑆𝑠𝑠 ) = 𝜕𝜕 𝜕𝜕𝜕𝜕 �𝑣𝑣𝐷𝐷 𝜕𝜕𝜃𝜃 𝜕𝜕𝜕𝜕 − 𝑞𝑞𝐶𝐶� − 𝑣𝑣𝜇𝜇𝑙𝑙𝐶𝐶 − 𝜌𝜌𝑏𝑏𝜇𝜇𝑠𝑠𝑆𝑆𝑠𝑠 + 𝑣𝑣𝛾𝛾𝑙𝑙(𝑧𝑧) + 𝜌𝜌𝑏𝑏𝛾𝛾𝑆𝑆(𝑧𝑧……… (2.2) 𝑆𝑆𝑠𝑠 = 𝑘𝑘𝑑𝑑C......................................................................................................................... (2.3) Assuming reversible equilibrium adsorption and steady state flow in a homogeneous soil, equation (2.2) reduces to 𝑅𝑅 𝜕𝜕𝜃𝜃 𝜕𝜕𝑠𝑠 = 𝐷𝐷 𝜕𝜕2𝜃𝜃 𝜕𝜕𝜕𝜕2 − 𝑞𝑞 𝜃𝜃 𝜕𝜕𝜃𝜃 𝜕𝜕𝜕𝜕 − 𝜇𝜇𝐶𝐶 + 𝛾𝛾(𝑧𝑧) ……..………………............................................... (2.4) R=1+ 𝜌𝜌𝑏𝑏𝐾𝐾𝑑𝑑 𝜃𝜃 …………………..………………………..................................................... (2.5) 𝜇𝜇 = 𝜇𝜇𝑙𝑙 + 𝜌𝜌𝑏𝑏𝐾𝐾𝑑𝑑𝜇𝜇𝑠𝑠 𝜃𝜃 …………...………………………....................................................... (2.6) 𝛾𝛾(𝑧𝑧) = 𝛾𝛾𝑙𝑙(z) + 𝜌𝜌𝑏𝑏𝐾𝐾𝑑𝑑𝛾𝛾𝑠𝑠 𝜃𝜃 …..…………………………………......................................... (2.7) 𝜇𝜇𝑙𝑙 , 𝜇𝜇𝑠𝑠 = First-order decay coefficients for degradation of the solute in the liquid and adsorbed phases respectively. 𝛾𝛾𝑙𝑙 𝑎𝑎𝑎𝑎𝑑𝑑 𝛾𝛾𝑆𝑆 = zero order production terms for the liquid and adsorbed respectively. R = retardation factor, 𝜇𝜇 and 𝛾𝛾 = combined first and zero order rate coefficients respectively. 𝑅𝑅 𝜕𝜕𝜃𝜃 𝜕𝜕𝑠𝑠 = 𝐷𝐷 𝜕𝜕2𝜃𝜃 𝜕𝜕𝜕𝜕2 − 𝑣𝑣 𝜕𝜕𝜃𝜃 𝜕𝜕𝜕𝜕 − 𝜇𝜇𝐶𝐶 + 𝛾𝛾(𝑧𝑧) ………………………….……………………................ (2.8) American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2017) Volume 38, No 1, pp 235-249 241 2.2 Assumptions • Nitrate transport occurs in the vertical direction only. • The soil is homogeneous and unsaturated. • Nitrate transport is affected by only convection, and hydrodynamic dispersion (combined dispersion and diffusion) processes. Any other solute process is negligible. • Nitrate is non-adsorbing solute 2.3 Initial and Boundary Conditions 𝐶𝐶(𝑧𝑧, 0) = 0 0 ≤ 𝑧𝑧 ≤ ∞, 𝑡𝑡 = 0 ……..………………………………………………. (2.9) 𝐶𝐶(0, 𝑡𝑡) = 𝐶𝐶0 𝑧𝑧 = 0, 𝑡𝑡 > 0…..………...……………………………………………… (2.10) 𝐶𝐶(∞, 𝑡𝑡) = 0 𝑧𝑧 → ∞, 𝑡𝑡 > 0 …...………………………….………………………….. (2.11) 2.4 Laplace Transform Solution of the CDE Applying the assumptions, equation (2.8) reduces to ∂c ∂t = D ∂2c ∂z2 − v ∂c ∂z ……………………………………………….………….……… ……… (2.12) Laplace Transform with respect to time is defined as 𝐿𝐿[𝑐𝑐(𝑧𝑧, 𝑡𝑡)] = 𝑐𝑐 �(𝑧𝑧, 𝑠𝑠) = ∫ 𝑐𝑐(𝑧𝑧, 𝑡𝑡)𝑒𝑒−𝑠𝑠𝑠𝑠∞ 0 𝑑𝑑𝑡𝑡 ……..………………………………………… (2.13) Applying Laplace transform with respect to time to equation (2.12) gave: 𝐿𝐿 �𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 � = 𝐿𝐿 �𝐷𝐷 𝜕𝜕2𝑦𝑦 𝜕𝜕𝜕𝜕2 � − 𝐿𝐿 �𝑣𝑣 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 �…………………………………….…………………………................ (2.14) The Laplace transforms are; [13] 𝐿𝐿{𝑓𝑓(𝑡𝑡)} = ∫ 𝑓𝑓(𝑡𝑡)𝑒𝑒−𝑠𝑠𝑠𝑠𝑑𝑑𝑡𝑡 = 𝐹𝐹(𝑠𝑠)∞ 0 …….………………………………………...……………………... (2.15) 𝐿𝐿{𝑓𝑓′(𝑡𝑡)} = 𝑠𝑠𝐹𝐹(𝑠𝑠) − 𝑓𝑓(0)…………...……………………………………………………………….… (2.16) 𝐿𝐿 �𝜕𝜕𝜕𝜕 𝜕𝜕𝑠𝑠 � = 𝑠𝑠𝑐𝑐̅ − 𝑐𝑐(𝑧𝑧, 0)………….…………………………………………………..........................................(2.17) Applying the initial condition, equation (2.17) becomes: 𝐿𝐿 �𝜕𝜕𝜕𝜕 𝜕𝜕𝑠𝑠 � = 𝑠𝑠𝑐𝑐̅ − 𝑐𝑐(𝑧𝑧, 0)……………………………………………………………………………………... (2.18) American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2017) Volume 38, No 1, pp 235-249 242 𝐿𝐿 �𝑣𝑣 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 � = 𝑣𝑣 ∫ 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 𝑒𝑒−𝑠𝑠𝑠𝑠𝑑𝑑𝑡𝑡∞ 0 …….…………………………………………………………………………. (2.19) Since the integral of derivative = derivative of integral, then by exchanging the order of integration and differentiation in eqn. (2.19) becomes: 𝐿𝐿 �𝑣𝑣 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 � = 𝑣𝑣 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 ∫ 𝑐𝑐(𝑧𝑧, 𝑡𝑡)𝑒𝑒−𝑠𝑠𝑠𝑠∞ 0 𝑑𝑑𝑡𝑡……………………………………………………………………….. (2.20) 𝐿𝐿 �𝑣𝑣 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 � = 𝑣𝑣 𝜕𝜕𝜕𝜕̅ 𝜕𝜕𝜕𝜕 ………………………………………………………………………………………………. (2.21) 𝐿𝐿 �𝐷𝐷 𝜕𝜕2𝜕𝜕 𝜕𝜕𝜕𝜕2 � = 𝐷𝐷 ∫ 𝜕𝜕2𝜕𝜕 𝜕𝜕𝜕𝜕2 𝑒𝑒−𝑠𝑠𝑠𝑠𝑑𝑑𝑡𝑡∞ 0 …….………………………………………………………………………. (2.22) 𝐿𝐿 �𝐷𝐷 𝜕𝜕2𝜕𝜕 𝜕𝜕𝜕𝜕2 � = 𝐷𝐷 𝜕𝜕2 𝜕𝜕𝜕𝜕2 ∫ 𝑐𝑐(𝑧𝑧, 𝑡𝑡) 𝑒𝑒−𝑠𝑠𝑠𝑠𝑑𝑑𝑡𝑡∞ 0 …….……………………………………………………………… (2.23) 𝐿𝐿 �𝐷𝐷 𝜕𝜕2𝜕𝜕 𝜕𝜕𝜕𝜕2 � = 𝐷𝐷 𝜕𝜕2𝜕𝜕̅ 𝜕𝜕𝜕𝜕2 …………….……………………………………………………………………….…….. (2.24) Putting eqn. (2.18), (2.21) and (2.24) into eqn. (2.14) yielded: 𝑠𝑠𝑐𝑐̅ = 𝐷𝐷 𝜕𝜕2𝜕𝜕̅ 𝜕𝜕𝜕𝜕2 − 𝑣𝑣 𝜕𝜕𝜕𝜕̅ 𝜕𝜕𝜕𝜕 …………………………………………………………………………..……………….. (2.25) Taking Laplace transform of the boundary conditions 𝐿𝐿[𝑐𝑐(0, 𝑡𝑡)] = 𝐿𝐿[𝑐𝑐0]……….………………………………………………………………………………...... (2.26) 𝑐𝑐̅ (0, 𝑠𝑠) = 𝑐𝑐0 ∫ 𝑒𝑒−𝑠𝑠𝑠𝑠𝑑𝑑𝑡𝑡 = 𝑐𝑐0 ∞ 0 �𝑒𝑒 −𝑠𝑠𝑠𝑠 𝑠𝑠 � 0 ∞ = 𝜕𝜕0 𝑠𝑠 ………….……………………………………………….... (2.27) 𝐿𝐿[(∞, 𝑡𝑡)] = 𝐿𝐿[0]……...……………………………………………………………………………............. (2.28) ∫ 𝑐𝑐(∞, 𝑡𝑡)𝑒𝑒−𝑠𝑠𝑠𝑠𝑑𝑑𝑡𝑡 = 0∞ 0 …………………………………………………………….…………………….... (2.29) 𝑐𝑐̅ (∞, 𝑠𝑠) = 0... ……………………………………………………..……………..……………………….… (2.30) Dividing through equation (2.25) by D and rearranging gives: 𝜕𝜕2𝜕𝜕̅ 𝜕𝜕𝜕𝜕2 − 𝑣𝑣 𝐷𝐷 𝜕𝜕𝜕𝜕̅ 𝜕𝜕𝜕𝜕 − 1 𝐷𝐷 𝑠𝑠𝑐𝑐̅ = 0……..………………………………………………………………………………… (2.31) Which is a second order ordinary differential equation with auxiliary equation 𝑚𝑚2 − 𝑣𝑣 𝐷𝐷 𝑚𝑚 − 𝑠𝑠 𝐷𝐷 = 0…………………………………………………….………………………………........ (2.32) American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2017) Volume 38, No 1, pp 235-249 243 The roots of equation (2.32) were determined from the quadratic formula. 𝑚𝑚 = −𝑏𝑏±�𝑏𝑏2−4𝑎𝑎𝜕𝜕 2𝑎𝑎 ………...……………………………………………….………..……………………….... (2.33) Hence equation (2.32) had the roots: 𝑚𝑚 = 1 2 �𝑣𝑣 𝐷𝐷 ± �𝑣𝑣2 𝐷𝐷2 + 4𝑠𝑠 𝐷𝐷 �…… …. ………………………………………………..…………………………. (2.34) The general solution of equation (2.31) is given as [13] 𝑐𝑐̅ = 𝐴𝐴𝑒𝑒𝐴𝐴𝐴𝐴 �𝑣𝑣𝜕𝜕 2𝐷𝐷 + 𝜕𝜕 2 �𝑣𝑣2 𝐷𝐷2 + 4𝑠𝑠 𝐷𝐷 � + 𝐵𝐵𝑒𝑒𝐴𝐴𝐴𝐴 �𝑣𝑣𝜕𝜕 2𝐷𝐷 − 𝜕𝜕 2 �𝑣𝑣2 𝐷𝐷2 + 4𝑠𝑠 𝐷𝐷 �……………...…………….………….. (2.35) Applying the lower boundary condition [𝑐𝑐̅ (∞, 𝑠𝑠) = 0]means as 𝑧𝑧 → ∞ there would be no concentration, hence A=0, and equation (2.35) reduced to: 𝑐𝑐̅ = 𝐵𝐵𝑒𝑒𝐴𝐴𝐴𝐴 �𝑣𝑣𝜕𝜕 2𝐷𝐷 − 𝜕𝜕 2 �𝑣𝑣2 𝐷𝐷2 + 4𝑠𝑠 𝐷𝐷 �……….……………………………………………………………….. (2.36) Applying the upper (surface) boundary condition, �𝑐𝑐̅(0, 𝑠𝑠) = 𝜃𝜃0 𝑠𝑠 � to equation (2.36) yields, 𝐵𝐵 = 𝜃𝜃0 𝑠𝑠 and eqn. (2.36) becomes: 𝑐𝑐̅ = 𝜃𝜃0 𝑠𝑠 𝑒𝑒𝐴𝐴𝐴𝐴 �𝑣𝑣𝜕𝜕 2𝐷𝐷 − 𝜕𝜕 2 �𝑣𝑣2 𝐷𝐷2 + 4𝑠𝑠 𝐷𝐷 �…………..……………………………………………………………. (2.37) Rearranging equation (2.37) produced: 𝜕𝜕̅ 𝜃𝜃0 = 1 𝑠𝑠 𝑒𝑒𝐴𝐴𝐴𝐴 �𝑣𝑣𝜕𝜕 2𝐷𝐷 � 𝑒𝑒𝐴𝐴𝐴𝐴�−𝑧𝑧� 𝑠𝑠+𝑣𝑣 2 4𝐷𝐷 𝐷𝐷 �……………....……………………………….…………………..… (2.38) Letting 𝑄𝑄 = 𝑠𝑠 + 𝑣𝑣2 4𝐷𝐷 ⇉ 𝑠𝑠 = 𝑄𝑄 − 𝑣𝑣2 4𝐷𝐷 𝜕𝜕̅ 𝜃𝜃0 = � 1 𝑄𝑄−𝑣𝑣 2 4𝐷𝐷 � 𝑒𝑒𝐴𝐴𝐴𝐴 �𝑣𝑣𝜕𝜕 2𝐷𝐷 � 𝑒𝑒𝐴𝐴𝐴𝐴 �−𝑧𝑧 �𝑄𝑄 𝐷𝐷 � 0.5 �……………… ……………………………………..……….. (2.39) From Shift theorem [14] 𝐿𝐿−1[𝐹𝐹(𝑠𝑠 + 𝑎𝑎)] = 𝑒𝑒−𝑎𝑎𝑠𝑠𝐿𝐿−1[𝐹𝐹(𝑠𝑠)]………………………………………………. ………… (2.40) American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2017) Volume 38, No 1, pp 235-249 244 𝐿𝐿−1[𝐹𝐹(𝑄𝑄)] = 𝐿𝐿−1 �𝐹𝐹 �𝑠𝑠 + 𝑣𝑣2 4𝐷𝐷 �� = 𝑒𝑒𝐴𝐴𝐴𝐴 �− 𝑣𝑣2𝑠𝑠 4𝐷𝐷 � 𝐿𝐿−1[𝐹𝐹(𝑠𝑠)]………...…….………………….............. (2.41) The Laplace inverse of equation (2.39) was written as: 𝜃𝜃(𝜕𝜕,𝑠𝑠) 𝜃𝜃0 = exp �𝑣𝑣𝜕𝜕 2𝐷𝐷 � 𝑒𝑒𝐴𝐴𝐴𝐴 �− 𝑣𝑣2𝑠𝑠 4𝐷𝐷 � 𝐿𝐿−1 � 1 𝑄𝑄−𝑣𝑣 2 4𝐷𝐷 𝑒𝑒𝐴𝐴𝐴𝐴 �−𝑧𝑧 �𝑄𝑄 𝐷𝐷 � 0.5 ��…………………………..………………… (2.42) Carslaw and Jaeger [15] have shown that the inverse Laplace expression on RHS of eqn. (2.42) could be written as: 𝐿𝐿−1 � 1 𝑄𝑄−𝑣𝑣 2 4𝐷𝐷 𝑒𝑒𝐴𝐴𝐴𝐴 �−𝑧𝑧 �𝑄𝑄 𝐷𝐷 � 0.5 �� = 1 2 𝑒𝑒𝐴𝐴𝐴𝐴 �𝑣𝑣 2𝑠𝑠 4𝐷𝐷 � ⎣ ⎢ ⎢ ⎢ ⎡𝑒𝑒𝐴𝐴𝐴𝐴 �− 𝑣𝑣𝜕𝜕 2𝐷𝐷 � 𝑒𝑒𝑒𝑒𝑓𝑓𝑐𝑐 � 𝜕𝜕 2√𝐷𝐷𝑠𝑠 − �𝑣𝑣2𝑠𝑠 4𝐷𝐷 � + 𝑒𝑒𝐴𝐴𝐴𝐴 �𝑣𝑣𝜕𝜕 2𝐷𝐷 � 𝑒𝑒𝑒𝑒𝑓𝑓𝑐𝑐 � 𝜕𝜕 2√𝐷𝐷𝑠𝑠 + �𝑣𝑣2𝑠𝑠 4𝐷𝐷 � ⎦ ⎥ ⎥ ⎥ ⎤ ………… (2.43) Substituting equation (2.43) into equation (2.42) gave 𝜕𝜕(𝜕𝜕,𝑠𝑠) 𝜕𝜕0 = 1 2 exp �𝑣𝑣𝜕𝜕 2𝐷𝐷 � 𝑒𝑒𝐴𝐴𝐴𝐴 �− 𝑣𝑣2𝑠𝑠 4𝐷𝐷 � 𝑒𝑒𝐴𝐴𝐴𝐴 �𝑣𝑣 2𝑠𝑠 4𝐷𝐷 � ⎣ ⎢ ⎢ ⎢ ⎡𝑒𝑒𝐴𝐴𝐴𝐴 �− 𝑣𝑣𝜕𝜕 2𝐷𝐷 � 𝑒𝑒𝑒𝑒𝑓𝑓𝑐𝑐 � 𝜕𝜕 2√𝐷𝐷𝑠𝑠 − �𝑣𝑣2𝑠𝑠 4𝐷𝐷 �+ 𝑒𝑒𝐴𝐴𝐴𝐴 �𝑣𝑣𝜕𝜕 2𝐷𝐷 � 𝑒𝑒𝑒𝑒𝑓𝑓𝑐𝑐 � 𝜕𝜕 2√𝐷𝐷𝑠𝑠 + �𝑣𝑣2𝑠𝑠 4𝐷𝐷 � ⎦ ⎥ ⎥ ⎥ ⎤ …………. (2.44) 𝜕𝜕(𝜕𝜕,𝑠𝑠) 𝜕𝜕0 = 1 2 �𝑒𝑒𝐴𝐴𝐴𝐴 �− 𝑣𝑣𝜕𝜕 2𝐷𝐷 � 𝑒𝑒𝑒𝑒𝑓𝑓𝑐𝑐 � 𝜕𝜕 2√𝐷𝐷𝑠𝑠 − �𝑣𝑣2𝑠𝑠 4𝐷𝐷 � + 𝑒𝑒𝐴𝐴𝐴𝐴 �𝑣𝑣𝜕𝜕 2𝐷𝐷 � 𝑒𝑒𝑒𝑒𝑓𝑓𝑐𝑐 � 𝜕𝜕 2√𝐷𝐷𝑠𝑠 + �𝑣𝑣2𝑠𝑠 4𝐷𝐷 � �…...…….. … (2.45) 𝜕𝜕(𝜕𝜕,𝑠𝑠) 𝜕𝜕0 = 𝜕𝜕0 2 �𝑒𝑒𝑒𝑒𝑓𝑓𝑐𝑐 � 𝜕𝜕−𝑣𝑣𝑠𝑠 2√𝐷𝐷𝑠𝑠 �+ 𝑒𝑒𝐴𝐴𝐴𝐴 �𝑣𝑣𝜕𝜕 𝐷𝐷 � 𝑒𝑒𝑒𝑒𝑓𝑓𝑐𝑐 �𝜕𝜕+𝑣𝑣𝑠𝑠 2√𝐷𝐷𝑠𝑠 � �………….……………………………..(2.46 2.5 Finite Difference Solution of the CDE The CDE was also solved numerically by fully implicit (backward in time) finite difference scheme. The model domain was discretized into grid points using space and time steps of 1cm and 1 hour respectively. Applying the assumptions, equation (2.8) reduces to ∂c ∂t = D ∂2c ∂z2 − v ∂c ∂z ……………………………………………….………….……. ……….. (2.47) Converting the CDE equation (3.24) into a difference equation gives: 𝜕𝜕𝜕𝜕 𝜕𝜕𝑠𝑠 = 𝜃𝜃𝑖𝑖 𝑗𝑗+1−𝜃𝜃𝑖𝑖 𝑗𝑗 ∆𝑠𝑠 ……………………………………………………........................................ (2.48) American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2017) Volume 38, No 1, pp 235-249 245 𝐷𝐷 ∂2c 𝜕𝜕𝜕𝜕2 = 𝐷𝐷 𝜃𝜃𝑖𝑖+1 𝑗𝑗+1−2𝜃𝜃𝑖𝑖 𝑗𝑗+1+𝜃𝜃𝑖𝑖−1 𝑗𝑗+1 (∆𝜕𝜕)2 ……………………………………………………………….. (2.49) 𝑣𝑣 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 = 𝑣𝑣 𝜃𝜃𝑖𝑖+1 𝑗𝑗+1−𝜃𝜃𝑖𝑖−1 𝑗𝑗+1 2∆𝜕𝜕 ……………………………………………………………………….. (2.50) Substituting equation (3.63), (3.64), and (3.65) into equation (3.24) 𝜃𝜃𝑖𝑖 𝑗𝑗+1−𝜃𝜃𝑖𝑖 𝑗𝑗 ∆𝑠𝑠 = 𝐷𝐷 𝜃𝜃𝑖𝑖+1 𝑗𝑗+1−2𝜃𝜃𝑖𝑖 𝑗𝑗+1−𝜃𝜃𝑖𝑖−1 𝑗𝑗+1 (∆𝜕𝜕)2 − 𝑣𝑣 𝜃𝜃𝑖𝑖+1 𝑗𝑗+1−𝜃𝜃𝑖𝑖−1 𝑗𝑗+1 2∆𝜕𝜕 ……………………………………………… (2.51) 𝐶𝐶𝑖𝑖 𝑗𝑗+1 − 𝐶𝐶𝑖𝑖 𝑗𝑗 = 𝐷𝐷(∆𝑠𝑠) (∆𝜕𝜕)2 �𝐶𝐶𝑖𝑖+1 𝑗𝑗+1 − 2𝐶𝐶𝑖𝑖 𝑗𝑗+1 + 𝐶𝐶𝑖𝑖−1 𝑗𝑗+1� − 𝑣𝑣∆𝑠𝑠 2∆𝜕𝜕 �𝐶𝐶𝑖𝑖+1 𝑗𝑗+1 − 𝐶𝐶𝑖𝑖−1 𝑗𝑗+1�……………………... (2.52) Let 𝐴𝐴 = 𝐷𝐷∆𝑠𝑠 (∆𝜕𝜕)2 , 𝑞𝑞 = 𝑣𝑣∆𝑠𝑠 2∆𝜕𝜕 𝐶𝐶𝑖𝑖 𝑗𝑗+1 − 𝐶𝐶𝑖𝑖 𝑗𝑗 = 𝐴𝐴�𝐶𝐶𝑖𝑖+1 𝑗𝑗+1 − 2𝐶𝐶𝑖𝑖 𝑗𝑗+1 + 𝐶𝐶𝑖𝑖−1 𝑗𝑗+1� − 𝑞𝑞�𝐶𝐶𝑖𝑖+1 𝑗𝑗+1 − 𝐶𝐶𝑖𝑖−1 𝑗𝑗+1�……………………………. (2.53) 𝐶𝐶𝑖𝑖 𝑗𝑗+1 − 𝐶𝐶𝑖𝑖 𝑗𝑗 = 𝐴𝐴𝐶𝐶𝑖𝑖+1 𝑗𝑗+1 − 2𝐴𝐴𝐶𝐶𝑖𝑖 𝑗𝑗+1 + 𝐴𝐴𝐶𝐶𝑖𝑖−1 𝑗𝑗+1 − 𝑞𝑞𝐶𝐶𝑖𝑖+1 𝑗𝑗+1 + 𝑞𝑞𝐶𝐶𝑖𝑖−1 𝑗𝑗+1……………………………. (2.54) −(𝐴𝐴 + 𝑞𝑞)𝐶𝐶𝑖𝑖−1 𝑗𝑗+1 + (1 + 2𝐴𝐴)𝐶𝐶𝑖𝑖 𝑗𝑗+1 + (𝑞𝑞 − 𝐴𝐴)𝐶𝐶𝑖𝑖+1 𝑗𝑗+1 = 𝐶𝐶𝑖𝑖 𝑗𝑗 ………………………………….. (2.55) When i =1 equation (3.70) reduces to: −(𝐴𝐴 + 𝑞𝑞)𝐶𝐶0 𝐽𝐽+1 + (1 + 2𝐴𝐴)𝐶𝐶1 𝑗𝑗+1 + (𝑞𝑞 − 𝐴𝐴)𝐶𝐶2 𝑗𝑗+1 = 𝐶𝐶1 𝑗𝑗………………………………… (2.56) Applying the boundary condition 𝐶𝐶(0, 𝑡𝑡) = 𝐶𝐶0 , 𝑧𝑧 = 0, 𝑡𝑡 > 0 equation ( ) becomes (1 + 2𝐴𝐴)𝐶𝐶1 𝑗𝑗+1 + (𝑞𝑞 − 𝐴𝐴)𝐶𝐶2 𝑗𝑗+1 = 𝐶𝐶1 𝑗𝑗 + (𝐴𝐴 + 𝑞𝑞)𝐶𝐶0 𝑗𝑗+1………………………………….. (2.57) When 𝑖𝑖 = 2 → 𝑎𝑎 − 1, equation (3.71) remains: −(𝐴𝐴 + 𝑞𝑞)𝐶𝐶𝑖𝑖−1 𝑗𝑗+1 + (1 + 2𝐴𝐴)𝐶𝐶𝑖𝑖 𝑗𝑗+1 + (𝑞𝑞 − 𝐴𝐴)𝐶𝐶𝑖𝑖+1 𝑗𝑗+1 = 𝐶𝐶𝑖𝑖 𝑗𝑗 ……………………………….. (2.58) n = number of grid points When i = n −(𝐴𝐴 + 𝑞𝑞)𝐶𝐶𝑛𝑛−1 𝑗𝑗+1 + (1 + 2𝐴𝐴)𝐶𝐶𝑛𝑛 𝑗𝑗+1 + (𝑞𝑞 − 𝐴𝐴)𝐶𝐶𝑛𝑛+1 𝑗𝑗+1 = 𝐶𝐶𝑛𝑛 𝑗𝑗 ……………………………...... (2.59) Applying the boundary condition 𝐶𝐶(∞, 𝑡𝑡) = 0, 𝑧𝑧 → ∞, 𝑡𝑡 > 0 equation ( ) becomes −(𝐴𝐴 + 𝑞𝑞)𝐶𝐶𝑛𝑛−1 𝑗𝑗+1 + (1 + 2𝐴𝐴)𝐶𝐶𝑛𝑛 𝑗𝑗+1 = 𝐶𝐶𝑛𝑛 𝑗𝑗 ……………………………………………….. (2.60) American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2017) Volume 38, No 1, pp 235-249 246 Assuming 𝐴𝐴1 = 𝐴𝐴 + 𝑞𝑞, 𝐴𝐴2 = 1 + 2𝐴𝐴 𝐴𝐴3 = 𝑞𝑞 − 𝐴𝐴, the above sets of simultaneous equations can be expressed in a tri-diagonal matrix notation as follows: ( )1 1 2 3 1 1 0 1 1 2 3 2 2 1 1 2 3 3 3 1 1 2 3 2 2 1 1 2 3 1 1 1 1 2 0 0 0 0 0 0 0 0 0 0 0 0 j j j j j j j j j n n j j n n j j n n x x C C q p C x x x C C x x x C C x x x C C x x x C C x x C C + + + + + − − + − − +    + +      −          −      ↓ ↓     =     ↓ ↓      −          −      −           � � ……….. (2.61) 3. Result and Discussion Figure 3.1: Analytical and Numerical simulation of nitrate transport with respect to soil depth Figure 3.2: Analytical and Numerical simulation of nitrate transport with respect to time 0 10 20 30 40 50 60 70 80 90 100 0 1 2 3 4 5 6 7 8 Depth(z),cm Co nce ntra tion ,g/c m3 Analytical Numerical 0 5 10 15 20 25 30 35 40 45 50 0 1 2 3 4 5 6 7 8 time,hrs con cen trat ion, g/c m3 Analytical Numerical American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2017) Volume 38, No 1, pp 235-249 247 Figure 3.3: analytical simulation of concentration, time and depth Figure 3.4: 3D numerical simulation of concentration, time and depth Figure 3.5: Plot of error in temporal numerical model against time Figure 3.5: Plot of error in spatial numerical model against depth 0 20 40 60 80 100 120 0 10 20 30 40 50 0 2 4 6 8 Time (hrs) Depth (cm) Co nc en tra tio n (g cm 3 ) 0 1 2 3 4 5 6 7 0 20 40 60 80 100 120 0 10 20 30 40 50 0 1 2 3 4 5 6 7 8 Time (hrs) Depth (cm) Co nc en tra tio n ( gc m3 ) 0 1 2 3 4 5 6 7 0 5 10 15 20 25 30 35 40 45 50 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 time,hrs Erro r 0 10 20 30 40 50 60 70 80 90 100 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Depth(z),cm Error American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2017) Volume 38, No 1, pp 235-249 248 Figure 3.1 and 3.2 compare the analytical and numerical simulation of nitrate transport in soil with respect to depth and time respectively. Figure 3.3 and Figure 3.4 show three dimensional analytical and numerical simulation of nitrate transport with respect to both depth and time respectively. Figure 3.5 and 3.6 depict the error associated with the numerical simulation with respect to time and depth respectively. The simulations show that although there are errors associated the finite difference numerical method, its results do not deviate sharply from the Laplace transform analytical simulations results. The errors in the finite difference simulations were as a result of truncation of Taylor series and therefore the error is mainly truncation error. The finite difference numerical simulation generated values which were greater than the values generated by the Laplace transform analytical simulation. The errors in both the temporal and spatial simulations initially grew to a point and started declining from that point. 4. Conclusion Modeling of solute transport in soil is a very complex process, hence the analytical models cannot capture most of the complexities of such process. Analytical models can only describe simplified process which is far from what pertains in reality. Numerical models are able to capture most the complexities of solute transport in soil hence it is very close to what happens in nature. However, numerical models unlike analytical models provide approximate and not exact solutions. There are errors associated with the numerical models. A comparison between Laplace transform model and finite difference model for nitrate solute transport in soil showed that the finite difference model gave a good approximation of the Laplace transform model. There were errors associated with the numerical model but it was mainly as a result of truncation error. The error associated with both the temporal and spatial simulations grew initially to some point and decline from that point onwards. The finite difference numerical model can be used to model the complex solute transport process in soil. However, there will be some errors but these errors will not cause too much difference between the numerical model solution and analytical model solution should it have been possible. References [1] Ann McCauley and Clain Jones (2005), Water and Solute Transport in Soils, Soil & Water Management, Module 4, pp. 1, Montana State University-Bozema. [2] Pimentel, D., H. Acquay, M. Biltonen, P. Rice, M. Silva, J. Nelson, V. Lipner, S. Giordano, A. Horowitz, M. D’Amore (1992), Environmental and economic costs of pesticide use. BioScience. 42:750-760. [3] Pampolino M. F. (2000), Measurement of nutrient availability and solute transport in soils with ion exchange resin capsules, Faculty of Agriculture, Hokkaido University, Vol. 70, Part. I, pp. 41-88. [4] Marshall T. J. and Holmes J. W. (1988), Soil Physics, 2𝑎𝑎𝑑𝑑 edition, Cambridge University Press, Cambridge, USA. [5] Abramowitz, M. and Stegun, I. A. (Eds.) (1972), "Laplace Transforms." Ch. 29 in Handbook of American Scientific Research Journal for Engineering, Technology, and Sciences (ASRJETS) (2017) Volume 38, No 1, pp 235-249 249 Mathematical Functions with Formulas, Graphs, and Mathematical Tables, 9th Printing, pp. 1019-1030, New York: Dover. [6] Contaminant Transport in the Unsaturated Zone Theory and Modeling (http//:www.pc- progress.com/Documents/Jirka/Handbook_of_GE_P [7] Or D. and Jon M. Wraith (2002), Unsaturated Flow and Transport, Course Notes, Utah State University (Logan), Montana State University (Bozeman). [8] Wexler E. J. (1992), Review of analytical solution for one, two, and three dimensional solute transport in groundwater systems with uniform flow, Techniques of Water Resources Investigations of the United States Geological Survey, Denver, Colorado, USA. [9] Lal R. and Shukla M. K. (2004), Principles of Soil Physics, Ohio State University, Marcel Dekker, Inc, New York, USA. [10] Jury, W.A. and D. R. Nielson. (1989), Nitrate transport and leaching mechanisms. In R.F. Follett (ed.), Nitrogen Management and Ground Water Protection (pp. 139-157). Elsevier Sci. Pubs. Amsterdam. [11] Kumar A. Kumar D., Kumar N. (2009), Analytical modeling of one-dimensional advection diffusion equation with variable coefficients in a finite domain, Journal of Earth System Science, Vol. 118, No. 5, pp. 539-549. [12] Toride, N. Leij F. J. and van Genuchten M., Th. (1995), The CXTFIT code for estimating transport parameters from laboratory or field tracer experiment, Research Report, No. 137, US department of agriculture, Riverside California, USA. [13] Stroud K. A. and Dexter J. B. (2007), Engineering Mathematics, Palgrave Macmillan, 6th edition, New York, USA. [14] Stroud K. A. and Dexter J. B. (2011), Advanced Engineering Mathematics, Palgrave Macmillan, 5th edition, New York, USA. [15] Carslaw, H.S. and Jaeger, J.C. (1959), Conduction of heat in solids 2nd edition, Clarendon Press, Oxford. 1.1 Laplace Transform analytical method 2.1 Solute Transport Equation 2.1.1 Convection-Dispersion Equation