AnalysisofaDirtyMotorPoolRichardJardine,MichaelKelley,JosephMyersAtopic of enduring concern to militaryoperators and logisticians is theenvironmental signature left bybothfield and garrison operations.Thispaperdetailsoneofthescenarioswhichweuse inourengineeringmathcoursetoexciteandmotivatefutureenvironmentalengineersaboutusesofmathematicsintheirdiscipline.Thisscenarioillustratestheuseofvectorcalculus,partial differential equations(PDE's),andnumerical methods,alongwithacomputeralgebra system (MathCad)and spreadsheet (Excel),inmodeling the advection and diffusion of an oil spill.Workingthroughthis scenarioexercisesthefollowingmathematical skills:Parameterization of space curves,using vector differential operators,and using the vector integraltheorems- Modeling using PDE's-Solvingthediffusion equationviaseparationofvariables-Making engineering value judgments (such as how much seepage is"appreciable",vs.nonzero)When and how to use the dominant eigenmode as a long-termapproximation-Refining amodel to incorporate neweffects-Numerically (via finite differences and a spreadsheet) solving morecomplexvariantsofthediffusionequation- Graphing solution curves and drawing inferencesScenarioYouare a Battalion ExecutiveOfficerstationed in Korea.Amongyourmanyduties,you are in chargeof vehiclemaintenanceand motor pool operations forthebattalion.TheAssistantDivisionCommanderforSupport (ADC(S))isflyingoveryourbattalionmotorpool areaonemorning,andhappenstonoticesomegrounddiscolorationnear theBrigadePOL (Petroleum, Oils,and Lubricants)TankFarm,forwhichyouhaveprimaryresponsibility.Helater callsyourcommanderandyouandtellsyoutoinvestigate.Yourinitial checkshowsthatone of the fuel storage tanks has started to leak at a seam. It is a large-areaground oil spillwith continued leakage(maintaining the surface at a constantlevelof contaminantconcentration).Thecontaminationisbeingadvectedbysurfacerunoff andis alsodiffusingdownwardtowardbedrocklevel.TheFacilityEngineersestimatethat itwilltake72hourstorepairthepipeandstoptheleak.ThelocalcivilianEnvironmentalOfficerimmediatelycallsanddemandstoknow41
41 Analysis of a Dirty Motor Pool Richard Jardine, Michael Kelley, Joseph Myers A topic of enduring concern to military operators and logisticians is the environmental signature left by both field and garrison operations. This paper details one of the scenarios which we use in our engineering math course to excite and motivate future environmental engineers about uses of mathematics in their discipline. This scenario illustrates the use of vector calculus, partial differential equations (PDE’s), and numerical methods, along with a computer algebra system (MathCad) and spreadsheet (Excel), in modeling the advection and diffusion of an oil spill. Working through this scenario exercises the following mathematical skills: - Parameterization of space curves, using vector differential operators, and using the vector integral theorems - Modeling using PDE’s - Solving the diffusion equation via separation of variables - Making engineering value judgments (such as how much seepage is “appreciable”, vs. nonzero) - When and how to use the dominant eigenmode as a long-term approximation - Refining a model to incorporate new effects - Numerically (via finite differences and a spreadsheet) solving more complex variants of the diffusion equation - Graphing solution curves and drawing inferences Scenario You are a Battalion Executive Officer stationed in Korea. Among your many duties, you are in charge of vehicle maintenance and motor pool operations for the battalion. The Assistant Division Commander for Support (ADC(S)) is flying over your battalion motor pool area one morning, and happens to notice some ground discoloration near the Brigade POL (Petroleum, Oils, and Lubricants) Tank Farm, for which you have primary responsibility. He later calls your commander and you and tells you to investigate. Your initial check shows that one of the fuel storage tanks has started to leak at a seam. It is a large-area ground oil spill with continued leakage (maintaining the surface at a constant level of contaminant concentration). The contamination is being advected by surface runoff and is also diffusing downward toward bedrock level. The Facility Engineers estimate that it will take 72 hours to repair the pipe and stop the leak. The local civilian Environmental Officer immediately calls and demands to know
whererunoffwill carrythecontamination,towhatdepthcontaminatedsoilwillhavetobeexcavated after theleak hasbeen stopped,and whatthe effect wouldbeonthespreadof contaminant iftheweatherisrainyoverthenextweek.TrackingAvectiveProcessesand RunoffYou decide to analyze the quickest acting process first: runoff.A quick survey ofthe terrain shows that the major runoff for the contaminant appears to be thestreamflowing from nearthespill (pointAon themapbelow)toward pointBneartheneighboringtown.Thevertical grid lines are1000mapart,thehorizontalgrid lines are75omapart,andcontour linesareatelevationsasmarked onthemapvelFigure1:Pathof streamfromatPointAtoB.Example 1:Mathematicallymodel thepath in orderto determine theflow of thestream.Solution:Webeginby using Mathcad toplot amodel the stream.We writeaparametricvectorfunctionforthepositionofamarkerparticle inthe streamas itmoves from A toB.A possibleparameterization of the streambed istheparametric equation r(t)=(e2t sin 7t+ sin15t)i+800t j+130t k, 0≤t≤1 (where t42
42 where runoff will carry the contamination, to what depth contaminated soil will have to be excavated after the leak has been stopped, and what the effect would be on the spread of contaminant if the weather is rainy over the next week. Tracking Avective Processes and Runoff You decide to analyze the quickest acting process first: runoff. A quick survey of the terrain shows that the major runoff for the contaminant appears to be the stream flowing from near the spill (point A on the map below) toward point B, near the neighboring town. The vertical grid lines are 1000m apart, the horizontal grid lines are 750 m apart, and contour lines are at elevations as marked on the map. Figure 1: Path of stream from at Point A to B. Example 1: Mathematically model the path in order to determine the flow of the stream. Solution: We begin by using Mathcad to plot a model the stream. We write a parametric vector function for the position of a marker particle in the stream as it moves from A to B. A possible parameterization of the streambed is the parametric equation ( ) ( sin 7 sin15 ) 800 130 , 0 1 2 t = e t + t + t + t £t £ t r i j k (where t
is measured in seconds and x,y,= are measured in meters).As graphicallyillustrated in Figure 2, this yields a good approximation to the path of the stream.Stream model0180-20060020X,Y,7Streammodel10050C200400600800X,Y,ZFigure2:ParameterizedrepresentationoftheAdvectingStreamExample2:Giventheabovepositionfunctionof the stream,withAatt=1andBat t=o,youmeasureand approximatethat thewater's velocityfield is infactv=0.li-0.2 j-0.1k.Find theflowof water(as defined in Finney/Thomas,Revised Printing,p.952)moving in theone-dimensional model of thestreambetweenAand B.Describewhat this meansphysically,includingthe signof theresult.43
43 is measured in seconds and x, y,z are measured in meters). As graphically illustrated in Figure 2, this yields a good approximation to the path of the stream. Figure 2: Parameterized representation of the Advecting Stream Example 2: Given the above position function of the stream, with A at t = 1 and B at t = 0 , you measure and approximate that the water's velocity field is in fact v = 0.1i - 0.2 j - 0.1k . Find the flow of water (as defined in Finney/Thomas, Revised Printing, p.952) moving in the one-dimensional model of the stream between A and B. Describe what this means physically, including the sign of the result. X,Y,Z Stream model 4 2 0 2 4 0 200 400 600 800 1005 X,Y,Z Stream model 4 2 0 2 4 0 200 400 600 800 50 100
Solution:Theflowof theadvectingstreamis calculated via the line integralJvdr for the given velocity field, yielding:CpJvdr = pfbv(t)-r(t)dt =-215.6m/ s.This isacombinedmeasurethatindicatesbothhowfarandhowfastthefluid ismoving.It isnegative,whichmeansthatthenetflowvelocityis inthedirectionfromAtoward B.The magnitude, obtained by dividing the total flowbythelength of the stream, is of about the correct order, as the Colorado River has aflowoffrom300to1000m/s.Wenowseek todeterminethe effect ofvorticity on contaminantmonitoringequipmentplaced inboreholesthat wedigatseveralspecified locationsalongthe length of the stream. Several boreholes have been dug near the spill siteand instruments placed in the wells to obtain data on the effect ofthe spill. (Seefiqurebelow).Eachholeiscylindricallyshaped,2metersdeepand20cmindiameter.In one oftheholes,the velocityfield ismeasured tobev=xi+yj-3zkm/s (withorigincenteredatthebottomof thehole)LandfillTestwellsKVAquiferFigure 3: Bore holes for contaminant measuring equipment.Example3:Atwhatrateiswaterenteringtheholefromabove?Solution:The flux of water into the hole is given byJJ v·ndA = J-kdA = JJ- 3zdA = -3.2[JdA = -6A= -0.1885m3 /s4The sign isnegative sincenetflux is opposite in directiontothe outward normalvector; i.e., netflux is downward throughthe upper surface.Example4:Are there any sources or sinks in the interior of the hole? Describe.44
44 Solution: The flow of the advecting stream is calculated via the line integral ò × C v dr for the given velocity field, yielding: ( ) ( ) 215.6 / . 1 0 d t t dt m s C r òv ×r = rò v ×r¢ = - This is a combined measure that indicates both how far and how fast the fluid is moving. It is negative, which means that the net flow velocity is in the direction from A toward B. The magnitude, obtained by dividing the total flow by the length of the stream, is of about the correct order, as the Colorado River has a flow of from 300 to 1000 m3 /s. We now seek to determine the effect of vorticity on contaminant monitoring equipment placed in bore holes that we dig at several specified locations along the length of the stream. Several boreholes have been dug near the spill site and instruments placed in the wells to obtain data on the effect of the spill. (See figure below). Each hole is cylindrically shaped, 2 meters deep and 20 cm in diameter. In one of the holes, the velocity field is measured to be v = xi + yj - 3zk m/s (with origin centered at the bottom of the hole). Landfill Test wells Aquifer Figure 3: Bore holes for contaminant measuring equipment. Example 3: At what rate is water entering the hole from above? Solution: The flux of water into the hole is given by òò × = òò × = òò- = - × òò = - = - A A A A dA x y z dA zdA dA A m s 3 v n , , 3 k 3 3 2 6 0.1885 . The sign is negative since net flux is opposite in direction to the outward normal vector; i.e., net flux is downward through the upper surface. Example 4: Are there any sources or sinks in the interior of the hole? Describe
Solution:Fromthevelocityfield of thefluid intheborehole,we calculatethatthe divergenceatall points is V·v=1+1-3=-1 per sec.Since thedivergenceis negative at all points, this indicates that there is a sink at every point in thehole.Thismightoccurif thedensityofthefluiddecreasedafterit entered thehole,perhapsduetowarmingbythe equipment inthehole.Example5:Atwhat rate is water entering thehole from all sides (including theporoussidesand bottom,aswell asthetop)?Solution:Thetotalflux offluid intotheholecanbecalculated bytransformingtheclosedsurface integral v·ndo intothevolumeintegral JjVvdv viaTStoke's Theorem, yielding fff-ldV =-1.V ; this gives a net flux into the hole of-.06275m2 / s.In another of thetest wells,ground waterrotates aroundthe center of the hole(within the plane z = 1) with velocity v= o(yi - xi), where , the angularvelocity,is 0.25rad/s.If the circulation in the well exceeds 15 m2/s,the datacollected by instruments inthewell will becorrupted.Example6:Will theempirical results obtained fromthatholebevalid?Solution:Theflow circulation around the midpointplaneof the cylindrical holecanbefound bytransformingtheline integral v·drintothedoubleintegralJJV×vdV via Gauss'Theorem. This yields:AJ·kdA=-20A=-0.0157m2/sThe circulation ismuch lessthanthe critical valueoftheequipment;therefore,collected datawillbereliable.Calculating DownwardDiffusionYou now turn your attention to the contamination of the ground underneath thespill duetodiffusion. Youbeginwiththefollowing assumptions:The spill is uniform over an area large enough such thatthere isdiffusion only in the vertical direction.-Thesoil ishomogeneouswithadiffusionrateof0.3ft/hr.-Thereisnoadvection.- The leak is such that it maintains a constant concentration of 5000 g/ft3overthesurface45
45 Solution: From the velocity field of the fluid in the borehole, we calculate that the divergence at all points is Ñ ×v = 1+ 1- 3 = - 1 per sec. Since the divergence is negative at all points, this indicates that there is a sink at every point in the hole. This might occur if the density of the fluid decreased after it entered the hole, perhaps due to warming by the equipment in the hole. Example 5: At what rate is water entering the hole from all sides (including the porous sides and bottom, as well as the top)? Solution: The total flux of fluid into the hole can be calculated by transforming the closed surface integral òò × s v nds into the volume integral òòòÑ × V vdV via Stoke’s Theorem, yielding dV V V òòò- 1 = - 1× ; this gives a net flux into the hole of - .06275 / 3 m s. In another of the test wells, ground water rotates around the center of the hole (within the plane z = 1) with velocity v = w ( yi - xj) , where w , the angular velocity, is 0.25 rad/s. If the circulation in the well exceeds 15 m2 /s, the data collected by instruments in the well will be corrupted. Example 6: Will the empirical results obtained from that hole be valid? Solution: The flow circulation around the midpoint plane of the cylindrical hole can be found by transforming the line integral ò × C v dr into the double integral òòÑ ´ A vdV via Gauss’ Theorem. This yields: òò × = - = - A 0,0, 2w kdA 2wA 0.0157m 2 /s. The circulation is much less than the critical value of the equipment; therefore, collected data will be reliable. Calculating Downward Diffusion You now turn your attention to the contamination of the ground underneath the spill due to diffusion. You begin with the following assumptions: - The spill is uniform over an area large enough such that there is diffusion only in the vertical direction. - The soil is homogeneous with a diffusion rate of 0.3 ft2 /hr. - There is no advection. - The leak is such that it maintains a constant concentration of 5000 g/ft3 over the surface
- There is bedrock beginning at a depth of 20 ft which is impervious tofuel.Example7:Determinehowmuchsoilmustbeexcavated.Solution: We begin by modeling the contaminant diffusion by the equationou/b,= 0.30*u/,2 , which we will hereafter abbreviate as u, = 0.3u , where u is0Tatthelocal concentration of contaminant,alongwith boundary conditionsu(0,t)= 5000 (constant concentration in the spill area), u,(20,t)=0(impermeableatthepermafrostlevel),andinitial condition u(x,O)=o (initiallyclean).Separationofvariables yields theinfinite series solutionu(x,t)=5000-20000/元Z1/2n-1Sin(2n-1)x/40)exp(-0.3(2n-1)元" /40)Plots of contaminant concentrationatfixed depths overtime and atfixed timesasfunctionsofdepthareshowninFigure4.SolutionatFixedDepthsoverTimeSolutionatFixed TimesasaFunction ofDepth:t2 :=0,.1..72x2:=0,0.1..20400040005,t2,50){x2,18,50)10,t2,50){x2,36,50)20002000{x2,54,50)15,t2,50)20,t2,50)(x2,72,50)105002012x2Figure 4: Contaminant concentration at fixed depths and fixed times.Theanswertothequestion“howmuchtoexcavate?"is:itdependsonhowlongyouwaitandwhatconcentrationofcontaminantisacceptable.Forexample,waiting72hoursmeansthatwewill havetoexcavatetobedrocklevel (20ft)ifwecannottolerateacontaminant concentrationgreaterthan5oog/m.Whatlevelofcontaminantconcentrationisunacceptable?Ifweareinterestedonlyinthephysical effectsofcontamination,thenchemistryandphysiologycansometimesgiveusananswer;thosedisciplinescanhelpusestimateatwhatlevelwill thecontaminant showupasamarginal healthhazard inourlocal46
46 - There is bedrock beginning at a depth of 20 ft which is impervious to fuel. Example 7: Determine how much soil must be excavated. Solution: We begin by modeling the contaminant diffusion by the equation 2 2 0.3 t u t u ¶ = ¶ ¶ ¶ , which we will hereafter abbreviate as t xx u = 0.3u , where u is the local concentration of contaminant, along with boundary conditions u(0,t) = 5000 (constant concentration in the spill area), u t x (20, ) = 0 (impermeable at the permafrost level), and initial condition u(x,0) = 0 (initially clean). Separation of variables yields the infinite series solution u x t n n ( , ) = - - = ¥ 5000 20000 å 1 2 1 1 p Sin((2n - 1)px 40) exp ( 0.3(2 1) 40) 2 2t - n - p . Plots of contaminant concentration at fixed depths over time and at fixed times as functions of depth are shown in Figure 4. Figure 4: Contaminant concentration at fixed depths and fixed times. The answer to the question “how much to excavate?” is: it depends on how long you wait and what concentration of contaminant is acceptable. For example, waiting 72 hours means that we will have to excavate to bedrock level (20 ft) if we cannot tolerate a contaminant concentration greater than 500 g/m3 . What level of contaminant concentration is unacceptable? If we are interested only in the physical effects of contamination, then chemistry and physiology can sometimes give us an answer; those disciplines can help us estimate at what level will the contaminant show up as a marginal health hazard in our local Solution at Fixed Depths over Time: Solution at Fixed Times as a Function of Depth: t2 0,.1. 72 x2 0,0.1. 20 f( 5,t2,50) f( 10,t2,50) f( 15,t2,50) f( 20,t2,50) t2 0 50 0 2000 4000 f( x2,18,50) f( x2,36,50) f( x2,54,50) f( x2,72,50) x2 0 10 20 0 2000 4000
drinking water orintheaquiferforthelarger community,andthat is where thescientistwould beinclinedtodrawtheline.However, fromapsychologicalstandpoint,local residentsmaytakedecadestoacceptthefactthattheirhomessitandtheirchildrenplayinthedirtabovecontaminants,eveniftheconcentrationiswellbelowthephysiologicalthreshold.Economically,propertyvalues may suffer for decades forthe same reason. Legally,owners anddevelopersmaycarrytheadditional riskofbeingheld liableforrandomhealthproblems,evenifnotclearlyrelatedtothepresenceofcontaminantsObviously,thedecisionastowhatconstitutesa"tolerable"concentrationgreatlyeffects the results toour question, andis asmuch (usuallymore)apoliticaldecision as it is a scientific one. (An important learning point for engineers!)Inmanysituations involvingdiffusion(suchas underlandfills),environmentalengineers areoftenconcernedonlywithsolutions atlargetimes (oftenontheorderofdecades),andnotwithtransientsovershorttimescales.Wecanderivesuchasimplerlongtermsolutiontooursituationabovebyrecognizingthathigherordermodesbecome increasinglyless importantovertime,yielding(forexample) the much simpler two term approximation for the contaminantconcentration at long times:u(x,t)=5000-20000/元Sin(元x/40)exp(-0.3元2t/1600)CalculatingDownwardDiffusionwithAdvectionWehave analyzed the purediffusion case;now we have to deal with theadditiontothemodelofanotherphysicalprocess:namely,advectionofthecontaminantdownwardduetorainfall.Example8:Determinewhetherrainfall will haveanyeffect onourresultsonthedownwarddiffusionofcontaminant.Solution: Analytic techniques appear harder in this case, so we turn to anumericalmethod.Topredictthedifferencesthatwouldbecausedbyrainwateradvection (downward)ifwehaverain,weconsideramodel that incorporatesbothdiffusionandadvection:u,=0.3u-0.0lux,alongwiththepreviousboundaryand initial conditions.Wedifferencetheequation using forward first differences in time and centeredfirstandseconddifferencesinspace,resultinginthealgorithm(whichweabbreviateastheFTCSalgorithm,forforward-time,centered-space):uj.kI = uj.x + 0.3N/x2 (uj-k -2ujx +uj-I,)- 0.01/2Ar (u j+1, - uj-1A). (1)八47
47 drinking water or in the aquifer for the larger community, and that is where the scientist would be inclined to draw the line. However, from a psychological standpoint, local residents may take decades to accept the fact that their homes sit and their children play in the dirt above contaminants, even if the concentration is well below the physiological threshold. Economically, property values may suffer for decades for the same reason. Legally, owners and developers may carry the additional risk of being held liable for random health problems, even if not clearly related to the presence of contaminants. Obviously, the decision as to what constitutes a “tolerable” concentration greatly effects the results to our question, and is as much (usually more) a political decision as it is a scientific one. (An important learning point for engineers!) In many situations involving diffusion (such as under land fills), environmental engineers are often concerned only with solutions at large times (often on the order of decades), and not with transients over short time scales. We can derive such a simpler long term solution to our situation above by recognizing that higher order modes become increasingly less important over time, yielding (for example) the much simpler two term approximation for the contaminant concentration at long times: u(x,t) = 5000 - 20000 pSin(px 40 ) exp( 0.3 /1600 2 - p t ). Calculating Downward Diffusion with Advection We have analyzed the pure diffusion case; now we have to deal with the addition to the model of another physical process; namely, advection of the contaminant downward due to rainfall. Example 8: Determine whether rainfall will have any effect on our results on the downward diffusion of contaminant. Solution: Analytic techniques appear harder in this case, so we turn to a numerical method. To predict the differences that would be caused by rainwater advection (downward) if we have rain, we consider a model that incorporates both diffusion and advection: t xx x u = 0.3u - 0.01u , along with the previous boundary and initial conditions. We difference the equation using forward first differences in time and centered first and second differences in space, resulting in the algorithm (which we abbreviate as the FTCS algorithm, for forward-time, centered-space): ( ) 2 0.3 ( 2 ) 0.01 j,k 1 j,k 2 j 1,k j,k j 1,k j 1,k j 1,k u u x t u u u x t u u + - - + - - D - + - D D = + D . (1)
Weimplementthis ina spreadsheet,asdemonstrated inthefigurebelow.Spatialgridpointsarenumberedhorizontallyindarkshade,temporalgridpointsvertically indark shade, initial condition horizontally in light shade, boundaryconditionsverticallyinlightshade810121416182006N400000000000OolO0000ol05000o000500037500o000o05000000000693.7528.1250o05000966.7969000075.93752.10937500o0137.21487.4882810.158203C050001202.473O500000001407.393207.379716.668020.6960940.0118650o050001586.837283.077329.77351.8426710.0622920.000890ol050001745.043361.861546.676473.8039550.1912160.0054286.67E-05050000.0004645.01E-06001885.426441.961267.099916.7484380.4482370.018961050002010.75910.802280.8885560.0497690.0018173.9E-053.75E-073.75E-07522.106590.688151150002123.30316.050193.28E-06601.399117.05311.5691770.1090810.005280.0001693.28E-06125000 2224.91322.539331.57E-05679.2159145.80382.5457460.2108030.0126820.000541.57E-0550005.51E-05132317.117755.1373176.564930.284653.8701440.3710650.026630.0014125.51E-051450002401.185828.8928208.986839.274585.5888010.6076630.0505710.0032010.0001570.0001571550000.0003852478.174900.3217242.751349.476567.741650.9394670.0888010.0065260.0003851650002548.972969.3429277.573560.8420510.36161.3858310.146430.0122360.0008460.0008461750002614.3270.00171035.932313.201373.3108813.474451.9660590.229320.0214460.00171850002674.8731100.107349.414486.8149317.099062.6989330.3439850.0355560.0031810.0031811950002731.151161.913386.0214101281221.247743.6023210.4974740.0562600056090.005609Figure5:Aspreadsheet(Excel)implementationofthefinitedifference(FTCS)solution.Wecanthengraphicallycomparetheeffectsofdiffusionaloneversusdiffusionwithconvection.56789101112131415187181920Mymoattsgoeooateeeaea0Tme HourgDepth (t)Figure6:Comparisonofdiffusiveanddiffusive-advectivesolutions.48
48 We implement this in a spreadsheet, as demonstrated in the figure below. Spatial grid points are numbered horizontally in dark shade, temporal grid points vertically in dark shade, initial condition horizontally in light shade, boundary conditions vertically in light shade. Figure 5: A spreadsheet (Excel) implementation of the finite difference (FTCS) solution. We can then graphically compare the effects of diffusion alone versus diffusion with convection. Figure 6: Comparison of diffusive and diffusive-advective solutions. Concentration vs Depth 0 1000 2000 3000 4000 5000 6000 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 Depth (ft) 13 25 38 50 72 13 25 38 50 72 At Times: t=0 Without Advection: With Advection: Concentration vs Time 0 500 1000 1500 2000 2500 3000 Time (Hours) 5 10 15 20 5 10 15 20 At 4 Depths No Advection: With Advection: t / x 0 2 4 6 8 10 12 14 16 18 20 0 0 0 0 0 0 0 0 0 0 0 0 1 5000 0 0 0 0 0 0 0 0 0 0 2 5000 375 0 0 0 0 0 0 0 0 0 3 5000 693.75 28.125 0 0 0 0 0 0 0 0 4 5000 966.7969 75.9375 2.109375 0 0 0 0 0 0 0 5 5000 1202.473 137.2148 7.488281 0.158203 0 0 0 0 0 0 6 5000 1407.393 207.3797 16.66802 0.696094 0.011865 0 0 0 0 0 7 5000 1586.837 283.0773 29.7735 1.842671 0.062292 0.00089 0 0 0 0 8 5000 1745.043 361.8615 46.67647 3.803955 0.191216 0.005428 6.67E-05 0 0 0 9 5000 1885.426 441.9612 67.09991 6.748438 0.448237 0.01896 0.000464 5.01E-06 0 0 10 5000 2010.759 522.1065 90.68815 10.80228 0.888556 0.049769 0.001817 3.9E-05 3.75E-07 3.75E-07 11 5000 2123.303 601.399 117.0531 16.05019 1.569177 0.109081 0.00528 0.000169 3.28E-06 3.28E-06 12 5000 2224.913 679.2159 145.8038 22.53933 2.545746 0.210803 0.012682 0.00054 1.57E-05 1.57E-05 13 5000 2317.117 755.1373 176.5649 30.28465 3.870144 0.371065 0.02663 0.001412 5.51E-05 5.51E-05 14 5000 2401.185 828.8928 208.9868 39.27458 5.588801 0.607663 0.050571 0.003201 0.000157 0.000157 15 5000 2478.174 900.3217 242.7513 49.47656 7.74165 0.939467 0.088801 0.006526 0.000385 0.000385 16 5000 2548.972 969.3429 277.5735 60.84205 10.3616 1.385831 0.14643 0.012236 0.000846 0.000846 17 5000 2614.327 1035.932 313.2013 73.31088 13.47445 1.966059 0.22932 0.021446 0.0017 0.0017 18 5000 2674.873 1100.107 349.4144 86.81493 17.09906 2.698933 0.343985 0.035556 0.003181 0.003181 19 5000 2731.15 1161.913 386.0214 101.2812 21.24774 3.602321 0.497474 0.05626 0.005609 0.005609
Inthisscenario(withanadvectivecoefficientofo.o1).weseethatadvectionduetorainwater only slightly (byabout10%)speeds thedownwardmovementofcontaminant.ExtensionsGiven thenumerical (FTCS)procedure weused in the previous requirement, weare inapositionto investigateotherissues and complications.Convergenceissues are present in both our analytic and numerical solutions; the analyticsolution suffers fromtruncationerrors when calculatingfromtheinfiniteseriessolution,whereas the numerical solution suffersfrom thepropagation ofroundofferrors and possible instability (Figure7)AnalyticNumericalFigure7:Truncationerror(analyticsolution)andinstability(numerical solution)Ontheleftisplottedthetruncationerrorintheanalyticsolutionvs.numberoftermsretainedintheinfinitesum.Ontheupperrightisthenumerical solutionofEqn1atagivendepthasafunctionoftime;notetheapparentlyunboundedtemporalerroroftheunstablesolution.Atlowerrightisthenumerical solutionofEqn 1at a given timeas a function of depth; note that although thefunctionobeystheboundaryconditions,wehaveoscillations(reachingnonphysicalnegativevaluesfortheconcentration)that indicatespatial instability49
49 In this scenario (with an advective coefficient of 0.01), we see that advection due to rainwater only slightly (by about 10%) speeds the downward movement of contaminant. Extensions Given the numerical (FTCS) procedure we used in the previous requirement, we are in a position to investigate other issues and complications. Convergence issues are present in both our analytic and numerical solutions; the analytic solution suffers from truncation errors when calculating from the infinite series solution, whereas the numerical solution suffers from the propagation of roundoff errors and possible instability (Figure 7). Figure 7: Truncation error (analytic solution) and instability (numerical solution). On the left is plotted the truncation error in the analytic solution vs. number of terms retained in the infinite sum. On the upper right is the numerical solution of Eqn 1 at a given depth as a function of time; note the apparently unbounded temporal error of the unstable solution. At lower right is the numerical solution of Eqn 1 at a given time as a function of depth; note that although the function obeys the boundary conditions, we have oscillations (reaching nonphysical negative values for the concentration) that indicate spatial instability
Othercomplicationscanalsobeinvestigatednumerically,suchastheeffectofnonhomogeneousgroundporosity(turningtheequation intoavariablecoefficient PDE),the effectof boundaries in a small-area spill, the breakdown ofcontaminant over time dueto chemical reactions,and theeffect of slowing orstopping theleak earlierbefore cleanup canbegin (introducing a time dependentboundarycondition).The finitedifferences and the spreadsheet canbe usedas aboveto investigatethe introductory combined diffusion-advection scenario.This leads naturally intomorerealistic,messiergroundwaterflowproblemsandmoresophisticatednumerical codes that are typically used in later hydrogeology and environmentalengineeringcourses.Exercises1.Look up Darcy'slaw (found inmost environmental engineeringtexts,suchasreferences [1] or [2]).Describe how this law can extend the ideas of this projecttomoregeneral advectionproblemsinatmospheric,reservoir,andwatertablemodeling.2.Inourseparationofvariablessolutionforthepurediffusionscenarioabove,comparetheratioofthefirstandsecondtermsofthesumwhent=0,t=10,t=20t=30.Repeatfortheratioofthesecondandthirdterms.Whatisthetrend?What does this implyabout thelongterm solutions to the diffusion equation?3.Our separation ofvariables solution forthepurediffusion scenario expressedthecontaminantconcentrationasafunctionofdepthandtime.Insteadofplotting curves at constant depth and time,we could instead justplot thisequationasasurface;engineerscallthisaresponsesurface.Describehowwecould use theresponse surface tographicallygenerate the curves of Figure4.4.Youare thefacility engineerin Berlin.The city ofDresden calls and asksyourassistanceinmakingcleanupestimatesforanold (nowunoccupied)Soviettank companymotorpool.ThegroundhasbeensaturatedwithPOL,andthe citygovernment wants to excavate the contaminated soil and replace it withclean fill in order tobuild low cost housing on the site.Thetank companyoccupiedthesitefor35years,andappearstohaveroutinelydumpedwastePOLduringtheentireperiod.Giventhat theaverage local diffusivity of the soilis 1.75m'lyr,to what depth will the soil needtobe excavated?50
50 Other complications can also be investigated numerically, such as the effect of nonhomogeneous ground porosity (turning the equation into a variable coefficient PDE), the effect of boundaries in a small-area spill, the breakdown of contaminant over time due to chemical reactions, and the effect of slowing or stopping the leak earlier before cleanup can begin (introducing a time dependent boundary condition). The finite differences and the spreadsheet can be used as above to investigate the introductory combined diffusion-advection scenario. This leads naturally into more realistic, messier groundwater flow problems and more sophisticated numerical codes that are typically used in later hydrogeology and environmental engineering courses. Exercises 1. Look up Darcy’s law (found in most environmental engineering texts, such as references [1] or [2]). Describe how this law can extend the ideas of this project to more general advection problems in atmospheric, reservoir, and water table modeling. 2. In our separation of variables solution for the pure diffusion scenario above, compare the ratio of the first and second terms of the sum when t=0, t=10, t=20, t=30. Repeat for the ratio of the second and third terms. What is the trend? What does this imply about the long term solutions to the diffusion equation? 3. Our separation of variables solution for the pure diffusion scenario expressed the contaminant concentration as a function of depth and time. Instead of plotting curves at constant depth and time, we could instead just plot this equation as a surface; engineers call this a response surface. Describe how we could use the response surface to graphically generate the curves of Figure 4. 4. You are the facility engineer in Berlin. The city of Dresden calls and asks your assistance in making clean up estimates for an old (now unoccupied) Soviet tank company motor pool. The ground has been saturated with POL, and the city government wants to excavate the contaminated soil and replace it with clean fill in order to build low cost housing on the site. The tank company occupied the site for 35 years, and appears to have routinely dumped waste POL during the entire period. Given that the average local diffusivity of the soil is 1.75 m2 /yr, to what depth will the soil need to be excavated?