WATER RESOURCES RESEARCH, VOL. 30, NO. 11, PAGES 3233-3235,NOVEMBER 1994 Commenton "A space-timeaccuratemethodfor solvingsolute transportproblems"by S. G. Li, F. Ruan, and D. McLaughlin Christopher Zoppou Australian CapitalTerritoryElectricityandWater,Canberra JohnH. Knight Centre for EnvironmentalMechanics,CSIRO, Canberra in which c(x, t) is the solute concentration,D(x) is the diffusioncoefficient,K(x) is the first-orderdecaycoefficient, Analyticaltechniquesfor the solutionof the advection- u(x) is the fluidvelocity,S(x, t) is the sink/sourceterm, x is diffusion equationare generallyrestrictedto simpleprob- the spatialcoordinate (0 -< x -< i), t is the time(t > 0), and lemswith constantcoefficients.Practicalproblemsusually I is the lengthof the computationaldomain. involve variable velocity and diffusion coefficients. The It will be shownthat this equationis, with certainrestricsolutionschemedescribed by Li et al. [1992], based on the tions,mathematicallyequivalentto the conservativeform of Laplacetransform,was developedwith the capabilityto the advection-diffusionequation. solvethe advection-diffusionequationwith spatiallyvariable Considerthe advection-diffusionequationwritten in concoefficients. The Laplace transform is used to evaluate the servative form temporal derivativein the advection-diffusion equationanalytically,therebyeliminatingthe effectsof the time deriva+ --(x) ...... o (2) tiveon accuracyand stability. Becausethe temporalderivOt Ox Ox Ox ativeis evaluatedusing the Laplace transform, there is no needfor the time-stepping that is associatedwith more where the sink/source and reaction terms have been netraditionaltechniques for the solution of the advection- glected. diffusion equation.Therefore it is a potentially more efficient An analytical solution to (2) will be given for the case technique. when the velocity, u(x), is taken to be a linear function of Theauthorsdo not provide numericalsimulationsinvolv- distance,suchthat u(x) = UoX,andthe diffusioncoefficient ing the spatially variable coefficient advection-diffusion a quadratic function of distance, D(x) = Dox2, whereDo Introduction Oc(x,t) O<c(x,t)u(x)• OtD Oc(x,t) ) equation.A reasonfor this might be that there are very few analyticalsolutionsto the advection-diffusionequation with spatiallyvariable coefficients. The purpose of this discussionis to present what we believeto be a new analyticalsolutionto a particularform of and u0 are constantsand 1 < x < I. Equation (2) now becomes Oc(x, t) •+ Ot theadvection-diffusion equationwith spatiallyvariablecoefficients and to highlighta potential mass conservation Expanding, O(c(x, t)x) uo Ox 0( Oc(x,t).)= 0 -OOTx x problem that could be encountered with the formulation Oc(x, t) usedby Li et al. The analytical solutioncan be used to Oc(x, t) + (UoX- 2D0x) validateschemesfor solvingthe advection-diffusion equa- • at Ox tionwithspatiallyvariablecoefficients. The efficiencyof the method proposedby Li et al. is demonstratedby comparing O2c(x, t) = D0x2 • Ox itsperformancefor the solutionof the advection-diffusion 2 - UoC( X,t) (3) equation with spatiallyvariablecoefficientswith that of a well-known time-stepping scheme. The last term in this equation is essentialfor the conserva- tion of mass.The additionalterm, 2Dxo, in the advection ConservativeForm of the Advection-Diffusion term only affectsthe distributionof the mass. Comparing(3) with (1) reveals that the form of (1) is sufficiently general,as it accommodatesdifferentcoefficients Li et al. describea numericalscheme,based on the (with differing physicalinterpretation)of (3). It follows from Laplace transform, for the solution of thefollowing formof thespatiallyvariablecoefficientadvection-diffusion equa- the derivationof (3) that any expressioncanbe usedfor u(x) and D(x) providedthat the concentrationprofile is smooth tion: and differentiable.It is possibletherefore to ensurethat (1) conservesmass if in (1) Equation Oc(x, t) Oc(x, t) 0( Oc(x, t).) --------+ Ot u(x) Ox =-- Ox D(x) -K(x)c(x, Ox t)+ S(x, t) Published in1994 bytheAmerican Geophysical Union. Paper number 94WR01492. Ou K( x) = -+ Kp(x) Ox (•) whereKp(x) is the first-orderdecaycoefficient of the physicalprocessbeingmodeled. 3233 3234 ZOPPOU AND KNIGHT: COMMENTARY If, in (1), the reaction and sink/sourceterms were ne- 120. • ........--•.ta• (•,• 7), •,t6 o glected,a naivemodelermightuse Oc(x,t) oc(x,t) o2c(x,t) Ot (UoX' Ox=DOx2 • Ox -------+ 2Dox) 2 (4) The majordifferencebetweenthis equationand (3) is the omission of the lasttermin (3), whichis necessary for the loo. ••,ml• •o. conservation of mass. This is a nonconservative form of the advection-diffusion equation,and caremustbe exercisedin 40. implementing the algorithmproposed by Li et al. so that conservation of mass is not violated. 20. 0• 'AnalyticalSolutionto a SpatiallyVariable Advection-Diffusion Equation . ß0 5. 10. 15. 20. Ananalytical solution to (3)canbeobtained using Laplace 1. Concentration profilefor the spatiallyvariable transforms. Thisis achievedwith the changeof variables, Figure advection-diffusion coefficient. y = In (x) (see,for example,Hildebrand[1962,p. 13]). Equation (3) becomes oc(x, t) •+ Ot Oc(x,t) o2c(x,t) (uo- Do)•= Do•Oy Oy2 UoC(X, t) Alsoillustrated inFigure1arethesimulated profiles using (3),(4),(5), and(6). Thefiniteanalytic/Laplace timemethod (5) proposed by Li et al. was usedto solve(3) and(4).In contrast to Li eta!., whousedthewell-known Crump[1976] whichisin theformof anadvection-diffusion equation with algorithm,the more robustnumericalLaplaceinversion constant coefficientsand a first-order reaction term. Condevelopedby de Hoog et al. [1982]was used. sidera problemwith the following initialandboundary Theresultsfor the numericalinversionof (6) areindistinconditions: c(x, 0) = 0, c(1, t) = Coandc(o•,t) = 0, which guishablefrom the exact results, which indicate that the is similarto the problem usedby Li et al. TheLaplace concentration profiledecaysexponentially.The resultsobtained forthesolution of(3)usingthefiniteanalytic/Laplace time scheme are more accurate than the results obtained d•(y, s) d26(y,s) (s + Uo)•(y, s) + (Uo- Do) usingthe well-knownthird-orderHolly and Preissmann = [1977]scheme, whichis a time-stepping scheme. To satisfy in whichif(y, s) is the Laplacetransformof the solute theCourantcriterion,500time stepswererequiredforthe of (5).Thisscheme required approximately 10times concentration in theLaplacespaces. Solving for 5(y, s), solution the computational time requiredby the finite analytic/ transformof (5) is •(y, s) =- Laplacetimescheme for the solutionof (3). co((uo -Do) y s exp 2D 0 Conclusions D•/2 4Dø + s+ u0 (6) The solution of the nonconservative form of the advec- tion-diffusion equation,equation(4), producederroneous Performing theinverse Laplacetransform, (3) thenhasthe results.Massis not conserved,and the predictedprofile followinganalyticalsolution bearsno resemblance to the analyticalsolution. The solutionmethodproposedby Li et al. doesnot provide anysignificant advantages overexisting methods for the solution of the advection-diffusion equationwithconstantcoefficients. However,forproblems involving spatially variable coefficients, it hassignificant advantages overconschemes providedcareis exercised +exp Do cftc• 2(Dot) 1/2 (7) ventionaltime-stepping c(x, t)=-•- erfc 2(D0t) •/2 (Uo In (x))[ln(x)+t(uo+Do)]) to ensure conservation of mass. in whicherfcis thecomplementary errorfunction. Theanalyticalsolutionpresented canbe usedto validate numerical schemes for solvingthe advection-diffusion equa- Hypothetical Example tionwithvariable coefficients. Theanalytical solution canbe readilyextended to multidimensional problems. Consider anexample whereu0 = 1.0,D0 = 0.005,1 _< x _<I = 20,andCo= 100,which corresponds toa range of References PecletnumbersbetweenPe = 40 and 2, wherePe = K. S.,Numerical inversion of Laplace transforms using a uoAx/D o andAx = 0.2. Theexactconcentration profile, Crump, givenby (7) at t = 2.5 for thisproblem,is illustrated in Figure 1. Fourier series approximation, j. Assoc. Cornput. Mach.,23(1), 89-96, 1976. de Hoog,F. R., J. H. Knight,andA. N. Stokes,An improved ZOPPOUAND KNIGHT: COMMENTARY 3235 method for numerical inversionof Laplacetransforms, SIAM J. J. H. Knight(corresponding author),CSIRO Centrefor EnvironI-lildebrand, F. B., AdvancedCalculus for Applications, Prentice- mentalMechanics,GPO Box 821, Canberra,ACT 2601, Australia. (email:kni•ht@enmech.csiro.au) HaH,EnglewoodCliffs, N.J., 1962. $ci. Star. Cornput.,3(3), 357-366, 1982. C. Zoppou,ACT ElectricityandWater, GPO Box 366,Canberra, Holly, F. M., andA. Preissmann, Accurate calculation oftransport ACT 2601, Australia. intwo-dimensions, J. Hydraul.Eng.Div. Am. Soc.Civ.Eng., IO$(HYII), 1259-1277,1977. Li, S. G., F. Ruan, and D. McLaughlin,A space-time accurate method for solvingsolutetransportproblems,Water Resour. Res.,28(9), 2297-2306, 1992. (ReceivedNovember2, 1993;accepted June8, 1994.)