Home

Here

Contents

1. eq 1 6 5 Updated concentrations of ions Moa Msop eq 1 9 Figure 2 1 Flow chart for speciation calculations modified from typically report the total concentration of an ion which may be present in many different complexes e g Ca t 2 Cat Moat F Moaont F Maggot T Magog T Mapes t Moart Te 2 9 Note that these are dissolved species they are not to be confused with minerals Because Ca is tied up in these complexes the activity of the free calcium ion is reduced The calculation of the concentrations of all possible species is a complex task because 1 for each ion in solution a mass balances similar to 2 9 has to be solved and 2 at the same time the activities of all of the species have to obey the equilibrium relationships that follow from the law of mass action equation 2 2 This requires an iterative procedure figure 2 1 As a first estimate it is assumed that the total concentrations of the ions equal the concentrations of the free uncomplexed ions It is then possible to calculate the ionic strength step 1 and use that to find the activity coefficients and correct 17 Chapter 2 for electrostatic effects From this follow the activities of the ions step 2 that are inserted into the mass action equations to obtain the activities of the aqueous complexes step 3 These are converted to concentrations using the already calculated activity coefficients step 4 and finally the concentrations of the unc
2. where C s is the single species aqueous phase solubility available from tabulated sources e g of the organic compound in question within a mixture of compounds with different physico chemical properties 7 is the activity coefficient of the i organic compound typically assumed to be unity and m is the mole fraction of the i organic compound within the NAPL mixture It is defined as Ch oni 4 15 Mi where C is the molar concentration of compound i in the NAPL phase and C is the total molar concentration of all organic compounds in the NAPL phase A high groundwater flow velocity among other fac tors might result in a kinetically limited dissolution of NAPL com pounds The simplest model that describes the concentration change of the it compound in groundwater is V dis wi C at me Ci 4 16 where C is the concentration of the it organic compound in the groundwater and w is a mass transfer rate coefficient that is a product of a mass transfer coefficient and the specific interfacial area between NAPL phase and water The combination of 4 14 and 4 16 ap plies to arbitrary dissolution rates Note that w approaches infinity for equilibrium dissolution In that case C equals C3 For the modelling exercise we use Eqn 4 16 to compute NAPL dissolution The equation was translated into a rate expression for 63 Chapter 4 the PHREEQC PHT3D database We assume that the NAPL source consists
3. specific reaction database into the model folder and rename it to pht3d_datab dat e Import the database and set a very simple water composition that contains only Na Cl Benz dissolved benzene Tolu Chemical kinetics dissolved toluene Ethy dissolved ethylbenzene and Xyl dissolved xylene Include also the species Benznapl Tolu napl Ethynapl Xylnapl and Lowsolubnap which together represent the immobile NAPL source For both the ambient water background and the recharge water define concentrations of 0 001 mol l for Na and Cl respectively Set the initial values initial concentrations of pH to 7 and set the pe to 4 0 Define a NAPL polluted source area near the water table 20 m gt z gt 19 m between x 15m and x 35m Select an initial concentration of 0 018 mol l for the 4 NAPL compounds repre senting BTEX Enter an initial concentration of 0 072 mol l for the compound Lowsolubnap to all the grid cells that sit above a level of 19 m Note that the concentration of these immobile NAPL compounds is defined in units of mol l Define the settings for advective transport Use the MMOC scheme in this and other PHT3D simulations that include free water ta bles Define a longitudinal dispersivity of 0 50m and a transverse ver tical dispersivity of 5mm Run the simulation and visualise the results Define some observation points and visualise the breakthrough curves of benzene and xylene Wh
4. 0m and y 7m for which a value of 1 is allocated This defines that the concentrations at these grid cells will not change during the en tire simulation The concentrations that iPHT3D will allocate to the model s grid cells are the ones defined as initial concentrations This means that the water composition that was defined as background water composition will also be used as inflow water composition To additionally define that the acidic leachate composition is used in the upper section we need to select Spatial Attributes gt PHT3D gt PH gt ph 4 Source Sink and define a zone line for the uppermost 3 m y 7m to y 10m of the boundary the depth zone at which the acidic leachate enters and allocate a Zone Value of 1 This induces that the concentrations that were previously defined for Solution 1 in the Solutions Tab will be applied to the water that enters the model domain in the upper section of the inflow boundary Mineral dissolution and precipitation 3 2 5 Running PHT3D To run PHT3D proceed as earlier in the exercise when executing MOD FLOW e Go to Parameters Chemistry and use the Pencil Button to write all the input files required to run PHT3D e Then go to Parameters gt Chemistry and use the Plume Button to start PHT3D 3 2 6 Visualisation of results After the end of the simulation visualise the results in various ways e Obtain a contour plot by navigating in the Results panel of ipht3d Selec
5. e To add the exchanger species to the reaction network select go to Parameters Chemistry and acitivate X under the Exchange Tab Add a background value of 0 001 as a first estimate for the CEC e Now rerun PHT3D e Once the model run is complete reuse the python script to plot again the breakthrough curves BTCs for the ammonium Ca and Na concentrations and repeatedly compare the results with the observed data 47 Chapter 3 e Continue this manual calibration of the CEC by rerunning the model with successively improved estimates until a good agree ment is achieved 48 Mineral dissolution and precipitation 3 4 Surface complexation Oxides hydroxides and organic matter carry a surface charge which enables them to sorb ions from solution The sorption behaviour de pends on the type of mineral and crystal morphology as well as the composition and pH of the solution A potential develops due to the surface charge and therefore the equation for the Gibbs free energy AG contains besides a chemical term a Coulombic term to account for the work needed to move ions to or from the surface It can be shown that for a dissociation reaction Appelo and Postma 2005 zFyY RT n 10 in which K is the apparent dissociation constant Kint is the intrinsic dissociation constant z is the charge of the ion F the Faraday constant 96 485 C mol wo the potential at the surface R the gas constant 8 314 J K mol and T the abs
6. sition the molality and the number of moles of each element are reported These should be equal to the concentrations that were entered under the SOLUTION keyword although small deviations may occur when the mass balance equation 2 9 could not be solved exactly The next group of data is Description of solution which lists some characteristics of the solution for example pH and ionic strength The molalities and activities of all the aqueous species are reported under Distribution of species Here for each element the molality is given it is the same number as under Solution com position and the molalities of each of the aqueous species of that element sum up to this molality mass balance In our example it can be observed that 99 84 of all Ca is present as uncomplexed Ca The remainder is present as the complexes CaF and CaOH The last block of data consists of the Saturation indices calculated by the program 2 2 2 PHREEQC Exercise 1 Cape Karoo mineral water According to the label on the bottle Cape Karoo mineral water derives from dolerite rock The chemical composition of the water reflects the mineralogy of the host rock Using the concentrations reported on the bottle label we 19 Chapter 2 RAL Sp y A KAR sS ERAI we Figure 2 3 Label of Karoo mineral water can 1 check the accuracy of the chemical analysis and 2 calculate the speciation of this water The PHREEQC inp
7. CaCO3 is abundant in many geological settings and constitutes the most important source for calcium Ca in groundwater Ca also derives from other minerals such as dolomite gyp sum and feldspars and also takes part in cation exchange reactions Magne sium Mg finally has a high concentration in seawater and is contained in some minerals the most important being dolomite Similar to Na and Ca it can be adsorbed to the exchange complex Considering the numerous sources and processes that influence solute concentrations investigators that try and study the chemistry of aquifers are faced with an overwhelming complexity Traditional methods that are applied to make sense of the vast pile of numbers that follow from the chem ical analyses of water samples involve plotting and classification of sam ples into groups e g in order to be able to discern regional trends and to identify chemical processes No matter how useful these methods are in deriving an idea of the reaction scheme that has given the water sample its composition they are incapable of determining whether this scheme is feasi ble from a chemical point of view In order to check whether a concept obeys basic chemical theory one needs a geochemical model that is based on the firm laws of thermodynamics Or as put it Quantitative models force the investigator to validate or invalidate ideas by putting real numbers into an often vague hypothesis and thereby star
8. Chemical Nature of Reactions Water Resour Res 19 5 1231 1252 Schwarzenbach R P Gschwend P M Imboden D M 1994 Envi ronmental Organic Chemistry John Wiley amp Sons Inc New York Stumm W Morgan J J 1981 Aquatic Chemistry An Introduction Emphasizing Chemical Equilibria in Natural Waters John Wiley amp Sons Inc New York Stuyfzand P J 1993 Hydrochemistry and hydrology of the coastal dune area of the western Netherlands Vrije Universiteit Amsterdam The Netherlands PhD thesis Thierrin J Davis G Barber C 1995 A Ground Water Tracer Test 83 Chapter 4 with Deuterated Compounds for Monitoring In Situ Biodegradation and Retardation of Aromatic Hydrocarbons Ground Water 33 3 469 475 Thierrin J Davis G B Barber C Power T R Patterson B M Lambert M 1992 Groundwater tracer tests within a BTEX con taminated zone in order to determine aquifer parameters and natural degradation rates of organic compounds Tech rep CSIRO Division of Water Resources van der Lee J de Windt L 2001 Present state and future directions of modeling of geochemistry in hydrogeological systems Journal of contaminant hydrology 47 265 282 Walter A L Frind E O Blowes D W Ptacek C J 1994 Model ing of multicomponent reactive transport in groundwater 1 Model development and evaluation Water Resour Res 30 11 3137 3148 Zheng C Wang P P 1999 MT3DMS A mod
9. Next run another simulation with the same sorption parameters To add a decay reaction e Select Models gt MT3DMS Simulation Settings Change the Type of reaction to First order irreversible reaction e Then go again to Models gt MT3DMS gt Chemical Reactions gt View Reset Matrix Now you can edit the first order decay rate coeffi cients for both dissolved and sorbed phases to 0 00027397 1 day 0 1 1 year Make sure that the value is applied to all grid cells in the model domain e Export the simulated concentrations for the breakthrough curve to an external ASCII file You can then use Excel to compare your calculated breakthrough curves under different sorption and decay parameters 13 Chapter 2 Calculation of a sample s equilibrium composition Speciation and saturation indices This chapter will be dedicated to equilibrium chemistry The kinetics of geo chemical processes i e the description of how and at what rate a system approaches equilibrium will be discussed later The reader is also referred to the standard textbooks on aqueous geochemistry for an elaborate treat ment of this subject as well as for a more detailed treatment of equilibrium chemistry e g Essential to the description of the equilibrium composition of a geochem ical system is the law of mass action According to this law the distribution of the species in the reaction aA bB amp cC dD 2 1 is given at equilibrium
10. Umar iS an as symptotic maximum specific uptake rate and Korg is the half saturation constant which is the substrate concentration at which the actual up take rate equals Umaz 2 Based on equation 4 1 the total uptake rate Um considers the dependency of the change of microbial mass on the actual microbial concentration X itself Cora Kora F Coty An additional potential growth limitation by electron acceptor availability might be incorporated into 4 2 leading to Cori Cea Um Umax 7 gt Korg TF Corg Kei A Cen where Cea is the electron acceptor concentration and Kea is the appro priate half saturation constant The complete mass balance equation for the microbial mass X describing the change of microbial concentration as a function of time includes both a microbial growth and decay term Ox OX growth i OX decay ot ot ot 4 1 X 4 2 Um Umax 4 3 4 4 52 Chemical kinetics with Ox growth Cor Cea Y max X YUm 4 5 Ot 7 Korg F Corg Kea T Cea l in which Y is a stoichiometric factor and OX decy N 4 6 Ot Vd 4 6 where Ugec is a decay rate constant During growth Um gt 0 both organic substrate and electron ac ceptors are consumed at rates that are proportional to v which will be illustrated in the following example of toluene degradation under sulfate reducing conditions 4 1 1 PHREEQC Example Kinetic toluene degradation The stoichiometry of the reactions that des
11. get 3 78 Chemical kinetics 20 ratio 1 get 4 30 moles parm 1 1000 1 ratio rate time 40 save moles end For each of the chlorinated ethenes the carbon isotope ratios 13C 1 C can be computed from the simulated concentrations of the light and heavy isotopes The calculation is done analogous to the calculation of measured concentration values The carbon isotope ratio is expressed as permil deviations from the Vienna Pee Dee Belemnite V PDB standard in the conventional no tation For the known simulated TC E and TCE concentrations the carbon isotope ratio 6 C rcp can be calculated from 8TCE Y TCE 8C O sta SBTC E C YC sua x 1000 4 19 with C C sta 0 0112372 4 20 Using these equations we can now determine the simulated 6 C values for all chlorinated ethenes e Open the spreadsheet pht3d_exercise_6 xls e In iPHT3D extract the select the various species to extract sim ulated concentrations at specific locations in the model e In iPHT3D Extract the simulated concentration values for Pecel and Pce_h at a location within the Western contaminant source x 850 z 6 5 and near the pumping well e Paste the values into the spreadsheet to calculate the 6 TCE for both locations e Extract the simulated concentration values for Tce_l and Tce_h at a location within the Eastern contaminant source x 70 m z 10 m e Calculate the 6 TCE for both locations e What
12. 55 y Position recovery well m 5 1 2 Flow rate injection well m3 day 500 1 2 Flow rate recovery well m3 day 500 Longitudinal dispersivity m 1 Transversal dispersivity m 0 1 that the geometry is defined for a half model due to the symmetry of the problem Neglect background groundwater flow and use for simplicity fixed head cells at the left x 0 m and the right x 200 m model boundary to allow exchange of water across these boundaries e Check your simulated heads in the injection well 24 22 m and in the recovery well 15 78 m 4 3 2 Setting up the nonreactive transport problem e If the simulated heads agree use MT3DMS to simulate a sim ple tracer experiment to estimate the travel times between the injection and the recovery well Set the tracer concentration in the backgroundwater to 0 and set the tracer concentration in the 70 Chemical kinetics injection well to 1 Include some observation points on the axis between injection and recovery well to allow the visualisation of breakthrough curves e What is your simulated travel time to the mid point between in jection and recovery well It should be approximately 20 days If your result differs try to debug the problem before proceeding e When does the recovery well receive 100 of the injection con centration 4 3 3 Setting up the reactive transport problem Once the tests with the nonreactive model indicate that the model is working properly pro
13. CSIRO Division of Water Resources Bethke C M 1996 Geochemical reaction modelling Oxford Univer sity Press Inc USA Broholm M M Jones I Torstensson D Arvin E 1998 Ground water contamination form a coal carbonization plant In Lerner D N Walton N R G Eds Contaminated land and groundwa ter Future directions Engineering Geology Special Publication The Geological Society London pp 159 165 Custodio E 1987 Hydrogeochemistry and tracers Studies and Re ports in Hydrology pp 213 269 Davis G B Barber C Power T R Thierrin J Patterson B M Rayner J L Wu Q 1999 The variability and intrinsic remedi ation of a BTEX plume in anaerobic sulphate rich groundwater J Contam Hydrol 36 3 4 265 290 Davison R M 1998 Natural attenuation and risk assesment of groundwater contaminated with ammonium and phenolics Univer sity of Bradford Bradford phD thesis Davison R M Lerner D N 2000 Proving natural attenuation with modeling at a site contaminated by ammonium and phenolics In Sililo O Ed Groundwater Past achievements and future chal lenges A A Balkema pp 737 740 Eberhardt C Grathwohl P 2002 Time scales of organic contaminant dissolution from complex source zones coal tar pools vs blobs J Contam Hydrol 59 1 2 45 66 Engesgaard P Kipp K L 1992 A Geochemical Transport Model for Redox Controlled Movement of Mineral Fronts in Gr
14. Natural Attenuation of remnant coking liquor contamination from a deep unsaturated zone In Groundwater Quality 2001 Conference pre prints Sheffield UK pp 331 334 Kindred J S Celia M A 1989 Contaminant transport and biodegradation 2 Conceptual model and test simulation Water Re sour Res 25 6 1149 1160 Lichtner P C 1996 Continuum formulation of multicomponent mul tiphase reactive transport In Lichtner P C Steefel C I Oelkers E H Eds Reactive Transport in Porous Media General Principles and Applications to Geochemical Processes Mineralogical Society of America Washington D C pp 1 79 Michaelis L Menten M L 1913 Die Kinetik der Invertinwerkung Biochemische Zeitschrift 49 333 369 Miller C T Poirier McNeill M M Mayer A S 1990 Dissolution of trapped nonaqueous phase liquids Mass transfer characteristics Water Resources Research 26 11 2783 2796 Monod J 1949 The growth of bacterial cultures Ann Rev Microb 3 371 394 82 Chemical kinetics Noorishad J Carnahan C L Benson L V 1987 Develop ment of the Non Equilibrium Reactive Chemical Transport Code CHMTRNS Tech rep Lawrence Berkeley Laboratory Earth Sci ence Division Parkhurst D L Appelo C 1999 User s guide to PHREEQC A com puter program for speciation reaction path 1D transport and in verse geochemical calculations Tech rep U S Geological Survey Water Resources
15. No additional input is needed here e In the Results Panel select the final time step of the simulation by going to Results Aquifer gt T step and selecting 0 24305 days e In the Results sub panel select one of the specific species that you want to visualise by going to Results Chemistry gt Species and selecting for example Ca e Then go to Results Observation and select the visualisation type Profile Then select the Profile option e in the following dialog you can either select a single or multiple species For pure visualisation the selection of multiple species makes only sense for all those species which prevail within the same concentration range However in this exercise we want to export the data for multiple species and therefore we select all of the simulated species e Once the species are selected a new dialog will open This new dia logue allows to make a selection among several types of concentration averaging techniques In the present case simply select the default option i e Value e Now a plot should appear that shows the simulated concentrations of all selected species along the previously defined linel i e concentra tion profiles along the column should show up e You can now use the Export button to export the simulated values into a text file 3 1 8 Comparison of simulation results If all steps of the model definition were performed correctly the results should agree closely with those ota
16. Simulation Settings Run the model Visualise the microbial fringe 67 Chapter 4 Table 4 7 Measured and equilibrated aqueous concentrations of the ambient water composition at the Eden Hill field site Aqueous component Cinit Crech mol 1 pH 5 14 pe 8 01 O 0 0 Neot 1 00 x 107 4 C 4 7 96 x 1073 C 4 0 S 6 7 85 x 10 4 S 2 0 Fe 2 5 20 x 1075 Fe 3 1 46 x 1075 Ca 8 42 x 1074 K 2 59 x 1074 Mg 5 67 x 1074 Na 5 25 x 1078 Cl 6 46 x 1073 Table 4 8 Rate constants for the reactive transport simulation of BTEX degradation under sulfate reducing conditions Rate Expression Parameter Nr Parameter Unit Value Benz il Umaz sul f benzene 1 day 0 0 Tolu 1 Umaz sul f toluene 1 day 5 0 Ethy il Umaz sulf ethylbenzene 1 day 2 0 Xyl 1 Uman aulf aylene 1 day 2 0 Benz 2 Obiogman mol 1 1 0 107 Tolu 2 Oiio max mol 1 1 0 1075 Ethy 2 Obiommas mol l 1 0 107 Xyl 2 Onecare mol l 1 0 1075 Benznapl 1 not used 1 Tolunapl 1 not used 1 Ethynapl 1 not used 1 Xylnapl 1 not used 1 Benznapl 2 w 1 day 0 025 Tolunapl 2 w 1 day 0 025 Ethynapl 2 w 1 day 0 025 Xylnapl 2 w 1 day 0 025 Sulfred 1 Udec 1 day 0 1 Sulfred 2 Cinit mol 1 1 0 1077 68 Chemical kinetics 4 3 PHT3D Exercise 5 Pyrite oxidation during deep well injection Managed aquifer recharge is increasingly used to enhance the sustain able development of water supplies Common recharge techniques in clude aquifer storage and
17. and available as btex_degradation phrq However you still need to insert some numbers to make it work Fill in the parameters from table 4 4 at the appropriate locations these are marked with signs Run the file and inspect the change of BTEX with time on the Chart tab Write down the order in which the 4 different compounds degrade 1 im fer The USER_GRAPH keyword block is set up in such a way that the concentrations of the individual BTEX compounds are plotted in the chart Modify it to have the concentrations of the biomass and nitrogen displayed instead e Explain the temporal concentration patterns of the biomass and nitrogen 58 Chemical kinetics 4 2 PHT3D Exercise 4 Microbially mediated petroleum hydrocarbon degradation under sulfate reducing conditions Numerical modelling is widely used to analyse and predict the risk associated with using natural attenuation as a remediation scheme With few exceptions most practical model applications to contami nated field sites involve a vertically averaged two dimensional simula tion model This assumes that vertical concentration gradients over the aquifer depth are negligible However as pointed out for example by transversal vertical mixing can be a critical physical factor for controlling the length of naturally attenuating contaminant plumes Unfortunately most contaminated sites are not well characterised with respect to the vertical distribution of dissolve
18. and kinetic pro cesses which will be the topic of the present chapter The PHREEQC cases that will be discussed were taken from 51 Chapter 4 4 1 Modelling bioremediation processes In order to be able to design active bioremediation systems and to understand passive bioremediation natural attenuation mechanistic descriptions that quantify microbial activity are needed Since both the rate of microbial growth and the rate of contaminant utilization are highly dependent on the amount of biomass available to catalyze the re actions such models must have the capability to predict both transient and spatial variations in biomass To model the production of biomass and the related consumption and production of other chemicals the key steps are 1 to formulate the rate expressions for the reaction ki netics and 2 to determine the stoichiometry of the biodegradation reactions and in order to describe the temporal variations In the macroscopic mathematical descriptions of microbial growth dynamics aimed at governing laboratory or field scale processes many of the complex interdependencies that are known at the microscopic scale are commonly neglected and described by empirical formulations based on the classical works of and The expression describing a specific bacterial growth rate Vsp observed in many batch experiments is v v Corg sp marty n Korg T Gorg where Corg is the concentration of the organic substrate
19. can be concluded from these values with respect to the most likely origin of the vinyl chloride VC Has it originated from the site of company W or the site of company E 79 Chapter 4 Table 4 10 Compilation of selected model input parameter Total simulation time days Time step length reactive transport days Groundwater recharge rate m day Prescribed heads Western boundary m Eastern boundary m River m Pumping rate m3 day Hydraulic conductivity coarse grained material m day Hydraulic conductivity clayey material m day Longitudinal dispersivity m Transversal dispersivity m Source 1 CPcenapl mol l Source 2 CTcenapl mol l CTolunapl mol l Ambient recharge and constant head inflow concentrations Coxygen mol l pH pe CChlorinated Ethenes mol l CToluene mol l 80 3650 36 5 0 001 13 7 14 5 10 8 7 5 40 0 1 0 5 0 005 4 x 10 4 13 5 Bibliography Appelo C A J Postma D 1993 Geochemistry Groundwater and Pollution A A Balkema Rotterdam Appelo C A J Postma D 2005 Geochemistry Groundwater and Pollution second edition A A Balkema Rotterdam Barber C Davis G B Thierrin J Bates L Patterson B M Pribac F Gibbs R Power T R Briegel D Lambert M Hosk ing J 1991 Final report for project on Assessment of the Impact of Pollutants on Groundwater beneath Urban Areas July 1989 to June 1991 Tech rep
20. cases the lines need to extend over the whole depth of the aquifer The value for both lines zones needs to be set to 1 Use for example bc_up and bc_down as Zone Names You will see the selected names displayed near the location of the lines This will define a the type of boundary condition to a fixed head boundary at the inflow and outflow model boundaries The hydraulic head at these grid cells will therefore remain unchanged over the whole simulation period It will remain at the same values that we will now define as initial heads Set initial hydraulic heads by going to Spatial Attributes MOD FLOW BAS6 bas 6 Initial Heads and set one_value to 10 cor responding to a head of 10m Then press ok to confirm the value In addition define a zone line at the upstream end of the model domain where you also defined that the type of boundary condition is 1 and allocate a hydraulic head of 12m over the entire depth of the aquifer Mineral dissolution and precipitation 3 2 2 Running MODFLOW You have now completed to set up the flow model and you are ready to run MODFLOW e Under Parameters Flow use the Pencil Button to write all the input files required to run MODFLOW e Under Parameters Flow use the Red Arrow button to start the MODFLOW simulation e Check the simulated head contours and assess whether the results are plausible 3 2 3 Data input for solute transport properties For the simulation of the
21. data to an ASCII text file for use with another program such as Excel To view a contour plot of the simulated concentrations after a simulation time of 3 years e Select Tools gt 2D Visualisation On the MT3DMS tab select Tracer and click OK Click on the right arrow next to Simulation Time until you have reached a value of 1095 days As you move forward in time 10 Single Species Transport E MT3DMS_EX1 PM5 Processing Modflow Pro File Value Options Help Layer gt Row Column Simulation Time t ge Y a a le i f rT GES z T 584i639 155 5 f 1 1 12 2D Visualization Solute Concentration MT3DMS Species 1 5 864385E 25 Figure 1 4 Calculated concentration distribution in layer 1 at year 3 under natural gradient conditions you should see the contaminant plume growing see Figure 1 4 for comparison e Select Options Environment and select the Contours Tab There you can for example define color ranges Click on the Fill Tab acti vate Fill contours and define the contour levels that will be displayed e Note that you can easily switch over to a cross sectional view 1 1 8 Effect of Solution Techniques Now let us examine the effect of transport solution techniques on the accu racy of simulation results Export the breakthrough curve at the observation point obtained using the TVD option to an ASCII file Then change the advection solution option to Upstream Fi
22. day With the selected model dimensions a steady state flow rate Qyey of 0 259 m d7 is required to achieve the pore velocity of 0 81 m d7 as defined by for a porosity of 0 32 A summary of the parameters that define the flow field and the non reactive transport simulation is given in Table 3 1 To implement this model setup into iPHT3D proceed as follows e In iPHTS3D create a new model and save it in a new and separate directory to avoid any mixup of different models and model results 27 Chapter 3 28 e Under Parameters Model select the Grid button and specify the model domain in the dialog that opens up Use the values as discussed above and listed in Table 3 1 Under Parameters select the Time button to define the simulation time of 0 2435 days under Total Simulation Time and the Step Size to a value of 0 01 day Go to Spatial Attributes gt MODFLOW DIS gt dis 6 Top of Layers TOP and set the top of the model to 1 m Don t forget to validate the input by clicking the ok button Go to Spatial Attributes gt MODFLOW DIS gt dis 7 Bottom of Layers BOT and set the value to 0 m the default value Set the boundary conditions for heads by going to Spatial Attributes MODFLOW gt BAS6 gt bas 3 Boundary Conditions and setting the value to 1 active cell for all cells After that select Zone and use the Zone Tools to add a line at the last cell which represents the downstream
23. initial and inflow concentrations you will enter the chemical composition of the respective aqueous solutions In this exercise the background contains water that is in chemical equilibrium with the mineral Calcite The composition of the initial water needs to be entered under the 31 Chapter 3 Solutions Tab in the Background column On the other hand the initial mineral composition must be entered under the Phase Tab with the value of the saturation index SI and the initial concentration in the assemblage Leave the SI at the default value of 0 as would be the case for most typical model applications The inflow water composition will be defined in the Solution 1 column under the Solutions Tab To additionally define that the inflow composition is used at the correct location i e the inflow boundary we need to select Spatial Attributes gt PHT3D PH gt ph 4 Source Sink and define a zone point near the inflow boundary and allocate a Zone Value of 1 This induces that the previously defined concentrations for Solution 1 as defined in the Solutions Tab will be used to specify the concentrations in the water that enters the column 3 1 5 Data input for solute transport For the simulation of solute transport of the mobile components you need to specify the model parameters that control advective and dispersive trans port It is assumed that the parameters that define the physical transport of the chemicals are the same for eve
24. input for solute transport properties For the simulation of the solute transport of the mobile components the model parameters that control advective and dispersive transport need to be entered e Go to Parameters Transport and select the Parameter button e In the dialog box that subsequently appears select ADV and se lect under adv 1 Solver Parameters the 3rd order TVD Scheme ULTIMATE as the Solution Scheme Leave the default param eters and click OK 43 Chapter 3 Select Spatial Attributes gt MT3DMS DSP enter a longitu dinal dispersivity value of 0 5 m and click OK Select Spatial Attributes MT3DMS DSP and specify that the ratio between transverse vertical dispersivity and the longitu dinal dispersivity is 0 01 This will define that the actual value is 0 005 m 5 mm Make a quick test and assess the conservative transport behaviour in this setting To do this select Spatial Attributes MT3DMS BTN gt btn 13 Concentrations and specify a small zone in the upper left corner of the model near the inflow boundary for which a concetration of 1 or any other concentration gt 0 is defined Under Parameters Transport use the Pencil Button to write all the input files required to run MT3DMS Under Parameters Transport use the Red Arrow button to start the MT3DMS simulation Inspect the results by navigating in the Results panel of ipht3d For example to see the concentration contours
25. recovery ASR infiltration ponds river bank filtration and deep well injection Following recharge the water qual ity of the injectant is typically altered by a multitude of geochemical processes during subsurface passage and storage Relevant geochem ical processes that affect the major ion chemistry include microbially mediated redox reactions mineral dissolution precipitation sorption and ion exchange The hydrochemical conditions and changes that oc cur under these circumstances in particular the temporal and spatial changes of pH and redox conditions are in many cases the control ling factor for the fate of micropollutants such as herbicides and phar maceuticals Similarly changes in mineralogical composition such as dissolution and precipitation of iron or aluminiumoxides may affect the mobility of trace metals as well as the attachment and subsequent decay of pathogenic viruses Laboratory and field scale experimental studies are aimed at investigating such processes under controlled con ditions and to eventually develop a better qualitative and quantitative understanding of their complex interactions both site specific and at a fundamental level carried out a reactive transport modelling study to analyse the data collected during a deep well injection experiment in an anaerobic pyritic aquifer near Someren in Southern Netherlands This exercise replicates some of the key processes that were iden tified to influence water qua
26. select Results gt and select Tracer Then select the first timestep under gt T step and observe the plume spreading by using the right arrow button to move to subsquent timesteps 3 3 3 Reactive transport Based on the flow and conservative transport model the next step is to define the input for the reactive transport 44 e Copy the database file that will used for this PHT3D exercise into the folder that contains all files for this exercise and make sure its name is pht3d_datab dat You can use the same dataase as in PHT3D Exercise 1 i e the standard PHREEQC database Go to Parameters gt Chemistry and select the Import database button iPHT3D will then read and interpret the PHT3D database file In this next step we define all initial concentrations of all aqueous components and also the two types of water compositions at the upstream model boundary which will be allocated to the inflow ing groundwater Activate the equilibrium components tick the checkboxes that are included and also enter the concentrations Mineral dissolution and precipitation provided in 3 8 Note that there is no need to include aqueous components that will not play a role in the simulations such as Mn 2 Mn 3 and Si However consider and activate all va lence states of redox sensitive components For example activate all valence states of nitrogen i e N 5 N 3 N 0 and Amm The composition of the ambient water initial wate
27. shield and reduce the reactivity of the ion This effect can be corrected for by using a so called activity coefficient that relates the activity of an ion to its concentration for example Ca MA 24 Ca t yoa2 ae Gh ites ai 2 6 Ca Tt where y is the activity coefficient which is multiplied with the concen tration of the ion mo 2 in this example divided by the standard state that for practical purposes conveniently equals 1 mole kg H20 According to the Debije Hiickel theory activity coefficients are a function of the ionic strength I of the solution r 1 2S gt m 2 2 7 Several empirical relationships exist to calculate the activity coefficients from the ionic strength e g the Davies equation at 25 C T log yi 0 5085 z e a oat 2 8 where z is the charge of ion i e g 2 for Ca Although these calculations can be awkward they are straightforward and can still be done by hand so the effect of electrostatic shielding is quite easily corrected for 2 2 Formation of aqueous complexes The second reason why the activity of an ion is lower than its concentration is that it forms aqueous complexes with other ions Laboratory analyses 16 Specation and saturation indices Chemical analysis xCa USO 1 lonic Strength I eq 1 7 2 Activities of ions Ca SO eq 1 6 1 8 3 Activities of complexes CaSO eq 1 2 4 Concentrations of complexes Measos
28. solute transport of the mobile components the model parameters that control advective and dispersive transport need to be entered e Go to Parameters Transport and select the Parameter button e In the dialog box that subsequently appears select ADV and select un der adv 1 Solver Parameters the 3rd order TVD Scheme ULTIMATE as the Solution Scheme Leave the default parameters and click ok e Select Spatial Attributes gt MT3DMS DSP enter a longitudinal dispersivity value of 2 5 m and click OK e Select Spatial Attributes MT3DMS DSP and specify that the ratio between transverse vertical dispersivity and the longitudinal dis persivity is 0 01 This will define that the actual value is 0 025 m 25mm 3 2 4 Data input for reactive transport In this step we define the reaction network and the various water and mineral compositions that play a role in this problem Both the initial and the inflow concentrations of of the aqueous components that will be used in the simulation are shown in Table 3 5 Also the initial concentrations of the minerals are given in Table 3 6 e Copy the database file that is going to be used for the PHT3D simula tion into the folder that contains all files for this exercise and name the file pht3d_datab dat You can use the same dataase as in PHT3D Exercise 1 i e the standard PHREEQC database 37 Chapter 3 38 e Go to Parameters gt Chemistry and select the Import database b
29. the pumping well 4 4 66 Carbon isotopes 0 0 20000 2 ee Acknowledgments Appendix I Appendix II 81 83 85 Introduction Knowledge of the chemistry of groundwater is a requirement for a number of practical purposes As groundwater is an important source for drinking water one has to ascertain that its quality is sufficient for consumption Quality requirements are equally important for other types of utilization such as irrigation or industrial purposes as well as for the protection of vulnerable ecosystems More recently pollution and clean up of aquifers has become a major topic in aqueous geochemistry Furthermore understanding of geochemical processes is needed for safety assessment studies e g for the storage of nuclear waste Clearly there are numerous practical applications for aqueous geochem istry Moreover geochemistry is an essential tool for understanding the hydrogeological systems that we study It provides information on the prove nance of groundwater on flow directions and on groundwater ages Water quality patterns in aquifers can be complex The input of dif ferent sources of water is the first of factors that adds to this complexity Sources include precipitation rivers lakes possibly polluted or saline due to strong evaporation seawater ascending deep groundwater and anthro pogenic sources such as wastewater or irrigation return flow Geochemical processes add to the complexity
30. 8 15 a graphical user interface for MODFLOW MT3DMS SEAWAT and PHT3D For information on how to order a fully func tional copy of the software contact the author of the program Wen Hsing Chiang wchiang pmwin net or see http www simcore com Course outline Many students and researchers that are introduced to reactive transport modelling for the first time are overwhelmed by the apparent complexity of it Mastering the subject requires not only understanding of the underly ing theory of groundwater flow solute transport and hydrochemistry but also of the many options that the modelling software packages offer It of ten takes quite some time before students are confident enough to build and apply their own models Although reactive transport modelling can indeed be complex one should not be scared away from it as mastering this subject provides powerful applications Even the simplest calculations of equilibrium speciation and saturation states already add to the knowledge of and insight in the chemical system that is being studied Each student in hydrochemistry should therefore have at least basic knowledge of aqueous geochemical models This course guide is a step by step introduction of the topic of reac tive transport modelling It is assumed that the reader has some knowledge of hydrochemistry which is indispensible for successful application of the models The first chapter introduces the basic theory of geochemical models and
31. Cells Click Edit to start editing the concentration for the species Tracer or any other name that you have chosen for the species Add a pollution source by setting the Flag and the concentration value at the grid cell positioned at Column1 Row 16 Layer 1 to a value of 1 for the first stress period Set the value to 0 for the second stress period Leave the editor and save your changes Define the initial concentrations by selecting Models MT3DMS gt Initial Concentrations Click Edit to start editing the initial concentra tion of Tracer or any other name that you have chosen for the species Use Reset Matrix to change the value for all initial concentrations to 0 Define the advection package by selecting Models gt MT3DMS gt Ad vection Select 3rd order TVD Scheme ULTIMATE as the Solution Scheme and click OK Define the dispersivity to be used by selecting Models MT3DMS Dispersion In the window that appears leave TRPT at a value of 0 1 Change TRPV to 0 01 for all three layers and then click OK In the editor use Reset Matrix to change the value for the longitudinal dispersivity to 1 0 m Make sure that you apply these values to all three layers Leave the editor and save your changes To visualise not only the results from the last time steps of each stress period but also some intermediate results select Models MT3DMS Output Control Click on the Output Times tab and configure the settings such th
32. Chemical equilibrium of water compositions Are aqueous so lutions that are provided as initial water composition s pre equilibrated with the equilibrium minerals that are included in the simulation e Chemical equilibrium of water compositions Are aqueous so lutions that are provided as initial water composition s pre equilibrated with the ion exchanger sites and have the concen trations on the exchanger been entered 89 Chapter 4 e Redox Have all redox states of redox sensitive components been included in the definition of the reaction module For example for sulphur not only S 6 but both S 6 and S 2 and or both C 4 and C 4 etc 90
33. EEQC Exercise 2 Fluoride rich waters Fluorite is a common mineral in volcanic rocks and an important source of fluoride in groundwater High fluoride concentration in drinking water are known to cause fluorosis a painful crippling disease That is why the WHO drinking water standard is 1 5 mg l although this may even be too high for arid areas where people consume a lot of water Table 2 1 Examples of fluoride rich waters Concentrations are in mmol l except pH A Maarum Denmark B Rajasthan India C Lake Abiata Kenya pH 7 8 7 3 9 62 Nat 19 1 47 9 194 Kt 0 36 0 15 4 91 Meg 1 19 0 79 0 02 Ca2t 1 05 0 68 0 042 i 5 67 17 4 53 9 HCO 17 8 14 8 138 s02 0 0 5 2 0 15 NO 0 03 7 8 s F 0 089 0 356 6 28 23 Chapter 2 Ca moll 24 0 009 0 008 0 007 0 006 0 005 0 004 0 003 0 002 0 001 0 0 0001 0 0002 0 0003 0 0004 0 0005 0 0006 0 0007 0 0008 0 0009 0 001 Fluoride mol l Figure 2 5 Plot of Calcium vs Fluoride concentrations Calculate the saturation index SI of the water samples in table 2 1 for fluorite CaF2 Ca t 2F7 K fluorite 107106 without making corrections for the ionic strength and complexes Results ST fluorite A SF teens ST ftuorite B E ORS SI fluorite Fee ee Calculate STfiuorite With PHREEQC The input file is available as fluorite1 phrq The concentrations of the elements are entered under SOLUTION Remember that NO is entered as N 5 and S
34. Investigations Report in review Plummer L N 1992 Geochemical modeling of water rock interaction Past present future A A Balkema pp 23 33 Prommer H Barry D A 2005 Modeling Bioremediation of Contami nated Groundwater American Society for Microbiology pp 108 138 Prommer H Barry D A Davis G B 1998 The effect of seasonal variability on intrinsic biodegradation of a BTEX plume In Her bert M Kovar K Eds Groundwater Quality Remediation and Protection pp 213 220 Proc GQ 98 Conf Tiibingen Germany 21 25 September 1998 IAHS Public 250 Prommer H Barry D A Davis G B 2002 Influence of transient groundwater flow on physical and reactive processes during biodegra dation of a hydrocarbon plume J Contam Hydrol 59 113 131 Prommer H Barry D A Zheng C 2003 MODFLOW MT3DMS based reactive multi component transport modelling Ground Water 42 2 247 257 Prommer H Davis G B Barry D A 1999 Geochemical changes during biodegradation of petroleum hydrocarbons Field investiga tion and modelling Org Geochem 30 6 423 435 Prommer H Stuyfzand P J 2005 Identification of temperature dependent water quality changes during a deep well injection experi ment in a pyritic aquifer Environmental Science and Technology 39 2200 2209 Rubin J 1983 Transport of Reacting Solutes in Porous Media Re lation Between Mathematical Nature of Problem Formulation and
35. O as S 6 Browse through the output file to find the saturation indices of each of the water samples ST fluorite A bes Ses ST fiuorite B Fees ST fluorite Fee Why is there a difference between hand and computer calculations Use figure 2 5 to make a plot of the concentration of Ca against the concentration of F Calculate the Ca concentration of water in equilibrium with fluorite at a F concentration of Use the values from the table above to plot the equilibrium line for fluorite in the graph 1 2 3 Specation and saturation indices Sample F Ca mol l 1 5 264 1077 mol l 1 0 mg l 2 7 895 1075 mol l 1 5 mg l 3 1 579 1074 mol l 3 0 mg l Samples that plot on the solubility curve are in equilibrium with flu orite Samples that plot above or below the curve are super and subsaturated respectively Will a 1 1 mixture of samples 1 and 3 be in equilibrium super or subsaturated with respect to fluorite In the hand calculations the difference between concentration and ac tivity is not taken into account A more accurate calculation is done with PHREEQC The basic set up of the input file is given in fluo rite2 phrq Open it in the editor fill in the missing numbers and run it Sample F Ca t mol l ratio hand computer 5 264 107 mol l 1 0 mg l m a s see 7 895 1075 mol l 1 5 mg l moane 0 aa 1 579 10 mol l 3 0 mg l eee Why does the ra
36. Period length and a second stress period of 2555 days e Set the boundary conditions for heads by selecting Select Grid gt Cell Status IBOUND Modflow Make sure that the value of IBOUND is or becomes 1 active cell at all cells except the cells at the upstream column 1 and the downstream boundary column 50 Change the value of IBOUND at the upstream and downstream boundary cells to 1 which will define those cells as fixed head cells The hydraulic heads at these grid cells will remain at the allocated values during the entire simulation To enter the data edit the values in Layer 1 Afterwards copy the values to Layers 2 and 3 Leave the editor and save your changes Single Species Transport e Select Parameters gt Initial amp Prescribed Hydraulic Heads Reset Matrix and set the value to 10 m and click ok Then change the value at all cells at the upstream boundary column 1 to 18 33 m Make sure that you apply these values to all three layers Leave the editor and save your changes e Select Parameters Horizontal Hydraulic Conductivity Reset Ma trix and set the value to 0 8640 m day Make sure that you apply these values to all three layers Leave the editor and save your changes e Select Parameters gt Vertical Hydraulic Conductivity Reset Matrix and set the value to 0 00864 m day Make sure that you apply these values to all three layers Leave the editor and save your changes e Select Parameter
37. Short Course Applied Reactive Transport Modelling University of Bordeaux France 24 28 March 2014 Prepared by Henning Prommer Henning Prommer csiro au Vincent Post vincent post flinders edu au Chunmiao Zheng czhengQua edu James A Davis mrwaterx gmail com Olivier Atteia olivier atteia ipb fr Contents Introduction 1 1 Single Species Transport 5 1 1 MT8DMS Exercise 1 Part 1 Single Species 3 D Transport Simulations ae 2 44 Aeesha Sd eo edge eee ed oe eee 5 1 1 1 Groundwater Flow 0 5 1 1 2 Solute Transport 2 0 5 1 1 3 Setting up the flow model 6 1 1 4 Running MODFLOW 8 1 1 5 Visualising hydraulic heads 8 1 1 6 Setting up the transport model 9 1 1 7 Breakthrough curves and contour plots 10 1 1 8 Effect of Solution Techniques 1i 1 1 9 Optional Numerical Experiment 12 1 2 MT3DMS Exercise 1 Part 2 Single Species 3 D Transport with with Sorption and First Order Decay 12 1 2 1 Linear equilibrium controlled sorption 12 1 2 2 First order Radioactive Decay or Biodegradation 13 Calculation of a sample s equilibrium composition Specia tion and saturation indices 15 2 1 Electrostatic shielding 2 16 2 2 Formation of aqueous complexes 16 2 2 1 Your first PHREEQC input file 18 2 2 2 PHREEQC Exercise 1 Ca
38. _l Dce_l Vc_l and Eth_l e From those concentration alone is it possible to determine whether the contamination in the well originated from Company E or from Company W e Check the oxygen concentrations at the pumping well use again the grid cell defined by Column 13 Layer 20 What would this suggest in terms of potential for reductive dechlorination 4 4 6 Carbon isotopes To analyse this modelling scenario further we make now use of the sim ulated carbon isotope ratios In the model each chlorinated ethene species was divided into two separate species that represent the light and heavy isotopes respectively That means that each dissolved com pound occurs twice in the reaction database to represent the molecules containing either C or C Degradation of the light isotopes typically proceeds faster than the reaction for the heavier isotopes The difference between the reaction rates depends on the enrich ment factor e The enrichment factor for a specific reaction may be taken from the literature or where possible determined by site specific experimental work The slight difference in reaction rates and the resulting C in the mother product is incorporated into the reaction rate expression for the heavier isotope For example in the case of Tce_h the rate expression is Tce_h i start 2 if tot O 0 gt le 7 then goto 40 5 if tot Tce_l tot Tce_h lt 1e 9 then goto 40 10 rate
39. ake place The transformation follows typically a sequential degradation pathway TCE DCE gt VC gt Ethene In our model this reaction pathway is incorporated as a sequence of simple first oder reactions The corresponding reaction rate expression is for exmaple for TCE Tce_l start 2 if tot O 0 gt le 7 then goto 60 5 if tot Tce_l tot Tce_h lt le 9 then goto 60 10 rate parm 1 tot Tce_l tot Tce_h 20 ratio tot Tce1 tot Tce_l tot Tce_h 30 moles ratio rate time 40 put rate 3 50 put ratio 4 60 save moles end Note that the incorporation of an if statement line 2 inhibits the reaction in the presence of oxygen The DCE produced in this reac tion is subsequently further degraded to VC the most toxic compound among the chlorinated ethenes While PCE could be transformed to TCE under reductive conditions it does not degrade under the aerobic conditions that prevail on the Western side of the mode domain 77 Chapter 4 4 4 5 Chlorinated ethenes at the pumping well As aresult of the transport and chemical reactions we can find a com plex mix of chlorinated ethenes in the pumping well e Open the model and inspect the simulation results at the pumping well Use a particular grid cell that is representing the pumping well Column 13 Layer 20 for comparison of the results e What are the concentrations of Pce_l Tce
40. at explains the different characteristics of the BTCs In the final phase of this exercise we replace the Na Cl solution with the measured and charge balanced water composition from the Table 4 6 Single species solubilities of BTEX compounds used for the sim ulations in PHT3D Exercise 4 Aqueous component ce mol l Benzene 0 022820 Toluene 0 005978 Ethylbenzene 0 00193 Xylene 0 001800 65 Chapter 4 Eden Hill site as given in and as listed in Table 4 7 The water is anaerobic and only sulphate will act as electron acceptor As discussed in the earlier PHREEQC example exercise we will include microbial growth and decay in the simulations Microbial growth will only occur where both toluene and sulfate are present simultaneously The rate expression used in the PHREEQC Exercise 5 was slightly modified and compared to Eqn 4 3 an additional biomass inhibition term was included Corg Cea Korg JE Corg Kea Cea This additional term bio was introduced to reflect a conceptual model or biologically mediated degradation reactions which suggest that with increasing biomass concentrations and thus increasing biofilm thickness degradation rates might become limited by the supply of re actants In order to avoid solving the diffusional transport of reactants at a microscopic i e within biofilm level a macroscopic formulation can be used to account for the rate limitation resulting from excessive biomass accumulation The biomass i
41. at results are printed every 50 days by setting the value for Interval to 50 Run the transport model by selecting Models MT3DMS gt Run and clicking OK Chapter 1 Time Series Curves Concentration Data Chart Axis Time Concentration Time Curve J FixBounds Logarithmic Lower Bound jo Upper Bound 3650 Reset Bounds Y Axis I FixBounds Logarithmic Lower Bound 0 Upper Bound 4973017 Reset Bounds Data Type M Calculated M Observed Use results of all observations gt w E 2 ki 6 ts o fa E S oO gt N 2000 Use results of the following Time BSNAM z Copy to Clipboard Save Plot As ok Cancel Help Figure 1 3 Concentration breakthrough curve at observation point 1 col umn 16 row 16 and layer 1 1 1 7 Breakthrough curves and contour plots To visualise the breakthrough curve for the contaminant at a specific location proceed as follows e Define the observation well by selecting Models gt MT3DMS gt Con centration Observations For the Observation Borehole add 77 5 as a value for X easting and similarely use 77 5 as value for Y northing Under Observation Data set the Proportion to 1 for Layer 1 Click OK e Visualise the breakthrough curve by selecting Models MT3DMS gt View Concentration Time Curve Select Tracer and click OK Click on the Chart Tab You can save the result as plot or you can go to the Data Tab export the
42. behind as a result of the cation exchange of ammonium In the modelling example the development of an ammonium plume and the subsequent flushing of ammonium contaminated groundwater by pristine background water is simulated The processes included in the example are advection dispersion cation exchange and the kinetically controlled oxidation of ammonium Dis persion ion exchange and nitrification act as attenuation processes for ammonium For simplicity the two dimensional reactive transport problem is set up for the same model domain as the previous exercise However the simulation period is now divided into two different phases The first phase represents the period of active contamination during which the plume grows successively while the second phase represents the period after the source was exhausted During this second phase clean groundwater will enter the model domain over the entire depth of the upstream boundary Furthermore we also include now the simulation of groundwater recharge in this model 3 3 1 Modification of the flow problem For convenience the flow model is similar to the one used in the previous exercise i e the spatial dimensions were modified from the original field setting The details of the flow model and its discretisation are listed in 3 7 To adapt the MODFLOW model from the previous exercise to the present case proceed as follows e Make a copy of the flow model that was readily prepared for this exerc
43. boundary The value of the zone needs to be set to 1 fixed head This will define a fixed head boundary i e the hydraulic head at this grid cell will remain at the value that will be defined as initial head Set initial hydraulic heads by going to Spatial Attributes MOD FLOW BAS6 gt bas 5 Initial Heads and set one_value to 1 The initial heads will be used by MODFLOW as initial estimates for solv ing the flow equations During the solution procedure the correct heads will be computed At locations for which a fixed head boundary condition was defined the heads will remain unchanged as discussed above Table 3 1 Flow and transport parameters used in PHT3D Exercise 1 Flow simulation steady state Total simulation time days 0 24305 Time step days 0 01 Grid spacing m 0 01 Model length m 0 50 Pore velocity m day 0 81 Hydraulic conductivity m day 1 Effective Porosity 0 32 Total Porosity 0 32 Dispersivity m 0 0067 Mineral dissolution and precipitation e Go to Spatial Attributes gt MODFLOW gt LPF gt Ipf 8 to specify the hydraulic conductivity to 1 m day e Go to Spatial Attributes gt MODFLOW WEL and add a zone Well and set the flow rate of the first cell i e the cell representing the inflow end of the model to 0 259 m3 day 3 1 2 Running MODFLOW You have now completed to set up the flow model and you are ready to run MODFLOW This step will provide PHT3D with the flow f
44. boundary m 10 Longitudinal dispersivity m 0 5 Transversal dispersivity m 0 05 3 2 PHT3D Exercise 2 Precipitation dissolution fronts in acid mine drainage This exercise is a second simulation problem that involves mineral disso lution and precipitation reactions as the principle geochemical reactions Compared to the first exercise it is somewhat more complicated as it con tains significantly more aqueous components and mineral reactions Those also occur now in a slightly more complex flow field It was first presented by and later also used as a benchmark problem by In contrast to the pre vious exercise this case includes now also redox reactions The simulations demonstrate the typical hydrogeochemical changes that occur when acidic mine tailings leach into an anaerobic carbonate aquifer Aqueous complexa tion and dissolution precipitation are all considered as equilibrium reactions If the original reaction network defined by is used the simulation includes 17 aqueous components 15 of which are transported 54 aqueous species and six minerals 3 2 1 Spatial discretisation and flow problem The problem is set up as a simple two dimensional flow and transport prob lem with dimensions and properties as defined in 3 4 To construct the MODFLOW model that underlies the transport problem proceed as follows e In iPHTS3D create a new model and save it again in a new separate directory 35 Chapter 3 36 e Under Pa
45. by C D K 2 2 AP BP where capital characters denote the species lowercase symbols indicate the stoichiometric coefficients and K is the equilibrium constant The quantities between the brackets denote the activity of a species Consider for example the dissolution of gypsum CaSO 2H20 Ca SO7 2H20 2 3 The equilibrium constant for this reaction is given by Kap nnn 10 2 4 ial CaSO4 j 2H O i which can be simplified to K gyps Ca SO7 10 ee 2 5 15 Chapter 2 since at low ionic strength see below the activity of water H20 approaches unity and the activity of a pure solid CaSO4 2H2O in this example equals one by definition This expression is referred to the solubility product for gyp sum Note that the equilibrium constant is temperature dependent Keyps 10 74 is valid at 25 C Application of the law of mass action requires that the activities of the species are to be known Generally however we are only provided with concentrations because those are analyzed in the laboratory which do not equal activities because of 1 electrostatic shielding and 2 the formation of aqueous complexes Geochemical models can take these effects into account which will be explained below 2 1 Electrostatic shielding Activities reflect the tendency of ions to react and form a precipitate 7 An ion in solution is surrounded by water molecules and other dissolved ions that act as a
46. ceed with the reactive transport model e Add the reaction module two files supplied for this exercise to the reaction module library See the previous exercise for the details e Set appropriate boundary conditions for the fluxes concentrations across the model boundary e Copy the reaction database that was developed for this problem into the folder where the iPHT3D model is located e Activate the relevant aqueous components that are needed to sim ulate the reactive transport of the chemicals listed in Table 1 of Also include Tmp Temperature For simplicity do not include DOC and ion exchange species From the minerals listed add only Pyrite kinetic version and Fe OH 3 a to the reaction network e Define the initial concentrations background water composition and mineral concentrations in the appropriate iPHT3D Tab Use the water composition listed in Table 1 of For pyrite use an initial concentration of 0 05 mol l bulk volume Note that the oxidation rate of pyrite is among other factors depending on the pyrite concentration Where pyrite is present at higher con centrations the reaction will proceed faster For the temperature Tmp enter a value of 0 017 which corresponds to 1 1000 of the temperature in Celcius e Define the water composition of the injected water as Solution1 Use the water composition measured at the 21 January 1997 as listed in Table 1 of For the temperature add a value of 0 0019 71 Cha
47. cribe the microbial degrada tion of organic compounds depends on the efficiency of the microorgan isms to divert electrons gained in the oxidation step to either biomass generation or towards electron acceptors Assuming that 10 of the carbon is incorporated into biomass the oxidation of toluene under sulfate reducing conditions can be described by the following overall reaction CrHg 0 14NH 4 15907 2 58H20 1 86Ht gt 6 30HCO 0 14C H70 N 4 15H9 4 7 Equation 4 7 shows that the complete mineralization of 1 mol of toluene consumes 4 15 mol of sulfate during growth and yields 0 14 mol of sulfate reducing bacteria represented by the genralized chemical formula for biomass C H7O2N Therefore X o 0 14Um 4 8 and OC sui f 4 15um 4 Ot 5v 4 9 These equations can be implemented in PHREEQC to compute the temporal development of toluene sulfate and the microbial mass in a batch experiment Tables 4 1 and 4 2 list the initial conditions and rate constants for the simulation 53 Chapter 4 The input file toluene_degradation phrq is available and can be opened in PHREEQC for Windows Two new aqueous species are defined in the input file that represent toluene species name Toluene and X species name Sulfred since X is already reserved for exchange species The rate expressions for toluene degradation and biomass decay are defined under the keyword RATES They are programmed in BASIC which ca
48. d i e 1 x 107 cm s is assumed and the vertical anisotropy ratio KV KH is 0 01 The resulting uniform seepage velocity is approximately 0 10 m day 1 1 2 Solute Transport The contaminant source is represented by a single cell with a constant con centration of 1 0 located at row 16 column 1 in the first model layer as shown on Figure 1 1 Clean water enters along the remaining portion of the upstream boundary Solute is free to exit from the constant head cells along the downstream boundary The boundary conditions along the other model faces are implicitly set as zero mass flux The transport parameters are listed in Table 1 1 and are assumed to be uniform for all model layers Chapter 1 Contaminant Source Figure 1 1 Illustration of the conceptual model and numerical model grid 1 1 3 Setting up the flow model To set up the flow model proceed as follows e In PMWIN create a new model and save it in a new and separate directory folder e Select Grid gt Mesh size to specify the model domain as described in Table 1 1 Set the model s Top elevation to 6 m e Select Grid Layer property and simply click OK We can leave it as is but PMWIN requires that all menu entries for specifying model input are clicked e Select Parameters Time and set the Simulation Time Unit to days Then define the simulation time Divide the total simulation time of 10 years into a first stress period of 1095 days
49. d contaminants and inorganic groundwater constituents within the aquifer Therefore the identification and quantification of the vertical mixing process and how to model it is generally rather difficult The present exercise is loosely based on the field e g and corresponding modelling studies e g that provided a detailed investigation on the fate of hydrocarbon compounds at a contaminated sels MP11 MP2 MP6 MP7 14 24 e e 34 e e 44 a Benzene n r gt 10 000 ug L z J 1 000 10 000 pg L E i ug l o o B 10 1000 ug L Q g x lt 10 ug L 2 5 1 e T e 24 e a e e e Q 34 e 3 44 b Toluene s a e e e g 5 a a 0 e 1 e 2 e 3 lt 10 mg L m 4 Sulphate 10 20 mg L 7 gt 20 mg L gt 5 ki 0 50 100 150 200 250 300 Distance along the plume m Figure 4 1 Measured benzene toluene and sulphate concentrations vertical cross section at the Eden Hill field site after 59 Chapter 4 site in Perth Western Australia At this site detailed information of hydrochemical parameters have been intensively recorded The data show that toluene ethylbenzene and xylenes were mineralised under sulfate reducing conditions in a seasonally varying groundwater flow field while benzene was shown to be persistent Figure 4 1 show
50. e field but pH and total carbon are known from the laboratory analysis The third option can be used to check the measured pH of a water sample in which both alkalinity and total car bon were analysed Options 4 and 5 are useful when the CO pressure of the water sample is known 87 Appendix II Common causes for model crashes One of the most frustrating things in learning a new code is getting stuck because the model crashes for unclear reasons Most of the time the solution is very simple So simple that overlooking it is easier than finding it The following checklist may be useful in tracking down the cause of the model failure e Model names Avoid using model names that contain blanks causes run time error 53 in PMWIN e PMWIN run time error 5 This one is a bit of a mystery Copying the entire folder with the model files to a different folder will probably resolve it e Units Are the values for aqueous concentrations and ion ex changer sites provided in mol l and in mol lb for minerals e Time stepping Has the temporal discretization that determines the number of PHREEQC reaction steps been defined The num ber of reaction steps should be selected such that solutes are not transported much further than one grid cell e Charge balance of water compositions Are aqueous solutions that are defined as initial water composition s and at model bound aries recharge wells constand head cells charge balanced e
51. ent options for specifying TIC alkalinity and pH are listed in appendix I e What is the percent error of the electrical balance of this chemical analysis Answer Because its absolute value is less than 5 we may conclude that the analysis is acceptable e The label lists no concentration for silica Si How would the calcula tion of the electrical balance change if we include this element e How much is the ionic strength equation 2 7 that PHREEQC has calculated Answer e What are the 2 most dominant Ca species and their concentrations Answer Ten easy concentration Do pesrik concentration 2 3 Saturation state The saturation index is an important parameter in aqueous geochemistry as it provides information on the minerals the groundwater has been in contact with and on which minerals are likely to precipitate or dissolve An expression analogue to equation 2 5 can be written that uses the actual activities of the species in the groundwater instead of the activities at equilibrium which is referred to as the ion activity product IAP TAP 55 Ca S027 2 10 A comparison of AP and K provides information on the saturation state of the groundwater for a mineral This is commonly done by calculating the saturation index ST for example for gypsum TAP ays Sloss loz ee 2 11 yps At equilibrium IA Pyyps equals Kgyps so SIgyps 0 Gypsum potentially dissolves in g
52. es reflect the changes in their concentration ratio Both Fe 2 and Fe 3 concentrations are strongly controlled by the pH of the solution Generally their solubility increases with a decreasing pH However as the increase of their solubilities is not proportional the 40 Mineral dissolution and precipitation Fe 2 Fe 3 ratio can change and so does the pe 3 2 8 Model variants e To see the effect of calcite buffering rerun the model with a much smaller initial calcite concentration and compare how far the low pH zone and dissolved aluminium are migrating in this case e Test the effect of calcite filled permeable reactive barrier PRB To do that add a small zone between x 40m and x 50 m and y 0m and y 10m where a high concentration of calcite is present but no other minerals After rerunning the model inspect the evolution of the low pH zone and of the aluminium migration 41 Chapter 3 3 3 PHT3D Exercise 3 Transport and nitrifica tion of ion exchangeable ammonium This modelling example is based on a field site contamination prob by product of the production of smokeless fuel has polluted ground water over several decades A reactive transport modelling study that integrated and reproduced the major processes that were believed to occur at the field site was presented by One of the key features ob served at the site is the strongly retarded migration of ammonium and the geochemical footprint that was left
53. ew results with the previous results model run without surface complexation load the previously saved break through curves ASCII files as observations for As 5 at the ob servation wells Chemical kinetics 4 4 PHT3D Exercise 6 Degradation of chlorinated ethenes and isotopes 4 4 1 Modelling Scenario In this exercise we investigate the sequential degradation of chlorinated ethenes and the use of stable isotopes to provide forensic evidence for the origin of elevated vinyl chloride concentrations in a pumping well The principle modelling scenario consists of an aquifer that is con taminated by two different contaminant spillages both positioned be low sites of former metal processing companies Company E and Com pany W The modelling scenario is set up for a vertical cross section The two different contamination sources consist of NAPL pools that are located on top of clay lenses One of the NAPL sources is assumed to consist of PCE tetrachlorethylene while the other is assumed to consist of a mixture of TCE trichlorethylene and toluene Groundwa ter is entering the model domain through the Eastern fixed hydraulic head 13 7 m a s 1 and Western boundary fixed hydraulic head 14 5 m a s l Under natural conditions the groundwater discharges to the river prescribed hydraulic head 10 8 m a s l1 However there is also an active pumping well located between the PCE source and the river The majority of the water that en
54. familiarizes the reader with PHREEQC by using simple example exer cises Each successive chapter discusses additional theory and the exercises increase in complexity After having worked through the entire course guide the student has enough theoretical background and practical experience to build his or her own models and apply these to real world problems Chapter 1 Single Species Transport 1 1 MT3DMS Exercise 1 Part 1 Single Species 3 D Transport Simulation The purpose of this exercise is to learn the use of MODFLOW and MT3DMS for analysis of 3 D transport of singe species contaminants under natural gradient hydraulic conditions The effect of different transport solution tech niques on plume migration are examined In addition the impact of equi librium sorption and first order decay will be simulated in the second part of this exercise 1 1 1 Groundwater Flow A sketch of the system being modeled is shown on Figure 1 1 The aquifer is confined with a uniform thickness of 6 m that is represented by three 2 m thick layers The model is discretized into 31 rows and 50 columns with a regular spacing of 5 m in the column and row directions Groundwater flow is assumed to be steady Hydraulic heads of 18 33 m and 10 0 m are assigned along the upstream and downstream faces of the model respectively so that there is a uniform ambient hydraulic gradient of 0 034 along the x axis A uniform horizontal hydraulic conductivity of 0 864 m
55. ffer ent electron donors This can be necessary when for example ben zene degrades more slowly in a BTEX benzene toluene ethylbenzene xylenes plume than the other compounds In this exercise you will model a case of a simultaneous uptake of BTEX constituents by one microbial group Different initial amounts of organic compounds are present Table 4 3 and the initial sulfate mass is increased compared to the first example Table 4 3 Initial concentration for simulation of a batch experiment of BTEX degradation Component Concentration mol l Benzene 4 0 1074 Toluene 3 0 1074 Ethylbenzene 1 0 1074 Xylene 2 5 1074 N 1 0 1074 s 6 5 0 1074 Sulfred 1 0 1078 Sulfred sulfate reducing bacteria Table 4 4 Rate constants for simulation of a batch experiment of BTEX degradation Component Constant Ksuiz mol l 1 0 107 Koenzene mol l 1 0 10 Ktotuene mol l 10 1075 Kethylbenzene mol l 1 0 1075 Kzylene mol l 1 0 1075 Umazsz sulf benzene 1 day 0 1 Umazsx sulf toluene 1 day 5 0 Umaz sulf ethylbenzene 1 day 0 5 Umaz sulf xylene 1 day 1 0 vaee 1 d 0 1 The corresponding reaction equations are for benzene 57 Chapter 4 CoHe 0 12NH 3 459077 2 64H20 1 38Ht gt 5 4HCO 0 12C3H7O2N 3 45HeS 4 12 and for ethylbenzene and xylene s CgH 0 0 16NH 4 85SO7 2 52H2O 2 34Ht gt 7 2HCO z 0 16C5H7O2N 4 85H2S 4 13 The PHREEQC input file for this exercise is already prepared
56. ich are released or con sumed per mole kinetic reactant In this example for case a 0 14 moles 55 Chapter 4 of sulfate reducing bacteria are produced for each mole of toluene re acted which removes 0 14 moles of C sH7O2N from the solution A negative stoichiometric coefficient indicates that a species is removed from the solution for each mole of kinetic reactant formed which is why the stoichiometric coefficient for toluene is 1 The keyword USER_GRAPH control the graphical output of this simulation If you go to the Chart tab in PHREEQC for Windows and start the calculations you can even see the chart being filled with data points during the calculation The BASIC statements GRAPH_X GRAPH_Y and GRAPH_SY allocate numbers to the x y and sec ondary y axes respectively The identifiers are basically self explanatory Notice that two cases are considered in this example case a and case b Case b differs from case a in that the effects of the different valence state of bacteria compared to the end product CO2 and the related geochemical changes are not considered All sulfate is consumed during bacterial growth but none during bacterial decay Run the model for case a first and look a the chart Notable is the lag period of several days before the degradation affects the aqueous concentrations of toluene and sulfate Its length depends largely on the initial bacterial concentration and on Umaz the maximum uptake rate values given
57. ield that is used to compute advective dispersive transport of mobile species components e Under Parameters Flow use the Pencil Button to write all the input files required to run MODFLOW e Under Parameters Flow use the Red Arrow button to start the MODFLOW simulation During the execution of MODFLOW a number of output files are cre ated including the specific file mt3d flo which contains the flow vectors of each cell for each time step in the simulation This file is required for the subsequent PHT3D simulation s and PHT3D will not run until the file mt3d flo exists If changes are made to the MODFLOW input parameters MODFLOW needs to be rerun such that the file mt3d flo is updated 3 1 3 PHTS3D reaction definition We continue now with the setup of the reactive transport model In some cases the next step would be to prepare a problem specific reaction module However for simpler problems such as the one in this exercise and many others that only include equilibrium reactions this is not the case All of the aqueous species components and minerals needed to simulate this LEA based reactive transport problem are already included in the original PHREEQC 2 database This means that we don t have to define our own set of equilibrium reactions but we can simply use a PHT3D database file that is equivalent to the original PHREEQC 2 database To specify this and to incorporate the correct reaction database proceed as follows e C
58. in tables 4 2 The removal of the initial toluene mass 0 1 mmol is reached after 14 days at a time when approximately 0 415 mmol of sulfate are depleted The microbial net growth then stops immediately Um 0 and the microbial mass is subsequently changing at the rate given by equation 4 6 thereby consuming the remaining 0 35 mmol of sulfate You can switch to case b by switching the signs in front of the formula identifier under KINETICS Now all sulfate is consumed during bacterial growth but none during bacterial decay Run the sim ulation and observe the resulting chart 4 1 2 PHREEQC Exercise 5 Kinetic BTEX degradation Of course the formulations for microbial growth above apply only to the uptake of a single substrate whereas contamination often involves numerous organic compounds The possibly simultaneous uptake of these substrates can be incorporated into the modeling approach de scribed above The model of states that the growth of the degrading microbial community is simply the sum of the growth rates arising from degradation of individual organic contaminants Ox Ok OL 5 oe ta X 4 10 56 Chemical kinetics where in analogy to equation 4 5 each of the growth terms can be derived from OXy org n ea Vas 2 ot org n T Coren Kea T Cea 4 11 The uptake rates v can differ between different substrates In this way it is possible to model varying degradation rates of di
59. ined by as illustrated in Figure 3 1 Generally the simulated concentrations including the positions of the mineral fronts agree very well However the PHT3D simulated pH near the inflow end of the model domain lies slightly above the pH simulated by MST1D 33 Chapter 3 x10 25 11 CI PHT3D 2 CI MST1D 10 L15 9 3 pH PHT3D o 1 8 pH MST1D 0 5 7 o 6 0 0 1 0 2 0 3 0 4 0 0 1 0 2 0 3 0 4 15X10 y 5x10 Mg PHT3D Ca PHT3D Mg MST1D 4 Ca MST1D T Ls 5 5 E X O05 oO 1 0 o o 0 1 0 2 0 3 0 4 0 0 1 0 2 0 3 0 4 x10 x10 25 25 Dolomite PHT3D 2 2 Dolomite MST1D 5 a5 iret Calcite PHT3D ka i 4 Calcite MST1D 4 o o q 0 5 0 5 I 1 o 0 o 0 1 0 2 0 3 0 4 0 0 1 0 2 0 3 0 4 Distance m Distance m Figure 3 1 Simulation results for PHT3D and MST1D Aqueous and mineral concentrations after 21000 s 34 Mineral dissolution and precipitation Table 3 4 Flow and transport parameters for the simulation Flow simulation steady state Total simulation period days 2000 Model length m 100 Model thickness m 10 Grid spacing Ax m 4 Grid spacing Az m 1 Porosity 0 35 Horizontal hydraulic conductivity m day 1 Vertical hydraulic conductivity m day 1 Piezometric head upstream boundary m 12 Piezometric head downstream
60. initial water composi tion is in chemical equilibrium and the aqueous solution is charge balanced If the aqueous solution is not in equilibrium equilibrium conditions will be adjusted within the first reaction step at the end of the first time step This can cause an undesired change for example of the solution pH or the disso lution precipitation of minerals To avoid this it is recommended to study and charge balance the aqueous solution with PHREEQC before entering the data into iPHTS3D Table 3 2 Aqueous concentrations used in PHT3D Exercise 1 Aqueous component Cint Cin flow mol 1 mol 1 pH C CHL 70 pe 4 0 4 0 C 4 1 23 x 1074 0 0 Ca 1 23 x 1074 0 0 Mg 0 0 1 0 x 1073 Cl 0 0 2 0 x 1078 Table 3 3 Mineral concentrations used in PHT3D Exercise 1 Mineral Cinit mol 171 Calcite CaCO3 3 906 x 10 5 Dolomite CaM g CO3 2 0 0 There are various ways and types of boundary conditions to define the hydrochemical composition of water that is entering the model domain dur ing a simulation In the present case the inlet inflow water composition is defined by providing the aqueous component concentrations for the injection well located at the first cell The inlet water contains 0 001 mol I of Mg and 0 002 mol I of Cl In contrast to the initial water composition no C 4 and no Ca is contained in the inflow water You can enter these values see also Table 3 2 into iPHT3D as follows To define both the
61. ise and save it in a new folder The model files are located in the folder 2d_flow_problem_haerens e Select Parameters and select the Time button and change the Total Simulation Time to 3000 days and select a Step Size to 40 days 42 Mineral dissolution and precipitation Table 3 7 Flow and transport parameters used for PHT3D Exercise 3 Flow simulation steady state Total simulation time days 3000 Time step size 40 Model length m 100 Model thickness m 10 Grid spacing Ax m 4 Grid spacing Az m 0 5 Porosity 0 35 Horizontal hydraulic conductivity m 1m d Vertical hydraulic conductivity m 1 m d Piezometric head upstream boundary m 12m Piezometric head downstream boundary m 10 m Longitudinal dispersivity m 0 5 Transversal dispersivity m 0 05 e To add groundwater recharge we first need to activate the MOD FLOW recharge package such that it is availale within iPHT3D To do this go to iPHT3D s top menu and select Add in gt Mod flow Modules You will see a list of all available MODFLOW packages and you will see which ones are already activated Add the recharge package by ticking the box near RCH Now go to Spatial Attributes MODFLOW RCH rch 2 Recharge Value to specify a rechare rate of 0 001 m day i e 1 mm per day Don t forget to press ok once you have entered the value e Now rerun MODFLOW to regenerate the flow file mt3d flo which will later be used by the transport model 3 3 2 Data
62. ive transport 2 000 304 Running PHV3D sa i e 220 4 4 8404 eo es 3 4 Surface complexation 2 00200 08 3 4 1 PHREEQC Exercise 4 calculation of a charged sur face COMPOSILION 6 4400 S maup hae ee Ee ee at Modelling kinetically controlled reactions 4 1 Modelling bioremediation processes 4 1 1 PHREEQC Example Kinetic toluene degradation 4 1 2 PHREEQC Exercise 5 Kinetic BTEX degradation 4 2 PHT3D Exercise 4 Microbially mediated petroleum hydro carbon degradation under sulfate reducing conditions 4 2 1 Setting up the flow problem 4 2 2 Simulating nonreactive transport 4 2 3 Simulating the dissolution from a multi component NAPL source gt aaacasa opia gaa b iaaa aiia 51 52 53 56 62 4 3 PHT3D Exercise 5 Pyrite oxidation during deep well injection 69 4 3 1 Setting up the flow problem 4 3 2 Setting up the nonreactive transport problem 4 3 3 Setting up the reactive transport problem 4 3 4 Seasonally changing redox zonation 4 3 5 Mobilzation of arsenic 0 4 4 PHT3D Exercise 6 Degradation of chlorinated ethenes and OLODE 4 reeni hand deo Ob BS aA SR e Ea g 4 4 1 Modelling Scenario 204 4 4 2 NAPL source zones 02 0008 69 4 4 3 Formation of dissolved plumes 4 4 4 Degradation reactions 04 4 4 5 Chlorinated ethenes at
63. lity changes during subsurface passage Pyrite oxidation will be defined as a kinetic process in which the reac tion rate depends on the water temperature For simplicity we include temperature as a seperate aqueous mobile component called Tmp The reaction rate expression expression has been programmed such that the value of the component Tmp is read and used during the computation of the reaction rate 4 3 1 Setting up the flow problem To simplify the original fully three dimensional model for this exer cise the flow model is defined as a two dimensional model for a single stratigraphic layer Aerobic surface water is injected through an injec tion well and extracted at an extraction well located 100 away from the injection well e Use the data supplied in Table 4 9 to set up the flow model Note 69 Chapter 4 Table 4 9 Flow and transport parameters used for the deep well injection problem Flow simulation steady state Total simulation time days 360 Time steps 180 Model length column direction m 200 Model width row direction m 80 Grid spacing Ax m 10 Grid spacing Ay m 10 Total Porosity 0 35 Effective Porosity 0 35 Horizontal hydraulic conductivity m 10 Piezometric head upstream left boundary m 20 Piezometric head downstream right boundary m 20 Depth target aquifer m below sealevel between 300 and 310 x Position injection well m 145 y Position injection well m 5 x Position recovery well m
64. n be interpreted by PHREEQC so that any type of user defined rate expression can be modelled RATES Toluene_sulf start 1 k_Tolu le 05 2k Sulf 1e 05 4 v_up_max 5 0 86400 10 mTolu TOT Toluene 22 mSulf TOT S 6 30 mSulfred TOT Sulfred 32 IF mTolu lt le 08 THEN GOTO 200 34 IF mSulf lt le 09 THEN GOTO 200 40 mon_Tolu mTolu k_Tolu mTolu 50 mon_Sulf mSulf k_Sulf mSulf 80 growth v_up_max mon_Tolu mon_Sulf mSulfred Table 4 1 Initial concentration for simulation of a batch experiment of toluene degradation under sulfate reducing conditions Component Concentration mol l Toluene 10 1077 N 1 0 1074 S 6 4 5 1074 Sulfred 1 0 107 T Sulfred sulfate reducing bacteria Table 4 2 Rate constants for simulation of a batch experiment of toluene degradation under sulfate reducing conditions Component Constant Ksuif mol 1 1 0 107 Ktotuene mol 1 1 0 107 Umaa 1 day 5 0 Udec 1 day 0 1 54 Chemical kinetics 100 rate growth 110 moles rate TIME 200 SAVE moles end Sulfred start 10 v_decay 0 1 86400 15 c_min le 08 20 mSulfred TOT Sulfred 30 IF mSulfred c min gt 0 THEN inhibit_fac_decay mSulfred c_init mSulfred 35 IF mSulfred cmin 0 THEN inhibit_fac_decay 0 40 IF mSulfred c min lt 0 THEN inhibit_fac_decay 0 190 rate v_decay mSulfred inhibit_fac_decay 200 moles rate TIME 220 SAVE mole
65. n when the effect of electrostatic shielding is taken into account as well Clearly neglecting the difference between total concentration and activ ity in the calculation of the saturation indices of seawater for minerals would result in overprediction of the values The following PHREEQC input file repeats the calculations by to calculate the saturation index for calcite in seawater using the seawater composition which is taken from example 1 in the user s manual of PHREEQC SOLUTION 1 Seawater units ppm pH 8 22 temp 25 0 Ca 412 3 22 Specation and saturation indices Mg 1291 8 Na 10768 0 K 399 1 Cl 19353 0 Alkalinity 141 682 as HCO3 S 6 2712 0 END You can run this input file and verify that the total concentrations of Ca t and CO amount to Mo 2 9 504 107 and Moo 3 826 1075 3 mol kgw respectively The corresponding activities are Ca t 2 380 1073 and C037 7 969 107 mol kgw respectively Inserting the activities in the expression for the saturation index equation 2 11 yields calcite 2 1 el TAP ca toe 380 10 3 7 969 10 N aes SIealcit log 1078 48 K calcite 2 12 which is the value that is reported by PHREEQC for calcite in under the heading Saturation indices Note that if the concentrations of the species had been used instead of the activities the calculated saturation index would be STeatcite 2 04 which is much too high 2 3 2 PHR
66. nhibition term is computed from bs 4 17 Um Umaribio TEETER Trio Obio mar tot 4 18 Obie sade The above equation was used as a base for the computation of the degradation rates with PHREEQC PHT3D Tolu start 10 mSulfred tot Sulfred 20 mSulf tot S 6 30 mTolu tot Tolu 35 if mTolu lt 1e 10 then goto 220 40 if mSulf lt le 10 then goto 220 68 v_up_max_Tolu parm 1 86400 84 theta_biomax parm 2 90 inhib_bio theta_bio_max mSulfred theta_bio_max 95 monSulf mSulf 0 0001 mSulf 100 monTolu mTolu 0 0001 mTolu 150 v_up_Tolu v_up_max_Tolu x monTolu monSulf 160 v_up_Tolu v_up_Tolu inhib_bio 190 rate v_up_Tolu mSulfred 200 moles rate time 220 SAVE moles 66 Chemical kinetics end To proceed and include the biodegradation reactions in the simula tions Add the additional species to the reaction network Define the background concentrations of the aqueous components and for the recharge water composition according to Table 4 7 Enter a value of le 07 mol l as initial concentration for the sulfate reducing bacteria Check if the boundary conditions for the transport are defined correctly such that the inflowing water at the upstream boundary has the same composition as the background and the recharge water Allocate the reaction rate parameters listed in Table 4 8 To enter the values select Select Models PHT3D gt
67. nite Difference and MOC respec tively and export the calculated concentrations to additional ASCII files Note that dispersion and sink source terms are always solved using the stan dard finite difference method and the GCG implicit matrix solver Use Excel to graph the three breakthrough curves The TVD and MOC solutions are nearly identical as both solutions have minimal numerical dispersion The Upstream Finite Difference solution on the other hand contains a signif icant amount of numerical dispersion In addition check the mass balance 11 Chapter 1 errors for each solution TVD and FD should have close to zero mass bal ance error while MOC has a large mass balance error at the beginning but diminishes quickly which is typical of a particle tracking based solution technique 1 1 9 Optional Numerical Experiment If you have time left you can freely experiment with the model you have set up For example if you increase the longitudinal dispersivity from 1 to 5 m and rerun the three cases do you see all three concentration curves become nearly identical If so why 1 22 MT3DMS Exercise 1 Part 2 Single Species 3 D Transport with with Sorption and First Order Decay 1 2 1 Linear equilibrium controlled sorption We use the second part of the exercise to introduce some simple reactions Initially we only consider sorption Proceed as follows e Make a copy of the folder in which the existing model is located e Sta
68. oceed as earlier in the exercise when executing MOD FLOW e Go to Parameters Chemistry and use the Pencil Button to write all the input files required to run PHT3D 46 Mineral dissolution and precipitation e Then go to Parameters Chemistry and use the Plume Button to start PHT3D After the end of the simulation visualise the ammonium plume and check if the simualtion results are plausible Plot also a breakthrough curve BTC for the ammonium con centration at x 50 m and z 7m To compare the results with observed data go to the main menu top row to Add in Batch and add the content of the pre prepared python script specIn data self importTabFile ex3BTC1 tat tlist core getTlist2 iper range len tlist species Amm Ca K Mg Na time conc labs core onPtObs BO iper Chemistry BTC species pylab figure nrows ncols 3 2 for i spec in enumerate species pylab subplot nrows ncols i col0 specIn index spec pylab plot data 0 data col0 o pylab plot tlist conc il pylab legend spec pylab draw This will show your simulated and the observed data in compari son To improve the agreement between observed and simulated data we now consider ion exchange reactions To include this process in the reation network and to estimate a suitable value for the cation exchange capacity CEC of the exchanger sites proceed as follows
69. of the four BTEX compounds and that the remaining fraction can be represented by a fifth compound that is assumed insoluble The computation of benzene dissolution within such a multicomponent NAPL mix can be programmed as follows Benznapl start Compute kinetic dissolution from multicomp NAPL 10 mBenznapl tot Benznapl 15 if mBenznapl lt le 10 then goto 200 20 solub_Benz 0 022820 25 mBenz tot Benz 32 mTolunap tot Tolunapl 34 mEthynapl tot Ethynapl 36 mXylnapl tot Xylnapl 38 mLowsolubnapl tot Lowsolubnap 40 m_napltot mBenznapl mTolunapl mEthynapl 41 m_napLtot m_napLtot mXylnapl mLowsolubnapl 42 if m_napLtot lt le 10 then goto 200 50 msolub_Benz mBenznap m_napltot solub_Benz 60 rate parm 2 24 3600 msolub_Benz mBenz 70 moles rate time 80 if moles gt m then moles m i suggest to remove that line that slows the computation can lead to strange results and is done by default 200 save moles end The above rate expression as well as those for the dissolution of toluene ethylbenzene and xylene have been included into the database for this exercise The values used for C are listed in Table 4 6 Note the much higher solubility of benzene which is considered the most toxic of the BTEX compounds 64 To start the reactive transport modelling e Copy the file pht3d_datab ex_eden_hill which contains the site
70. olute temperature The intrinsic dissociation constant Kint applies to the chemical binding of the ion and the surface The PRHEEQC database contains a compilation by Dzombak and Morel 1990 of laboratory determined values of intrinsic dissociation constants for hydrous ferric oxide The apparent dissociation constant Ka however varies with the potential and thus the charge at the surface In order to calculate the K the Wo needs to be known PHREEQC uses a double layer model to relate the charge density on the surface os with ionic strength J equation 2 7 and Wo 7 log K log King 3 1 3 4 1 PHREEQC Exercise 4 calculation of a charged surface composition This exercise demonstrates the use of the SURFACE keyword which is used in PHREEQC for surface complexation calculations The input file has already been prepared and is called kd_surf phrq The definition of the surface is as follows SURFACE 1 Hfo_w 2e 4 600 0 088 Hfo_s 5e 6 equilibrate 1 The names Hfo_w and Hfo_s correspond to the names of the mas ter species in the PRHEEQC database that refer to the weak and the 49 Chapter 3 strong surface sites respectively The number that follows is the num ber of surface sites in moles For the weak sites the surface area 600 m g and mass of ferrihydrite 88 mg are also specified Unless these values are explicitly typed for Hfo_s as well they apply to both the weak and the strong sites The inp
71. omplexed ions are updated from the mass balance equations step 5 The newly found values are used to obtain a more reliable value for the ionic strength and iteration continues until the result no longer changes significantly 2 2 1 Your first PHREEQC input file To do calculations you need to tell the program what it has to do This is done via input files that contain a list of instructions So called keywords are used that have an intuitive meaning The keyword SOLUTION for example is used to define the measured composition of our water sample Within the keyword data blocks identifiers are used to set various options For example the identifier units can be set to mg l if the reported concentrations are in mg l Let s look at an example to make things more clear Say that we have some pure water in which we have dissolved 1074 moles of the mineral fluorite so the only dissolved ions are Ca and F and we want to know the distribution of species in the water sample This is then what your input file looks like SOLUTION 1 units mmol kgw these are the default units Ca 0 1 F 0 2 END This input file tells PHREEQC that we have a solution which is given the number 1 the significance of these numbers will become clear later The units of the concentrations are specified in mmol kgw The concentra tions are simply entered after the element names The input datablock is terminated by the keyword END which tells PHREEQC to sto
72. on of dissolved plumes When groundwater passes through these zones dissolution mass trans fer to the aqueous phase will occur therefore creating a PCE plume on the Western side of the river On the Eastern side of the river the NAPL source creates a TCE as well as a toluene plume During the simulation period of 10 years these plumes grow initially but eventually become stable given that we assume steady state flow and constant source release conditions This means that chemical con centrations at a particular location are not changing any further after a specific simulation time 4 4 4 Degradation reactions Besides the plumes of the primary contaminants there are also plumes of the transformation products developing These result from the degra 76 Chemical kinetics dation reactions that may proceed in suitable hydrochemical and mi crobial conditions Under the oxic conditions found in this river valley aquifer the toluene plume released from the NAPL source undergoes aerobic degra dation thereby consuming oxygen On the Eastern side of the model domain this creates locally anaerobic conditions between the NAPL source and the groundwater discharge locations e Visualise the contours of the toluene plume Tolu after 10 years 3650 days simulation time e Visualise the oxygen concentration contours after 10 years 3650 days simulation time Where these anaerobic conditions prevail reductive dechlorination of TCE may t
73. opy the database file that is going to be used for the PHT3D simu lation into the folder that contains all files for this exercise and name the file pht3d_datab dat e Go to Parameters Chemistry and select the Import database button iPHT3D will then read and interpret the PHT3D database file 29 Chapter 3 e Then click on the Chemical Reaction button This will open up a new dialog box e Within the new dialog box go to the Solution Tab and select the aque ous components that you want to include in the model In this exercise the simulation needs to include C 4 Ca Cl Mg pe and pH e Then click the Phases Tab and select i e activate the two minerals Calcite and Dolomite These steps are the only steps needed to define the set of chemical reac tions that will be considered in the PHT3D simulation From this informa tion iPHT3D will create the interface file that controls the communication between the MT3DMS part and the PHREEQC 2 part of the PHT3D sim ulator This file pht3d_ph dat contains the records which define that the reaction network contains 6 equilibrium aqueous components including pH and pe and two equilibrium minerals The order in which the components are listed determines which species number will be allocated to each entity included in the simulation For each of the first MCOMP 2 C 4 Ca Mg Cl of the NCOMP components advective dispersive transport steps will be simulated by the appropriate MT3DMS r
74. or the recharge water in the zone between x 15m and x 30m e Run the simulations and evaluate the plausability of the results 4 2 3 Simulating the dissolution from a multi component NAPL source The first part of the reactive transport simulation considers only the mass transfer process in a NAPL containing source zone towards the dissolved i e aqueous phase and the subsequent advective dispersive transport of the dissolved mass In the source zone dissolution of NAPL compounds e g acts as a time varying contamination source for the passing groundwater The rate at which mass is trans ferred from the NAPL into the aqueous phase is a function of e Interfacial area between the NAPL phase and the aqueous phase 62 Chemical kinetics e Extent and morphology of the source in particular the maximum cross sectional area perpendicular to the main groundwater flow direction e Groundwater flow velocity e Solubility of individual hydrocarbon compounds e Composition of the NAPL source mole fraction Notwithstanding the relative importance of each of these factors the concentration of individual hydrocarbon compounds in the ground water will at the local scale within the contamination source zone very often reach an equilibrium concentration that is equal or close to the multi component solubility of the compound This multi component solubility Ci ame is described by Raoult s law Cw C8 ym 4 14
75. oundwater 81 Chapter 4 Flow Systems A Case of Nitrate Removal by Oxidation of Pyrite Water Resour Res 28 10 2829 2843 Fetter C W 2001 Applied Hydrogeology Prentice Hall USA Garrels R M Thompson M E 1963 A chemical model for sea water at 25 C and one atmosphere total pressure American Journal of Science 260 57 66 Grathwohl P Klenk I D Eberhardt C Maier U 2000 Steady state plumes mechanisms of transverse mixing in aquifers In John ston C Ed Contaminated Site Remediation from source zones to ecosystems Proc 2000 CSRC Melbourne Vic 4 8 Dec 2000 pp 459 466 Guerin M Zheng C 1998 GMT3D Coupling multicomponent three dimensional transport with geochemistry In Poeter E P Zheng C Hill M C Eds MODFLOW 98 pp 413 420 Proc International conference Colorado School of Mines Golden CO Haerens B 2004 Reactive transport modelling of a groundwater con tamination from a former coking plant Departement of Geografie Geologie KU Leuven Belgium PhD thesis Jones I Davison R M Lerner D N 1998 The importance of understanding groundwater flow history in assessing present day groundwater contamination patterns a case study In Lerner D N Walton N R G Eds Contaminated land and groundwater Fu ture directions Engineering Geology Special Publication The Geo logical Society London pp 137 148 Jones I Lerner D N 2001
76. outines For component numbers MCOMP 1 pH and MCOMP pe no transport step is carried out The entities MCOMP 1 NCOMP here Calcite and Dolomite are immobile and no transport simulation will be carried out for them 3 1 4 Initial and inflow concentrations The next step in setting up this problem is to specify the initial concen trations that define the hydrogeochemistry of the aquifer at the start of the simulation Time 0 The information provided here will be trans lated by ipht3d into the the basic transport package file which is called pht3dbtn dat in PHT3D The initial concentrations that need to be en tered into iPHT3D as defined in the paper by are listed in Table 3 2 for aqueous components and in Table 3 3 for minerals Note that aqueous concentrations are always defined in units of mol l t In contrast the unit for the initial concentrations of min erals is NOT mass per volume of water i e mol but is defined as mass per bulk volume i e mol Le defined their mineral concentrations as mass per mass of soil i e mol en and provided the bulk density 1800 kg m of the soil Therefore their initial concentrations for calcite of 2 176 x 107 mol ko translates to 3 906 x 1075 mol 171 which needs to be used in PHT3D volume 30 Mineral dissolution and precipitation If a simulation problem like the one we are simulating here is a pure equilibrium problem it should be warranted that the
77. ow expand the last simulation and include additionally the two aqueous components As 5 and As 3 as well as the kinetic minerals Fe OH 3 a and Arsenopyrite in the reaction network 74 e Set the initial concentration for Arseonpyrite to 0 01 mol lbutk e Set the reaction rate constants for Fe OH 3 a to 2 x 107 and define that the ratio of Arsenopyrite Pyrite dissolution is 0 001 i e parm 1 0 001 for Arsenopyrite e Run the simulation and plot a breakthrough curve for As 5 and As 3 at the extraction well and at a location half way between injection and the extraction well Save the breakthrough curves from the observation wells to ASCII files e Now include surface complexation reactions by copying the file postfix phrq into the folder in which your model is located This file will be read by PHT3D whenever it is present In this way some instructions and definitions can be quickly and simply sup plied outside the Visual Modflow environment Get the file from the course CD or download it from http www pht3d org course ex5 postfix phrq Open it with PHREEQC for Windows to inspect the input in structions e Run the model and inspect again the As 5 and As 3 break through curves at the extraction well Note convergence is only achieved if the diffuse_layer option is included under the KEY WORD SURFACE in the file postfix phrq As a result the model execution time will increase significantly e To compare the n
78. p reading the input file and perform the calculations for this input datablock The sign indicates a comment All text that follows on the same line is ignored when the input file is read To start the calculations click Calculations gt Start in the main menu or press the green arrow button on the toolbar The program will now show a window with information on the progress of the calculations When the calculations are finished the caption of the button on the bottom end of this window changes from Cancel to Done figure 2 2 Press it to go to the output editor where you will find the results of your calculations 18 Specation and saturation indices FREEGE what de yow want to mute to rssi passid Eaves Siam ienser T alej E l Oxo ireas Oeatace Cred Cha aOLUTICN 1 unite yaol kgw these sre the default units kos fie 0 Piogen Fle Ptessgeunttied Dupa fie C Progam Fes Phasegoluarttled on Dotshsie tie C Prog on Flay Ptessac Piveene dy End ot Run 4 Nodied Iremt Figure 2 2 Screenshot of the input editor and the progress window after the calculations have finished The panel to the right of the input editor shows the keywords and identifiers that can be used in PHREEQC The output file contains a lot of data and it is easy to get lost in all the numbers The results of the speciation calculations are reported under the heading Beginning of initial solution calculations Under Solution compo
79. pe Karoo mineral water 19 2 3 Saturation state s kes e ee ra pip ee 21 2 3 1 Calcite saturation state of seawater o oo 22 2 3 2 PHREEQC Exercise 2 Fluoride rich waters 23 Mineral dissolution and precipitation 27 3 1 PHT3D Exercise 1 Transport and mineral reactions 27 3 1 1 Spatial discretisation and flow problem 27 3 1 2 Running MODFLOW 29 3 1 3 PHT3D reaction definition 29 3 1 4 Initial and inflow concentrations 30 iii 3 1 5 Data input for solute transport 36 Running PHV3D io 23 442 6 s ebb Pe as 3 1 7 Visualization of simulation results 3 1 8 Comparison of simulation results 3 2 PHT3D Exercise 2 Precipitation dissolution fronts in acid mine drainage w sis 6 isora Sek ee a le gw ee ee aS 3 2 1 Spatial discretisation and flow problem 3 2 2 Running MODFLOW 3 2 3 Data input for solute transport properties 3 2 4 Data input for reactive transport 32 0 Running PHT3D eeren 404 4644 be 4 he es 3 2 6 Visualisation of results 04 3 2 7 Discussion of simulation results 3 2 8 Model variants 2 2 0008 3 3 PHT3D Exercise 3 Transport and nitrification of ion ex changeable ammonium 00 0000 8 0 3 3 1 Modification of the flow problem 3 3 2 Data input for solute transport properties 3 3 3 React
80. pter 4 1 9 C For pe which is not given in Table 1 use a value of 13 5 Attriute this water composition to the injection well under Spatial Attributes PHT3D PH gt ph 4 Source Sink and define a zone a single point in this instance at the well location Make sure that the reaction rate parameters 1 5 for the kineti cally controlled pyrite oxidation are set to values of 16 67 0 5 11 and 115 respectively Select the TVD scheme as advection package and define the dis persivities Run the reactive transport model Inspect the results How far did oxygen penetrate into the aquifer at the end of the simulation time 4 3 4 Seasonally changing redox zonation In the following part of the exercise we attempt to mimic the effects of a seasonally changing injection water temperature This will result in a dynamically changing redox zonation To simulate this effect proceed as follows 72 e To effectively consider the transiently changing water composi tions iPHT3D allows to upload ascii files that contain multiple water compositions For this exercise upload the pre prepared file that contains the time varying injectant water compositions To link the uploaded water compositions to the injection well modify the previous steady state input that was made for the Zone that is attributed to the injection well Go to Spatial Attributes gt PHT3D PH ph 4 Source Sink select the Zone and replace the
81. r composition at the start of the simulation needs to be entered under the Solutions Tab in the Background column and a second time in the Tab for Solu 2 This solution can also be used to define the recharge water composition The water composition of the contaminated water will be entered in the Tab for Solu 1 If the concentration of an aqueous component is 0 it is not necessary to enter the value Calcite is the only mineral that is considered in this simulation It must be activated and the corresponding initial concentrations must be entered under the Phase Tab The allocation of the appropriate water compoitions as model boundary condition can be done by selecting Spatial Attributes MT3DMS BTN gt btn 12 Boundary Condition and defin ing a zone line for the lower section of the inflow upstream boundary between y 0m and y 7m for which a value of 1 is allocated This defines that the concentrations at these grid cells will not change during the entire simulation The concentrations that iPHT3D will allocate to the model s grid cells are the ones defined as initial concentrations This means that the water com position that was defined as background water composition will at these locations also be used as inflow water composition To additionally define that the two different water compositions that sequentially enter the boundary s upper section we need to select Spatial Attributes PHT3D PH gt ph 4 Sou
82. rameters Model select the Model button and specify the attributes of the new model Select the values Xsection for Dimension Confined for Type and MODFLOW Family for Group Under Parameters gt Model select the Grid button and specify the details of the model domain and the spatial discretisation in the dialog that opens up Define a vertical transect that is 100m long and 10m thick The model domain may be discretised into grid cells of 4m length and 1m thickness Note that the model discretisation can be easily changed again later Select Parameters and select the Time button Set the Total Simula tion Time to 2000 days and select a time Step size of 20 days Go to Spatial Attributes gt MODFLOW LPF gt Ipf 8 to specify the hydraulic conductivity to 1 m day Go to Spatial Attributes gt MODFLOW LPF gt Ipf 9 to specify the vertical hydraulic conductivity to 1 m day Go to Spatial Attributes MT3DMS BTN gt btn 11 and specify the effective porosity to be 0 35 Implement hydraulic boundary conditions by defining the hydraulic heads at the upstream and the downstream boundary of the model First define the type of boundary condition by going to Spatial At tributes MODFLOW BAS6 gt bas 3 Boundary Conditions and setting the value to 1 active cell for all cells After that select Zone and use the Zone Tools to add one line at the inflow end of the model and one line at the effluent end In both
83. rce Sink and define a zone line for the uppermost 3 m y 7m to y 10m of the boundary In the dialog that comes up with the Zone selection enter two lines under Zone Value In the first line enter 0 1 and in a separate second line enter 1000 2 This means that Solution as defined in the Solutions Tab will be used until day 1000 and Solution 2 will be used until day 3000 To correctly define the recharge water composition select Spatial Attributes PHT3D PH gt ph 5 Recharge and define a zone line across the uppermost layer In the dialog that comes up 45 Chapter 3 Table 3 8 Concentrations of aqueous components and initial mineral con centration in PHT3D Exercise 3 Background and Contaminated water flushing water Chackgr C flush Ceont mol lt mol O 0 2 51 x 1074 0 Na 8 62 x 1074 1 30 x 1078 K 1 24 x 1074 1 30 x 1074 Ca 1 83 x 1073 1 50 x 1074 Mg 1 38 x 107 5 00 x 1075 Amm 0 6 87 x 1073 Cl 1 74 x 1078 3 23 x 107 S 6 9 89 x 10 4 1 56 x 1073 N 5 8 88 x 1074 0 N 3 0 N 0 0 0 C 4 2 82 x 1073 2 92 x 1078 C 4 0 0 pH 7 9 8 3 pe 13 5 0 Cinit mol z Calcite 0 1 with the Zone selection enter 2 under Zone Value to define that Solution_2 will be used for recharge The pollution source will be active during the first 1000 days while during the 2000 days following the pollution event the inflowing water is similar to the composition of the ambient water 3 3 4 Running PHT3D To run PHT3D pr
84. rid using the discretisation suggested in Table 4 5 e Select Parameters and select the Time button and set the Total Simulation Time to 600 days and set the Step Size to 5 days 60 Chemical kinetics Table 4 5 Flow and transport parameters used for PHT3D Exercise 4 Flow simulation steady state Total simulation time days 600 Time step size 5 d Model length m 150 Model thickness m 4 Aquifer top m 20 m Aquifer bottom m 16 m Grid spacing Ax m 5 Grid spacing Az m 0 20 Porosity 0 2 Horizontal hydraulic conductivity m day 20 Vertical hydraulic conductivity m day 20 Flux at upstream boundary m day 0 05 Groundwater recharge m day 0 001 Piezometric head downstream boundary m 19 m Longitudinal dispersivity m 0 1 Transversal dispersivity m 0 01 e Define a fixed head boundary condition at the downstream end Consider that such a fixed head boundary condition will not work for grid cells that are positioned above the water table Therefore make sure that the fixed head condition is only set for the depth zone below 19 m i e for locations below the water table position at the effluent end e Allocate suitable initial hydraulic heads e Enter a groundwater recharge rate of 0 001 m day as descried in the previous example e Define a specified flow boundary with a flow rate of 0 05 m day at the upstream boundary This can be implemented by adding a well The required injection rate equals the flow rate m
85. rm from C 4 CO3 released during dissolu tion and the F e 2 rich inflow solution At locations where calcite is com pletely dissolved siderite becomes the buffering mineral dissolving until it is also entirely removed Finally gibbsite Al OH 3 precipitation is the buffer 39 Chapter 3 Table 3 5 Aqueous concentrations used in the Walter et al Example Aqueous component Cini Ctailing mol I mol 1 pH 6 96 3 99 pe 1 67 7 69 C 4 3 94 x 1078 4 92 x 1074 S 6 7 48 x 1078 5 00 x 107 S 2 z Fe 2 5 39 x 1075 3 06 x 107 Fe 3 2 32 x 1078 1 99 x 1077 Mn 2 4 73 x 1075 9 83 x 1076 Ca 6 92 x 107 1 08 x 107 Mg 1 96 x 107 9 69 x 1074 Na 1 30 x 107 1 39 x 107 K 6 65 x 1075 7 93 x 10 4 Cl 1 03 x 107 1 19 x 1074 Al 1 27 x 1077 4 30 x 1073 Si 1 94 x 10 3 2 08 x 1078 Table 3 6 Mineral concentrations Walter et al Example Mineral Cinit mol 1 Calcite CaCOs3 1 95 x 107 Siderite FeCO3 4 22 x 107 Gibbsite Al OH s3 2 51 x 1078 amorphous Fe OH 1 86 x 1078 Gypsum CaS04 2H20 0 0 amorphous SiO2 4 07 x 107 mechanism at locations where both calcite and siderite are completely re moved The three different buffering mechanisms lead to three distinct levels of pH Also three distinct pe levels evolve with fronts from more reduced to more oxidised water The fronts have the same positions as those in the pH plot As the pe of the solution is completely controlled by the redox cou ple Fe 2 Fe 3 the pe chang
86. roundwater that is subsaturated i e STsyps lt 0 or precipitates from groundwater that is supersaturated i e STgyps gt 0 21 Chapter 2 100 Uncomplexed B Activity 80 60 4 40 4 of total concetration 20 4 i E Ca Mg Na K S04 HCO3 CO3 Species Figure 2 4 Effect of complexation and electrostatic shielding expressed as percentage of total concentration Based on calculations by 2 3 1 Calcite saturation state of seawater The basic principles of calculating a water s composition at equilibrium were presented in a landmark paper by and the theoretical framework presented therein continues to be applied in present day geochemical models Using a combination of equilibrium equations from the law of mass action and mass balance relations as explained in the previous section they calcu lated the distribution of dissolved species speciation in seawater at 25 C and 1 013 10 Pa 1 atm pressure Their results have been summarized in Figure 2 4 which shows the effects of both 1 the formation of aque ous complexes and 2 electrostatic shielding Sodium Na for example hardly forms any aqueous complex but due to electrostatic effects its activ ity amounts to only 75 of its total concentration More than 90 of the carbonate ion CO3_ is tied up in a complex and it can be seen that its activity amounts to a mere 2 of its total concentratio
87. rt to add the reactions by selecting Models MT3DMS gt Chem ical Reactions Click Edit and start editing the input data Under the Value Tab click Reset Matrix and select Linear equilibrium isotherm as sorption type Add a value of 3 x 1074 m3 kg for Kg which corresponds to 3 x 1077 L mg Click OK and leave the editor The value to enter for Kg obviously depends on the mass and the length units that are used for the model Those have to be consistent Therefore the value will depend on which units will be used for the definition of the bulk densitity With retardation being defined as R 1 2K 1 1 a bulk density of 1700 m3 kg will lead to a dimensionless retardation factor of 2 7 To enter the bulk density e Select Parameters Bulk density Cell by Cell Go to Value gt Reset Matrix and set the value to 1700 m kg Make sure that the value is applied to all three layers e Also you need to select the TVD option again under the MT3DMS SSolution MethodS menu 12 Single Species Transport e Moreover if you have changed dispersivity values make sure to change them back to 1 0 0 1 and 0 01 m for the longitudinal transverse and vertical dispersivities respectively e Run the MT3DMS simulation and check the concentration distribution in plan view and cross sections e Export the simulated concentrations for the breakthrough curve to an external ASCII file 1 2 2 First order Radioactive Decay or Biodegradation
88. ry species component e Go to Parameters Transport and select the Parameter button e In the dialog box that subsequently appears select ADV and the 3rd order TVD Scheme ULTIMATE as the Solution Scheme and click OK e Select Spatial Attributes gt MT3DMS DSP enter a longitudinal dispersivity value of 0 0067 m and click ok e Leave the values for the dispersivity ratios and diffusion coefficients as they are and click OK e Go to Spatial Attributes gt MT3DMS gt BTN gt btn 11 and specify the effective porosity to be 0 32 3 1 6 Running PHT3D To run PHT3D proceed as earlier in the exercise when executing MOD FLOW e Go to Parameters Chemistry and use the Pencil Button to write all the input files required to run PHT3D e Then go to Parameters Chemistry and use the Plume Button to start PHT3D 32 Mineral dissolution and precipitation 3 1 7 Visualization of simulation results To visualise the simulation results and to compare them with the results we need to extract the simulated concentrations in the form of concentration profiles for the end of the simulation time This can be achieved by creating a so called observation zone This option can be used to display the results graphically but also to export the results and save them in newly created data files e Select Spatial Attributes gt Observation OBS obs 1 gt Zone and create a zone as a horizontal line that maybe named linel
89. s Effective Porosity gt Reset Matrix and set the value to 0 3 Make sure that you apply these values to all three layers Leave the editor and save your changes Table 1 1 Flow and solute transport parameters used for MT3DMS Exercise i Parameter Value Flow simulation type steady state Simulation time days 3650 Model extent column direction m 250 Model extent row direction m 155 Model thickness m 6 Model top elevation m 6 Grid spacing Ax m 5 Grid spacing Ay m 5 Number of columns 50 Number of rows 31 Number of layers 3 Horizontal hydraulic conductivity m day 0 864 Vertical hydraulic conductivity m day 0 00864 Porosity 0 0 3 Piezometric head upstream boundary m 18 33 Piezometric head downstream boundary m 10 Longitudinal dispersivity az 1 0m Horizontal transverse dispersivity ary 0 1m Vertical transverse dispersivity ary 0 01 m Effective diffusion coefficient D 0 Chapter 1 T MT3DMS_EX1 PM5 Processing Modflow Pro Berk File Value Options Help Se LET 1 42E 01 1 53E 01 1 25E 01 1 17E 01 1 08E 01 1 33E 01 1 25E 01 1 17E 01 7 08E 01 1 25E 01 1 08E 01 217 0543 150 0089 5 1 1 44 2D Visualization Hydraulic head 18 33 Figure 1 2 Hydraulic head distribution under natural gradient conditions 1 1 4 Running MODFLOW You have now completed to set up the flow model and you are ready to run MODFLOW This s
90. s end Try if you can understand the syntax here In lines 1 2 and 4 of the toluene degradation rate block the constants from table 4 2 are specified Then the total concentrations of the aqueous components Toluene S 6 sulfate and sulfate reducing bacteria are read from PHREEQC s internal memory lines 10 22 and 30 Then in lines 32 and 34 a check is carried out to make sure that the concentrations of toluene and sulfate are sufficiently high Otherwise the rest of the code is skipped which means that the rate becomes zero Lines 40 50 and 80 implement the actual rate expression according to equation 4 5 The number of moles that react are calculated by multiplying the rate times the duration of the time step the TIME parameter in line 110 and passed back to PHREEQC in line 200 with the SAVE statement The rate expression for the decay of the sulfate reducing bacteria is constructed in a similar way under the identifier Sulfred Note that an extension of equation 4 6 is used in line 30 to limit the decay rate by multiplying with X Xinit X where Xini is the initial concentration of sulfate reducing bacteria The KINETICS keyword actually incorporated the rate expres sions into the PHREEQC simulation The names of the rate expres sions are listed Toluene_sulf and Sulfred and the identifier steps takes care of the temporal discretization The formula identifiers control the stoichiometric proportions of aqueous species wh
91. s mea sured concentrations of benzene toluene and sulphate for a vertical cross section along the contaminant plume The observed long thin plumes within the relatively homogeneous sand aquifer indicate low transversal dispersivities providing a significant challenge to achieve an accurate numerical description In the exercise you will prepare a simple vertical cross sectional model along a flow path of the contaminant plume within the uncon fined aquifer The main purpose of the exercise is to demonstrate the e incorporation of microbial kinetics into reactive transport models e finge controlled degradation of contaminant plumes e simulation of multi component NAPL Non Aqueous Phase Liq uid dissolution with the PHREEQC PHT3D framework and the effect of an ageing NAPL source e numerical problems that may be associated with advection dom inated reactive transport The exercise has been designed such that all the major phenom ena that characterise the behaviour of petroleum hydrocarbon plumes occur within the relatively short simulation period In reality NAPL concentrations near the water table might be higher and the mixture of compounds within the NAPL phase can also differ from the one used in the present example 4 2 1 Setting up the flow problem A cross sectional model will be set up Groundwater flow in the model domain is driven by both inflow from the upstream boundary and by groundwater recharge e Set up the model g
92. since they alter the waterSs composition as it travels through the subsurface Mineral dissolution and precipitation ion exchange transfer of solutes between solution and solid and redox processes reactions that involve electrons are the 3 main categories of chemical pro cesses that determine water quality Mixing of different water types which is more a physical than a chemical process further exerts great influence on the composition of groundwater Figure 1 illustrates the chemical processes that influence the concentra tion of the major ions of groundwater as it flows from the recharge area towards the sea Sulfate S027 derives from mineral sources such as pyrite and gypsum as well as from seawater and may be lost by conversion to HS7 e g by oxidation of organic matter Bicarbonate HCO3 is formed in the recharge area when COg is produced by aerobic respiration and decay of organic matter in the soil zone Its concentration is further influenced by mineral equilibria e g calcite dolomite siderite and redox reactions that involve sources of organic carbon in the aquifer itself Sodium Nat is contained in some silicate minerals that may dissolve but mainly originates Figure 1 Schematic representation of chemical processes that influence the concentration of major ions in coastal areas from seawater Cation exchange is the most important chemical process that affects its concentration Calcite
93. single value with 01 30 2 60 3 90 4 120 5 150 6 180 7 210 8 240 9 270 10 Chemical kinetics 300 11 330 12 Note that the heat transport approximation in this exercise is not very accurate as it ignores the retardation that occurs as a result of heat transfer between the injectant and the sediment matrix While neglected here it is possible to use MT3DMS PHT3D to exactly approximate the relevant heat transport equation 4 3 5 Mobilzation of arsenic In the last part or the exercise we are going to model the release of arsenic during oxidation of arsenopyrite and its subsequent sorption to ferrihydrite HfO In a first step we will model the dissolution of Arsenopyrite The rate expression that is used for Arsenopyrite will simply be linked to the rate of the computed pyrite oxidation HP TE TP PE Arsenopyrite Ha PT TE TP PP start 6 moles parm 1 get 1 200 save moles end The stoichiometric ratio at which arsenopyrite will be dissolved compared to pyrite is determined by parm 1 For example if parm 1 is set to 0 001 the dissolution of 1 mmol pyrite will be accompanied by the dissolution of 1 umol Arsenopyrite This is achieved by including put 1 in the rate expression for pyrite Corr eee renee ean one yrite TERT TTT TET TET ETE TE start 90 moles moles_oxy moles_nitr moles_0 100 if moles m then moles m 150 put moles 1 200 save moles 73 Chapter 4 end We do n
94. t a specific timestep and a species for which you want to see your results For example to see the pH contours after 1000 days select Results Aquifer Tstep and select 1000 and then go to Results gt Chemistry Species and select pH You can easily move forward and backward in time by clicking Tstep once more and then using your keyboard arrows to move forward or backward e Plot a breakthrough curve at a selected location to visualise how con centrations change over time First define the location for which you want to plot breakthrough curves by going to Spatial Attributes gt Observation OBS gt obs 1 Zone and use your mouse to position the cursor and clicking at the location where you wan to create an observation point In the dialog that comes up name the zone BTC1 No specific input is required for the Zone Value field The coordinates of the location selected by the mouse click can be modified if required 3 2 7 Discussion of simulation results Only a few aspects of the geochemical evolution observed in the model sim ulation are described in the following For a detailed description see The acidic inflow solution is initially buffered by calcite CaCO3 main taining the pH at an almost neutral level 6 5 7 In this zone gypsum CaSO4 2H2O is formed from calcium released during calcite dissolution and sulphate 5 6 from the inflow solution Cin 5 6 gt Cbackgr s 6 Simi larly siderite F eCO3 can fo
95. tep will provide MT3DMS with the flow field that is used to compute advective dispersive transport e Select Models MODFLOW gt Run e A window appears that lists all the files that will be generated There are also a number of options which you can leave as they are Click OK to start the calculations First PMWIN will generate the nec essary MODFLOW files and then it will call MODFLOW to run the simulation model e When the calculation has finished hit any key to close the command window and return to PMWIN 1 1 5 Visualising hydraulic heads e Visualise results by selecting Tools 2D Visualisation On the MOD FLOW tab select Hydraulic head and klick OK You should get a steady state head distribution like shown Figure 1 2 Single Species Transport 1 1 6 Setting up the transport model Set the boundary conditions for concentrations by selecting Select Grid Cell Status ICBUND Trabsport models Leave the editor with out making changes The default value for ICBUND of 1 active cell will remain at all grid cells Select Models gt MT3DMS Simulation Settings Define a new species by adding a name like Tracer to the Description list Make it active by placing a tick mark in the box behind it and then click OK Leave the default Type of Reaction No kinetic reaction is simulated Define the initial concentrations by selecting Models gt MT3DMS gt Source Sink Concentration Constant Head
96. ters the model domain through the fixed head boundary discharges to the well during active pumping However some groundwater infiltrates to the river On the other hand some river water is also drawn towards the pumping well The model setup is shown in Figure 4 2 West River TCE Toluene Clay Figure 4 2 PHT3D Exercise 6 Model Setup 4 4 2 NAPL source zones The contamination sources are modelled by defining two discrete NAPL contaminated zones On the Western side of the river a PCE source is defined by allocating an initial concentration of 1 mol l of Pcenapl to a discrete source zone situated on top of the clay lense Similarly the 75 Chapter 4 Eastern source is defined by allocating 1 mol l of Tcenapl and 1 mol l of Tolunapl to a discrete zone The rate expression that controls the kinetically controlled release of PCE is it Pcenapl start 10 mPcenapl tot Pcenapl 15 if mPcenapl lt 1e 10 then goto 200 20 solub_Pce 1 28 131 39 aqueous solubility of PCE 22 frac_of_solub parm 2 Fraction of single species sol ubility 25 mPce tot Pce_1 tot Pce_h 50 msolub_Pce frac_of_solub solub_Pce 60 rate parm 1 msolub_Pce mPce 65 if mPce gt msolub_Pce then rate 0 70 moles rate time 80 if moles gt m then moles m 200 save moles end In the reaction database similar rate expressions are defined for TCE and toluene 4 4 3 Formati
97. ting the thought process along a path that may result in acceptance rejection or modification of the original Introduction hypothesis Geochemical models were originally developed in the 1960 s to calcu late the speciation dissolved ions The original models were soon improved and extended to include chemical reactions that alter the water composi tion such as mineral equilibria and cation exchange By virtue of the increase in computer computational power it became possible to couple geo chemical models with hydrological models to calculate how the water com position changes as it travels through the subsurface Today s models have reached a level of sophistication that allows us to simulate real world processes to understand and explain field observations Elaborate treatment of the basic concepts and applications of aqueous geochemical modeling can be found in for example and This course manual provides a treatment of the basic principles of the models as well as exercises for MT3DMS PHREEQC 2 and PHT3D Extensive use is made of the following packages e PHREEQC for Windows a graphical user interface for PHREEQC which has the advantage that it s not necessary to switch between programs to create input perform calculations or view output The program can be installed on any computer that has Windows95 or higher It can be obtained for free from http pfw antipodes nl download html e PMWIN
98. tio between hand and computer calculated concentra tions change How much fluorite precipitates dissolves in the 1 1 mixture of samples 1 and 3 25 Chapter 3 Mineral dissolution and precipitation 3 1 PHT3D Exercise 1 Transport and mineral re actions The present exercise is the first PHT3D application example Building on the experience gained from the first PHREEQC examples it combines a simple one dimensional flow mass transport simulation with mineral pre cipitation dissolution reactions The case simulated in this exercise was originally presented by for a model verification of their MST1D code against the CHEMTRNS model by It involves a one dimensional model domain in which an aqueous water composition that is in equilibrium with two minerals calcite and dolomite is successively replaced i e flushed by water of a different chemical compo sition leading to multiple precipitation dissolution fronts Dolomite is not present initially but is formed temporally 3 1 1 Spatial discretisation and flow problem In order to exactly reproduce the discretisation chosen by the total length of the model domain must be set to 0 5 m and subdivided into 50 grid cells of 0 01 m length Each of these cells has a width of 1 m and also a height of 1 m height With this discretisation the model has 50 columns 1 row and 1 layer The total simulation time is 0 24305 days The temporal discretisation time step length is set to 0 01
99. ular three dimensional multi species transport model for simulation of advection disper sion and chemical reactions of contaminants in groundwater systems documentation and user s guide Tech rep U S Army Engineer Research and Development Center Contract Report SERDP 99 1 Vicksburg MS 84 Acknowledgments Acknowledgments Thanks to Bruno Haerens Janek Greskowiak Phil Ham Christopher Schmidt Boris van Breukelen and many others for their various con tributions and for testing the modelling examples 85 Appendix I Options for entering TIC alkalinity and pH in PHREEQC There are a number of options available in PHREEQC to enter the carbon concentration alkalinity and the pH The best option depends on the problem at hand The options are 1 Enter pH and alkalinity Total moles of carbon 4 is calculated by the program 2 Enter pH and total carbon Alkalinity is calculated by the pro gram 3 Enter alkalinity and total carbon pH is calculated by the pro gram 4 Enter pH estimate of total carbon and a partial pressure of COs Program will adjust total carbon until desired Poo is attained 5 Enter total alkalinity or total carbon estimate of pH and a partial pressure of CO If possible the program will adjust the pH to attain the desired Peo The first option is used when both pH and alkalinity were deter mined in the field The second option should be used when alkalinity was not titrated in th
100. ultiplied by the height of the aquifer e to make sure that the recharge is always applied to the highest ac tive grid cell Go to Parameters Flow and select the Parameter button e In the dialog box that subsequently appears select RCH and se lect under rch 1 Flags for recharge the highest active option 61 Chapter 4 e To prevent drying grid cells activate MODFLOWS re wetting ca pabilities use the same dialog and go to LPF and select under Ipf 6 wetting active or inactive the option active under lpf 7 calculation of wetting set the wetting factor WETFCT to a value of 0 1 and select the option h BOT WETFCT x THRESH Set the wetting treshold TRESH to a value of 01 e Now rerun MODFLOW to regenerate the flow file mt3d flo which will later be used by the transport model e As acheck your simulated head at the upstream boundary should be 19 56 m If your simulated head does not agree with this value revisit the previous steps and try to debug the problem 4 2 2 Simulating nonreactive transport Before the reactive transport simulations are started we will briefly check the setup of the model with a tracer transport simulation Per form the following steps e Define the settings for advective and dispersive transport e Set the initial concentration background to 0 e Decide what type of boundary conditions would suit the problem e Define a contaminant source by defining a concentration value of 1 f
101. ut file contains BASIC statements that control the output to a text file in spreadsheet format The distribution coefficient of Zn is calculated for both the weak and the strong sites as well as an overall distribution coefficient The file is called kd_surf prn and can be opened in the Grid tab of PRHEEQC for Windows after the calculation has finished What is the value of the overall distribution coefficient for Zn Answer Kg Test the effect on the Ky of e doubling the Zn concentration Answer Kg e doubling the Mg concentration Answer Kg e doubling the amount of ferrihydrite Answer Kg e changing pH from 7 to 5 Answer Kg 50 Chapter 4 Modelling kinetically controlled reactions In the preceding chapters we have only considered equilibrium chem istry i e situations where chemical reactions occur instantaneously For reactive transport problems the simplifying assumption of instan taneous reactions is warranted as long as the reaction time scale is fast compared to the transport time scale which is known as the local equilibrium assumption LEA e g For many applications however the kinetics of chemical reactions have to be considered Examples include the slow dissolution of silli cate minerals microbially mediated decay of organic pollutants and the oxidation of reduced mineral substances The versatility of PHREEQC allows for simultaneous modelling of both equilibrium
102. ut file is given below Only the pH 6 8 has not been entered yet Look up the identifier to specify pH and enter it in the input file SOLUTION 1 Cape Karoo still units mg l enter pH on this line Na 29 K 0 9 Mg 2 2 Ca 34 Cl 46 S 6 6 8 N 0 29 F 0 72 Fe 0 08 Alkalinity 83 as CaCO3 END The concentrations here are entered as mg l Because PHREEQC recal culates all concentrations to moles it is important to specify what molecular weights must be used for recalulating mass based units such as mg l In this case alkalinity was reported as CaCO3 and the statement as CaCO3 after the concentration ensures that the recaculation is done using the molec ular weight of CaCO3 Redox sensitive species can have different valence states Elemental sul fur in sulfate S037 has a valence of 6 and therefore the element name S 6 indicates the sulfate ion The same holds for nitrogen and iron but be cause the label does not specify what form they are in we omit the valence state If we knew however that nitrogen for example was present in the 20 Specation and saturation indices form of NO3 we would have entered N 5 because in NO nitrogen has a charge of 5 Another special case is HCO which can not be entered directly but is specified here as Alkalinity note that element names are case sensitive In natural waters were the alkalinity comes from carbonate species this is a valid assumption The differ
103. utton iPHT3D will then read and interpret the PHT3D database file In this next step we define all initial concentrations of all aqueous components and minerals and also the concentrations at the upstream model boundary which will be allocated to the inflowing groundwa ter The composition of the ambient water initial water composition at the start of the simulation needs to be entered under the Solutions Tab in the Background column Activate all relevant species by ticking the appropriate boxes for pH pe C 4 S 6 S 2 Fe 2 Fe 3 Mn 2 Ca Mg Na Cl Al and Si and enter the concentrations The set of minerals that are considered in this simulation must be activated and the corresponding mineral composition must be entered under the Phase Tab In addition to the background water composition we will also need to define the composition of the acidic tailings water that enters the model domain in the upper section of the upstream boundary This water composition will be defined in the Solution 1 column under the Solutions Tab For the lower section of the inflow boundary it is as sumed that the inflowing water is the same as the background water composition The allocation of the appropriate water compositions as model bound ary condition can be done by selecting Spatial Attributes gt MT3DMS BTN gt btn 12 Boundary Condition and defining a zone line for the lower section of the inflow upstream boundary between y

Download Pdf Manuals

Related Search

Here here we go again lyrics hereditary heretic heredity here movie here comes the sun heretic definition hereditary meaning hereinafter heresy definition here\u0027s johnny hereditary movie hereby here comes the guide hereditary angioedema hereditary hemochromatosis here comes the sun lyrics heretic movie hereditary spherocytosis hereditary hemorrhagic telangiectasia herencia here to slay heredia costa rica heretic game heretic streaming

Related Contents

HDMI DVD PLAYER - Craig Electronics Inc.  User Manual - e-catalog  E KERN 572/573/KB/DS/FKB/FCB/KBJ  Severin KA4767 coffee maker  Mode d`emploi Protectiv Plus:Mode d`emploi  Verbatim Portable Hard Drive Combo USB User's Manual  iPad サポータ養成講座  V7 Vertical Flip Case for iPhone® 6 Plus - White  manual de instruções no break senoidal kns  "service manual"  

Copyright © All rights reserved. Failed to retrieve file