Home

User Manual of the Multicomponent Variably - PC

1. 1 ecce esee e eee eren ee eee eene eee teen nest ene sesta 7 2 2 2 Initial and boundary conditions for the solute transport equation 9 2 3 Heat transport in the vadose zone eere eee eee sees sees seen esent tn sten sten n tena enn 9 2 3 1 The heat transport equation e eee e eere nue seen en nee enar ra a a eaae e een asse eae 9 2 3 2 Initial and boundary conditions for heat transport 10 ZA Geochemical redactions aie ROREM A IR RR DO da OTR OL 10 2 4 1 Homogeneous aqueous reactions c ee ecce eere ee eene enne ee ee ene sesta ses een ae 11 2 4 Heterogeneous ion exchange reactions ee eee ceres eerte ee eren enean etn 11 2 4 3 Heterogeneous mineral dissolution precipitation 12 LAA Kinetic reactions uie dre eee obi eissoes sosoca EBERT EUR VEN U R RR BSc PERRO sesse BERE 13 2 5 Multicomponent reactive transport e sessi esee eese eese eene eene eene ee tone tn setate setae s een a 13 2 6 Coupling procedures sso he x te XR nR VE REPRE RE VY REX YR S GRE Re re eR gut 15 Del Model SIPUCDEEO o dr Doi ape I ER a AE Un del PO ad Qus Aa caa dri e ERREUR 18 2 7 1 CPHREEGQC 12er d epe DDR ENG Ie Rp M obe EU MN EE 18 2 52 JHYDRUSSID 4 docete viu dietro soos oidor oe isese risoo eso odote o DSSS ssion 19 3 Description of data input eessoesscocsssccssecesc
2. 8 Sor oO 9 EIE GT CE 10 gt Sem 4 hn Cont Someone inp Freundlich term atat no check mole balance SorCont 12 Log k AE sie m Joe unes S 13 SURFACE 1 14 equilibrate 1 15 no edl true iG Sos S 1 8 S Total amount of adsorption sites Cont its aqueous speciation and the surface species Sor are defined on lines 1 6 with the PHREEQC keywords SOLUTION MASTER SPECIES SOLUTION SPECIES and SURFACE MASTER SPECIES The definition of the mass action equation of the Freundlich isotherm Eq 4 4b is given on lines 10 12 with the PHREEQC keyword SURFACE SPECIES Two options are included 1 no check since the mole balance of the reaction equation line 10 is not fulfilled due to nonlinearity of the Freundlich isotherm nrzl and 2 mole balance to ensure balancing of 1 Sor and 1 Cont The association coefficient K is defined by K Sor which is approximately equal to the total amount of adsorption sites S if S is chosen very large e g 10 mol I put S 10 in line 16 As an example for K7 5 p 1 5 and S 10 we obtained a K value log_k of 99 1249 Verification Problem 3 Steady state transport of nonlinearly adsorbed contaminant STADS In this problem we consider saturated steady state water flow and single component Cont transport through a soil column of 1 m length for a period of 1000 d Very low initial concentrations of Cont are assumed to be present 10 mol I in all cases Transpor
3. Calculate new water flow properties Compute new water flow state variables p T7 H Calculate new heat transport properties Compute new heat transport state variables p H T based on H Transport step Compute new chemical state variables rL no chemical reactions Ctransport based on H Calculate new solute transport properties pt H Geochemical step Compute new geochemical state variables no transport cat Gt based on Ctransport Calculate new geochemical properties pere 1 End of simulation Figure 2 1 Schematic of the modelling approach of the coupled HP1 model 18 2 7 Model structure In this section we summarize the changes implemented in both codes While the Fortran routines of HYDRUS 1D are compiled into Hydrus dll the c functions of PHREEQC are compiled into PHREEQC dll 2 7 1 PHREEQC The original PHREEQC code contained the following files Advect c Basic c Cll c Integr c Inverse c Kinetics c Main c Mainsubs c Model c P2clib c Parse c Phgalloc c Prep c Print c Read c Readtr c Spread c Step c Structur c Tidy c Transp c and Utility c In the coupled HP1 we modified mostly only the main c routine and added a new file Hydr tr c main c From the main routine main Phreeqc we call first the read text file function that reads from the Species in file the names of components to be transported by HYDRUS 1D and stores t
4. sess 51 Table 4 11 Overview of aqueous equilibrium reactions and cation exchange half reactions with corresponding equilibrium constant esee 52 Table 4 12 Chemical composition of the initial pore water and the inflowing alkaline SODITON cae E aa ae ec e E ea A eL EE 53 Table 5 1 Soil hydraulic and other properties of six soil horizons sss 59 Table 5 2 Chemical composition of the rain water from Stolk 2001 and in the initial CETT S TE 60 vi vii LIST OF FIGURES Figure 2 1 Schematic of the modelling approach of the coupled HP1 model 17 Figure 4 1 Outflow concentrations for Verification Problem 1 Full and dashed lines are results for physical equilibrium and nonequilibrium transport respectively obtained with HYDRUS 1D Dots are results obtained with 1 l M a CS Ne ae Oe Sd Ml neat ee 26 Figure 4 2 Outflow concentrations for Verification Problem 2 The full line was generated with HYDRUS ID while dots were obtained with HP 27 Figure 4 3 Depth profiles of Cont after 250 500 and 750 d for different simulations defined in Table 4 1 for Verification Problem 3 Tests consider effects of grid size for PHREEQC a and effects of maximum time steps for HP1 b 31 Figure 4 4 Depth profiles of Cont after 250 500 and 750 d for different simulations defined in Table 4 1 Verification Problem 4 with a first order decay coefficient of O 2d
5. OPEN REPORT SCK CEN BLG 998 05 DJa P 36 User Manual of the Multicomponent Variably Saturated Flow and Transport Model HP1 Description Verification and Examples Version 1 0 Diederik Jacques and Jirka im nek Waste and Disposal SCK CEN Mol Belgium Department of Environmental Sciences University of California Riverside USA June 2005 SCK CEN Waste amp Disposal Department Boeretang 200 2400 Mol Belgium OPEN REPORT OF THE BELGIAN NUCLEAR RESEARCH CENTRE SCK CEN BLG 998 User Manual of the Multicomponent Variably Saturated Flow and Transport Model HP1 Description Verification and Examples Version 1 0 Diederik Jacques and Jirka im nek Waste and Disposal SCK CEN Mol Belgium Department of Environmental Sciences University of California Riverside USA June 2005 Status Unclassified ISSN 1379 2407 SCK CEN Waste amp Disposal Department Boeretang 200 2400 Mol Belgium SCK CEN Belgian Nuclear Research Centre Boeretang 200 2400 Mol Belgium Phone 32 14 33 21 11 Fax 32 14315021 http www sckcen be Contact Knowledge Centre library sckcen be RESTRICTED All property rights and copyright are reserved Any communication or reproduction of this document and any communication or use of its content without explicit authorization is prohibited Any infringement to this rule is illegal and entitles to claim damages from the infringer withou
6. Adler 2001 p 55 60 In this report the same model is used in both CRUNCH and HP1 51 Table 4 9 Mineralogical reactions Log K molar volume volume percent and moles present in Opalinus Clay and secondary minerals considered mineral Reaction LogK Molar Volume Moles volume 99 Primary minerals Kaolinite ALSiO4OH 6H 2AI 2SiO 5H O 6 36 99 520 23 2 31 Illite Ko4Mgo5sAL 3Sis Oi OH 8 H 0 6 K 0 25 Mg 2 3 AI 3 5 SiO 5 HO 8 51 z 19 0 38 Quartz SiO SiO 3 91 22 680 11 5 5 06 Calcite CaCO3 Ca COT 8 51 36 934 13 5 3 655 Dolomite CaMg CO Ca Mg 2 COT 16 71 64 365 2 0 0 31 Gypsum CaSO 2H 0 Ca SO 2 HO 4 48 74 69 0 02 0 00268 Secondary minerals Hydrotalci Mg4AlO OH 14 H 4Mg 2 Al 10 7544 207 57 0 00 0 0 te H O Mg Si6O1s 0H 2 6H O 8 H 4 Mg 6SiO 29 96 285 60 0 00 0 00 Sepiolite 11 H0 Used as input in CRUNCH Used as input in HP1 a lower reactive surface area during dissolution Mmo and Mm are the initial and current numbers of moles in HP1 and volume percentages in CRUNCH of mineral m n is the reaction order of the OH dependency of the reaction rate ko is the intrinsic rate constant of mineral m mol m sec K is the equilibrium constant of the reaction and Q is the ion activity product Table 4 10 lists the parameters of Eq 4 7 Due to the stiffness of the system a combination of high and low rate constants not all dissol
7. Rain water 5 25 69 32 64 4 6 8 0 0 0 Initial soil solution A 3 4 69 727 64 4 97 8 0 80 50 2 5 E 3 5 69 521 64 4 65 8 0 38 24 1 2 Bhl 3 6 69 398 64 4 39 8 0 43 23 1 1 Bh2 3 8 69 287 64 4 33 8 0 41 21 1 0 BC 4 4 69 235 64 4 55 8 0 63 33 1 6 Cl 4 4 69 252 64 4 76 8 0 33 21 1 0 C2 4 5 69 131 64 4 27 8 0 20 14 0 7 Br is used as charge balance 61 The initial water content was estimated by assuming a constant flux of 0 1206 cm day corresponding to the long term actual infiltration rate Initial concentrations of Na K Mg were assumed to be equal to the concentration in the precipitation while initial concentrations of Ca Cd and Zn were obtained from Seuntjens 2000 Table 5 2 Pb concentrations were arbitrarily set to be 20 times smaller than the Zn concentrations no data available The simulation was carried out using the following numerical guidelines a spatial discretization of 1 cm a performance index of 0 1 and a maximum time step of 0 125 5 1 2 Analysis of the simulation results 5 1 2 1 Accuracy analysis To assess the accuracy of the simulation the mass balance was evaluated in terms of the following two variables as a function of time 91 9 att fata 5 1 N cells ud Q E Xe xn Je Q M Yee os 5 2 i l where j 1 Nm Oln and QZ are the cumulative inflow and outflow rates of the jth master species moles L respectively q and iy are solute flux boundary conditions at the bott
8. and the initial and boundary conditions DISCRETIZATION INITIAL CONDITIONS and BOUNDARY CONDITIONS the flow and transport properties FLOW and TRANSPORT and the output OUTPUT The keyword POROSITY defines the porosity which is fixed in this simulation Numerical information for solving the equations are defined in RUNTIME 4 3 2 3 Comparison between CRUNCH GIMRT CRUNCH OS and HP1 Simulations were performed with two different codes HP1 and CRUNCH where CRUNCH was applied in two different modes the global implicit option in the CRUNCH model CRUNCH GIMRT and the sequential non iterative option in the CRUNCH model CRUNCH OS The coupled HP1 model runs in a sequential non iterative mode For the latter model different maximum time steps corresponding to different Courant numbers 0 125 and 0 065 were used The maximum allowed Courant number of the CRUNCH OS simulation was set to 0 125 Figure 4 11 shows distributions of pH and selected elements after 0 3 0 7 and 1 1 y Outflow curves are shown in Figure 4 12 In general correspondence between the different modelling results is relatively good Very good agreement was obtained for the propagation of both high pH and elements that are only influenced by cation exchange Na and K the latter not shown Figure 4 11a b and Figure 4 12b The pH near the inlet as obtained with HP1 was slightly larger than the pH of the CRUNCH simulations The positions of the propagation waves of Si
9. index number of secondary species index number of secondary exchange species index number of master species index number of master exchanger index number of solute in a decay chain intrinsic rate constant of mineral m unsaturated hydraulic conductivity empirical adsorption coefficient 75 M L7 K L L ML T K equivalents LTY LTY L L L L ML T L T7 M L Ms 76 Kj adsorption coefficient in mass per unit of volume water MC 9T Kj equilibrium constant of reaction for ith secondary exchange species K equilibrium constant of reaction for ith aqueous secondary species K equilibrium constant of reaction for ith mineral K saturated hydraulic conductivity L T m parameter in the water retention curve m index number of a mineral M amount of the jth master species in the flow region at time t M L Mi amount of the jth master species in element i at time t M L Mi amount of the jth master species in the flow region at the beginning of the M L simulation Mi amount of the jth master species in element i at the beginning of the simulation M L t pore connectivity parameter n parameter in the water retention curve n reaction order of the OH dependence of the reaction rate np empirical Freundlich exponent Nije moles of ith secondary exchange species on j th master exchan
10. parm 1 TOT water MOL Cont 6 30 moles rate TIME 7 40 SAVE moles 8 end 9 KINETICS 1 10 degradation ii wexinituke Come ip 90s 2 o ASE 33 product of the stoichiometric element coefficient i e 1 in the formula option multiplied by the moles of reaction during a particular time step is positive or consumed i e when the product is negative during the kinetic reaction Since in this example the coefficient for the element Cont is 1 line 11 and the reaction progress is negative the concentration of Cont will decrease with the formation of the imaginary phase degradation Note that we could write the input also with a negative coefficient and a positive reaction progress i e dissolution of the phase However in that case we could reach complete dissolution of the phase and consequently termination of the decay reactions The last option under KINETICS in Box 4 2 is parms for the purpose of defining parameters in the rate expression Verification Problem 4 Steady state transport of nonlinearly adsorbing contaminant with first order decay STDECAY This verification problem is the same as the previous one i e Verification Problem 3 but with the inclusion of first order decay All transport and simulation parameters are equal to those of Verification Problem 3 see Table 4 1 However we did not perform the P 3 and HP 1 simulations because the former P 3 was computationally inefficient and the latt
11. 0 0 02 0 04 0 06 0 08 Depth in clay core m Depth in clay core m CRUNCH GIMRT _ HP1 Cr 0 125 CRUNCH OS HP1 Cr 0 065 Figure 4 11 Distribution of pH and selected elements versus depth mol l after 0 3 0 7 and 1 1 y of infiltration of a high pH solution in an Opalinus Clay core 57 0 0006 A 0 0004 0 0002 0 02 04 06 08 1 0 02 04 06 08 1 Time y Time y 0 04 11 e f 0 03 10 I 2 0 02 L9 0 01 8 0 0 02 04 06 08 1 0 02 04 06 08 1 Time y Time y CRUNCH GIMRT HP1 Cr 0 125 CHUNG OS sssr HP1 Cr 0 065 Figure 4 12 Outflow curves of selected elements mol l during infiltration of a high pH solution in an Opalinus Clay core 58 59 5 Example This section considers an example that combines transient variably saturated water flow solute transport and geochemical reactions The example is used to illustrate some of the capabilities of HP1 5 1 Long term transient flow and transport of major cations and heavy metals in a soil profile HEAVMET 5 1 1 Problem definition and governing chemical reactions This example considers similar processes as those discussed in Verification Problem 8 The problem involves the transport of major cations Na K Ca and Mg and heavy metals Cd Zn and Pb in a l m deep multi layered soil profile subject to atmospheric boundary conditions While steady state water flow was considered in Verification Problem
12. 0 37 and 0 15 as a function of depth The bottom boundary condition for water flow is free drainage HYDRUS 1D was used to calculate the steady state water content profile corresponding to these boundary conditions The dispersivity and diffusion coefficient were taken to be 0 05 m and 9 2 10 m s respectively Table 4 7 Overview of aqueous equilibrium reactions and corresponding equilibrium constants data from phreeqc dat database Parkhurst and Appelo 1999 Nr Aqueous speciation reaction Log K 1 H O OH H 14 2 Na H O NaOH H 14 18 3 K HO KOH H 14 46 4 Ca HO CaOH H 12 78 al O ooo Mg t MO al io MIRRORS NER Cd Pb Zn 6 X HO XOH H 10 08 7 71 8 96 7 X 2 H O X OH 2 H 20 35 17 12 16 90 8 X 3 H O X OH 3 3H 33 30 28 00 2840 9 x 4 H O X OH 4 H 47 35 39 70 41 20 10 X CT XC 1 98 1 60 0 43 11 xX 2 Cr XCl 2 60 1 80 0 45 12 X 3 CT XC 2 40 1 70 0 5 13 X 4Cl xcly 1 38 0 2 46 Table 4 8 Log K parameters for multi site exchange complex Y NaY KY MgY CaY CdY gt PbY gt ZnY gt exchanger 1 0 0 3 0 4 0 2 0 2 0 05 0 2 HY HYa HYb HYc HYd HYe HYf 1 65 3 3 4 95 6 85 9 6 12 35 The value for NaY is taken from Appelo et al 1998 Values for the other complexes are taken from the phreeqc dat database Parkhurst and Appelo 1999 and adapted relative to the K for NaY Values taken from Appel
13. 1 3 72 1 0 10 C2 75 100 1 56 0 08 0 021 0 39 2 1 4 33 1 0 10 60 Cumulative P E P E cm 1972 1974 1976 1978 1980 1982 Time year Precipitation P Sa Evaporation E SS SS Potential surface flux P E e Actual surface flux Figure 5 1 Cumulative precipitation evaporation and surface flux into the soil between 1972 and 1982 Figure 5 1 shows precipitation P potential evaporation E and potential flow across the upper boundary P E Since the amount of water available in the soil profile for evaporation was smaller than that required for potential evaporation the actual evaporation rate will be smaller than the potential evaporation rate Consequently the actual surface flux was larger than the potential surface flux see Figure 5 1 For the considered 10 year period the actual surface flux around 45 cm year was somewhat larger than expected for this region of Belgium This is because no water uptake by plant roots was considered in this example Element concentrations in the precipitation incoming water were obtained from Stolk 2001 for station 231 located in Gilze Rijen the Netherlands close to the Belgian border The composition was based on the average of 13 measurements during 1999 Table A5 in Stolk 2001 Table 5 2 Table 5 2 Chemical composition of the rain water from Stolk 2001 and in the initial soil solution 3 Concentration pH Cl Br Na K Ca Mg Cd Zn Pb umol T
14. Al and Mg Figure 4 11c e and Ca similar to Mg not shown were identical for the different models although a slight difference occurred in the concentration of these elements The reason for this is not clear it might be due to difference in the precipitation kinetics of the secondary minerals hydrotalcite and sepiolite Precipitation of these minerals in HP1 is governed by thermodynamic equilibrium whereas kinetic precipitation is assumed in the CRUNCH simulations Due to this equilibrium assumption Mg concentrations as calculated with HP1 should be slightly smaller than those obtained with CRUNCH The Si or Al concentrations are consequently also slightly larger Other minor differences between the two models e g the exact activity coefficients may have contributed also to this difference 55 Compared to cation exchange problems see section 4 3 1 problems involving kinetic dissolution precipitation of minerals generally require somewhat smaller time steps if a sequential non iterative coupling method is used Nevertheless acceptable results were obtained compared to those using the global implicit method 56 0 25 0 2 0 15 ss z 0 1 0 05 IO 0 0 02 0 04 0 06 0 08 0 0 02 0 04 0 06 0 08 Depth in clay core m Depth in clay core m 0 016 0 012 0 008 0 004 0 0 0 02 0 04 0 06 0 08 0 0 02 0 04 0 06 0 08 Depth in clay core m Depth in clay core m 8x10 6x10 0 4x10 2x10 0x10 0 0 02 0 04 0 06 0 08
15. IRAE TRA 49 Figure 4 11 Distribution of pH and selected elements versus depth mol l after 0 3 0 7 and 1 1 y of infiltration of a high pH solution in an Opalinus Clay COV Oe PINK TUN E 56 Figure 4 12 Outflow curves of selected elements mol l during infiltration of a high p H solution in an Opalinus Clay COFe i ee o e oed Re eb read delto 57 Figure 5 1 Cumulative precipitation evaporation and surface flux into the soil between 19 72 and FUGA iaaeaie e b i ta asd san D etiani s PE TEE EEES Eosin Enea 60 Figure 5 2 Time series of inflow minus ouftlow dashed line and change in the total amount present in the soil profile dots for Cl Na Ca and Cd 62 viii Figure 5 3 Relative mass balance errors amp as a function of time for Na Ca and Cd 63 Figure 5 4 Space time plots of a water content b amount of Cl mmol 1000 cm soil c pH and d amount of aqueous Cd mmol 1000 cm soil in the A horizon qe Zem depth osi despues diad est rta dete van Bea cated RRNGU ee d qus 64 Figure 5 5 Time series of water content Cl concentration mmol 1000 cm soil pH and Cd concentration mmol 1000 cm soil in the third horizon Bh1 horizon t 2 DCMU ACTN ase a e EU ei esu ie esu TI d n unt 65 1 Introduction 1 1 Background The migration of many naturally occurring elements and contaminants in the subsurface is affected by a multitude of complex interactive physical chemical
16. SO4 4 9 11 Ca H O CaOH H 12 78 12 Ca HCO 7 CaCO H 7 0017 13 Ca HCO CaHCO 1 0467 14 Ca CT CaCl 0 6956 15 Ca 2 CI CaCl 0 6436 16 Ca SO CaSO 2 111 17 Cl H HCI 0 67 18 K H O KOH H 14 46 19 K CI KCl 1 4946 20 K SO KSO 0 8796 21 K H SO KHSO 0 8136 22 Mg H O MgOH H 11 44 23 Na H O NaOH H 14 18 24 Na HCO NaHCO 0 1541 25 Na HCO NaCO H 9 8144 26 Na CI NaCl 0 777 27 Na SO NaSO 0 82 28 SO H HSO 1 9794 29 SOF 2 H H SO 1 0209 30 SiO 2 H O H5SiO 2H 22 96 31 SiO H O HSiO H 9 9525 32 SiO Na H O NaHSiO H 8 3040 Cation exchange half reactions X is the exchange site 33 Na X NaX 0 00 34 K X KX 0 757 35 Ca 2X CaX 0 604 36 Mg 2 X MgX 0 505 An overview of the considered aqueous equilibrium reactions and the ion exchange reactions is given in Table 4 11 The total cation exchange capacity of Opalinus Clay is 120 meq per kg of rock The initial concentrations and concentrations in the infiltrating water are given in Table 4 12 4 3 2 2 Database and model input HP1 The chemical reactions and their equilibrium constants are stored in the database phreeqc dat This problem uses the following keywords SOLUTION 53 Table 4 12 Chemical composition of the initial pore water and the inflowing alkaline
17. c o o O O 0 02 0 0 20 40 60 80 100 0 20 40 80 100 Depth cm Depth cm 0 1 0 1 0 08 2 Contc 0 08 Contc o o E E c 0 06 c 0 06 2 o 0 044 S 0 04 o o c c o Oo O O 0 02 0 02 0 0 0 20 40 60 80 100 00 Depth cm Depth an HYDRUS 1D HP1 1 At 30000 s 2 left HP1 3 At 86400 s o 0 4 right e HP1 2 At 10000 s o 2 left HP1 4 Atsa 86400 s 0 2 right Figure 4 5 Distribution of Conta Contb and Contc versus depth after 250 500 and 1000 d for different simulations as defined in Table 4 3 Verification Problem 5 39 4 2 Multicomponent transport during steady state flow In this section we test different keywords of the PHREEQC program while using the coupled HP1 model The following specific reactions are tested 1 cation exchange EXCHANGE and 2 equilibrium and kinetic dissolution precipitation of minerals EQUILIBRIUM PHASES The key words TRANSPORT and KINETICS were tested in the previous section 4 2 1 Cation exchange reactions Verification Problem 6 Transport of heavy metals subject to multiple cation exchange CATEXCH In this problem the transport of ten ions Al Br Ca Cd Cl K Mg Na Pb Zn through a soil column is modelled Initial and inflow concentrations of the ions are given in Table 4 4 The cation exchange capacity is equal to 0 011 mol cell The soil core has a length of 8 cm and is discretized into 40 cell
18. ea eit taie Cd eer fioqutit feda beau uai ERR 34 Figure 4 5 Distribution of Conta Contb and Contc versus depth after 250 500 and 1000 d for different simulations as defined in Table 4 3 Verification Problem 3 ve veo ado dpa ctae LN M LU ne EE 38 Figure 4 6 Distribution of pH dissolved Ca Cl Zn and Cd concentrations exchangeable Ca Na Zn and Cd concentrations versus depth after three days of infiltration for Verification Problem 6 sse 41 Figure 4 7 pH Ca Cl Zn and Cd concentrations in the effluent for Verification Problem 6 Full line HP1 dashed line PHREEGQYO sse 42 Figure 4 8 Distribution of pH Si Al amorphous SiO and Gibbsite versus depth after 150 days of infiltration for Verification Problem 7 Full line HP1 dashed Ing DEHRBEQU atis orat s t sd tui bu d e MA een vip fied bate 43 Figure 4 9 Selected results for cation and heavy metal transport with a pH dependent cation exchange complex a and b pH and Cd concentrations in outflow at 50 cm depth respectively c and d distributions of pH and Cd concentrations versus depth after 0 3 0 5 and 0 7 y respectively e and f distribution of the fraction of deprotonated sites 1 H sites and sites with Cd versus depth after 0 3 0 5 and 0 7 y respectively sss 48 Figure 4 10 Effect of pH and Cl concentration on Cd concentrations of leaching water CEE TE 20 em depth iae qe etre ea aout Cages Pe qu
19. input files should be placed in the same folder except the database file for which the path and name is defined in the species in file J The simulation is started by activating the executable Output files are created in the same folder as the input files The content of the output files created by HYDRUS 1D are described in section 11 of Sim nek et al 1998 The output file of PHREEQC contains user defined output data the SELECTED OUTPUT key word in the input file see Parkhurst and Appelo 1999 Part of the created output can be viewed through the interfaces 24 25 4 Verification problems In this section we describe different test examples of the coupled HP1 model We present examples having different levels of complexity The following problems were solved with HP1 for verifying the numerical correctness of the coupling procedure 1 transport problems with single or multiple components subject to sequential first order decay will be compared with simulations using HYDRUS separately ii multicomponent transport problems for steady state flow will be compared with simulations using PHREEQC separately and iii several more complicated problems will be compared with simulations using the CRUNCH model The CRUNCH model was a result of further development of the GIMRT OS3D codes of Steefel and Yabusaki 1996 see http www csteefle com Several of the verification examples as well as additional applications for transient flow were pr
20. of a given master species on concentrations or activities of the master species on concentrations or activities of secondary aqueous species or on external factors such as temperature and biomass It is possible to define rather complex sequential and parallel kinetic pathways e g Steefel 2000 This flexibility in PHREEQC is possible due to an embedded BASIC interpreter which permits one to define rate expressions in a general way in the input file see the section Numerical method and rate expressions for chemical kinetics in Parkhurst and Appelo 1999 Typical examples of homogeneous kinetic reactions are oxidation reduction reactions radioactive decay and degradation reactions Another type of homogeneous kinetic reaction is the production or consumption of a particular master species the consumption production of the specific secondary species is also possible Only mineral dissolution precipitation can be considered as a heterogeneous kinetic reaction The BASIC interpreter in PHREEQC allows one to use a broad range of reaction rate expressions 2 5 Multicomponent reactive transport Multicomponent reactive transport system may be viewed as consisting of Nm transport equations for the master species obc Oc Ogc j o wj af C Qm n x x SC R 2 23 14 j 1 Nm and Ns transport equations for the secondary aqueous species 00c wOC 096 at glen S ax Sc F R 2 24 i 1 Nsa The te
21. of data input 3 1 Input data The following separate input files are required to run the coupled HP1 model HYDRUS 1D selector in contains the following blocks A Basic Information Water Flow Information Time Information Root Growth Information Heat Transport Information EGER rn y Solute Transport Information G Root Water Uptake Information profile dat contains the following block H Nodal Information except the initial total aqueous concentrations atmosph in contains the following block I Atmospheric Information PHREEQC xxxxxxx xxx is the database file containing thermodynamic data for aqueous species pure phases and exchange reactions e g phreeqc dat minteq dat phreeqc in contains information about the geochemical reactions HP1 species in contains a list of master species to be transported The data and format of the different input files and blocks are described in section 10 of the HYDRUS ID manual im nek et al 1998 for input files related to HYDRUS 1D selector in profile dat and atmosph in and in Parkhurst and Appelo 1999 for input files related to PHREEQC We recommend that users refer to original manuals in order to complete the data input Specific guidelines for the coupled HP1 code are e The species in file contains on the first line the path to the geochemical database and then a list of master species to be transported An example of the species in file is 22 c program files p
22. present in the soil solutions resulting in the formation of aqueous complexes with the heavy metals 45 Table 4 5 Soil hydraulic properties and cation exchange capacities of five soil layers Horizon Layer 6 6 a n K l Cation exchange thickness cm cm day capacity cm eq 1000 cm soil A 13 0 065 0 476 0 016 1 94 93 0 5 0 0183 E 10 0 035 0 416 0 015 3 21 311 0 5 0 0114 Bhl 5 0 042 0 472 0 016 1 52 39 0 5 0 0664 Bh2 5 0 044 0 455 0 028 2 01 860 0 5 0 0542 Bh C 17 0 039 0 464 0 023 2 99 1198 0 5 0 0116 Table 4 6 pH and solution concentrations used in the simulation umol I Solution pH Na K Ca Mg Br Cl Cd Pb Zn 0 28 cm depth 8 5 4019 120 98 5 780 0 08 2 5 50 28 50 cm depth 8 5 4540 120 98 5 780 0 00 0 0 Applied water 3 5 127 5 120 98 5 780 0 00 0 0 Concentration of Na is adjusted to obtain the desired pH The soil profile is assumed to contain five distinct layers with different soil hydraulic properties and cation exchange capacities Table 4 5 gives thicknesses of the different horizons parameters of van Genuchten s equations for the water retention and hydraulic conductivity functions and the total cation exchange capacities The higher exchange capacities of the Bhl and Bh2 horizons reflect their enrichment with immobilized organic matter Flow is assumed to be steady at a constant flux of 0 05 m day 18 25 m year which causes the soil profile to be unsaturated water contents vary between
23. solution Element Initial composition Composition of the of the pore water inflowing alkaline solution mol kgw mol kgw K 5 42 10 1 61 107 Na 3 28 10 6 52 10 Ca 2 16 10 2 24 10 Mg 3 57 10 1 00 10 Cl 2 64 10 1 00 10 S 9 20 10 1 00 10 Al 1 00 10 1 00 10 Si 4 40 10 1 00 10 C 3 96 10 4 21 10 pH 7 90 13 20 1 From Table 2 Part III Adler 2001 From Table 6 Part III Adler 2001 MASTER SPECIES SOLUTION SPECIES EXCHANGE MASTER SPECIES EXCHANGE SPECIES PHASES and RATES The last keyword is used to define the specific rate reactions for the dissolution of primary minerals except calcite with Eq 4 7 translated into the BASIC programming language for each mineral for details see the PHREEQC 2 manual Parkhurst and Appelo 1999 The input file for water flow and solute transport can again be easily created with the graphical interface of HYDRUS 1D Nine components i e Ca K Na C S Cl Mg Si and Al are transported in this problem Only parameters related to transport are defined in HYDRUS 1D i e the dispersivity and the aqueous diffusion coefficient This information is stored in input files selector in and profile dat Note that these files also contain information for the steady state inflow flux and concentrations of the nine components in the infiltrating water The link between HYDRUS 1D and PHREEQC is defined in the input file species in which conta
24. the components primary species or master species in the terminology of PHREEQC and ii the secondary species The number of master species Nm equals the total number of species minus the number of reactions if the reactions are written stoichiometrically independent none of the reactions is a linear combination of the others The number of secondary species N is then defined as the number of species minus the number of master species Thus each reaction can be written in the canonical form 11 Yvan a 2 14 where i 1 Ns Ns is the number of secondary species 47 and A are the chemical formula of the master and secondary species respectively and v are the stoichiometric coefficients in the reaction In the remainder of this section we present the mass actions relations for different reactions aqueous speciation ion exchange mineral precipitation dissolution and the activities of the aqueous and exchange species 2 4 1 Homogeneous aqueous reactions For an aqueous complexation reaction Eq 2 14 is written as Nn l m l 2 55 Ar A 2 15 where the superscript indicates the liquid phase i 1 Nsa where N is the number of aqueous secondary species For equilibrium conditions the mass action equation is K ye jire 2 16 j l where K is the equilibrium constant for reaction 2 15 y is the activity coefficient of the ith aqueous complex and 77 is the activity coefficient of the jth master
25. 1 i tarit o f P f it i 4o10 M s amp dn te le on ff E eoad uj 9 8 60 6 pb on d A i MW h di zt n Hh MN if KU P m 3 Son n M r1 h MI bora d roy bot a m EL n ui i pt i voa dna d ia i D i g HZ ag ie tae 3 1i ie i t n HE le oppi sagene arag S 7 Chie dal pear HEN ooo M uti iu ur de BA EU T A ic D s a i IM hy d yon Wy S ys a Bw LM amp ooy y Ge F pum f u o vTWm 8 oni 4 0x10 3 0x10 1972 1974 1976 1978 1980 1982 1972 1974 1976 1978 1980 1982 Time year Time year 2x10 0x10 0x10 x amp 1x10 E 2x10 E o o E E 4x10 o Q Q 2x10 6x10 8x10 3x10 1972 1974 1976 1978 1980 1982 1972 1974 1976 1978 1980 1982 Time year Time year igure 5 2 Time series of inflow minus ouftlow dashed line and change in the total amount present in the soil profile dots for Cl Na Ca and Cd 100 5 3 ud ft h where Mi and Mi are the total amount of the jth master species in the ith cell at time and at the beginning of the simulation respectively Figure 5 3 shows relative mass balance errors for Na Ca and Cd The errors were consistently smaller than 1 indicating that mass balance errors were small in this HP1 simulation using the non iterative sequential approach r N cells max 3 1 i J sn Mj Mj t i for solving the coupled transport reaction equations 5 1 2 2 Simulation results As an illustration two dimensional time dept
26. 1237 1247 73 Prommer H 2002 A reactive multicomponent transport model for saturated porous media User s Manual Version 1 0 Contaminated Land Assessment and Remediation Research Centre University of Edinburgh UK 132 p Schaerlaekens J D Mallants J Sim nek M Th van Genuchten and J Feyen 1999 Numerical simulation of transport and sequential biodegradation of chlorinated aliphatic hydrocarbons using CHAIN 2D Hydrol Process 13 2847 2859 Sim nek J and D L Suarez 1993 Modeling of carbon dioxide transport and production in soil 1 Model development Water Resour Res 29 487 497 Sim nek J and D L Suarez 1994 Two dimensional transport model for variably saturated porous media with major ion chemistry Water Resour Res 30 1115 1133 Sim nek J and D L Suarez 1997 Sodic soil reclamation using multicomponent transport modelling ASCE J Irrig and Drain Engineering 123 5 367 376 Sim nek J and A J Valocchi 2002 Geochemical Transport In Methods of Soil Analysis Part 1 Physical Methods Chapter 6 9 Eds J H Dane and G C Topp Third edition SSSA Madison WI 1511 1536 Sim nek J M Sejna and M Th van Genuchten 1998 The HYDRUS 1D software package for simulating the one dimensional movement of water heat and multiple solutes in variably saturated media Version 2 0 IGWMC TPS 70 International Ground Water Modeling Center Colorado School of Mines Golden CO 202 p Sim n
27. 1987 A 0 b b O b 0 2 12 where bi b2 and b are empirical parameters ML T K H 2 3 2 Initial and boundary conditions for heat transport To solve Eq 2 9 initial and boundary conditions must be specified The initial temperature as a function of space at time zero must be defined Possible boundary conditions again include first type or Dirichlet type boundary conditions defining a prescribed boundary temperature and third type or Cauchy type boundary conditions defining a prescribed heat flux through the boundary At an impermeable boundary q 0 or at a boundary where water flows out of the domain a third type boundary condition reduces to a second type Neumann type boundary condition At the soil surface the temperature can be represented using a sine function Kirkham and Powers 1972 m 2at 7 T Ty Asin 221 1 2 13 where p is the period of time necessary to complete one temperature cycle T Tavg is the average temperature at the soil surface K and A is the amplitude of the sine wave K 2 4 Geochemical reactions In general species present in a system are related to each other by a set of reaction equations It is possible to write the various reaction equations in terms of a limited set of independent components The latter group permits one to define the stoichiometry of the system and are independent of each other Morel and Herring 1993 The species are thus divided in two groups i
28. 2 HP1 3 HP1 4 Input simulation parameters Atinit sec 864 864 864 864 864 Atmin sec 0 0864 0 0864 0 0864 0 0864 0 0864 tina Sec 86400 30000 10000 86400 86400 2 2 2 0 4 0 2 AER 1 1 1 1 1 Simulation results N 9070 9960 13585 9761 11775 At min sec 864 864 864 864 864 Atmax Sec 86400 30000 10000 86400 57600 Pe max 1 09 1 20 1 17 1 08 1 04 Cr pa 0 47 0 99 0 89 0 47 0 24 Note that Courant numbers are calculated differently R is included in the denominator Number of time steps for the complete simulation 37 An overview of different numerical criteria is given in Table 4 3 The HP1 simulation parameters temporal and spatial steps for this particular problem had no significant influence on the accuracy of the simulation results For purposes of computational efficiency it is however always better to use the stability criterion based on the performance index o This stability criterion ensures that small time steps are used when flow is rapid but allows for larger time steps when flow is slow This is done by modifying time steps such that the Courant number obeys the inequality given by Eq 4 6 Note that HP1 requires a lower performance index five times smaller than HYDRUS 1D to obtain similar results Conta 0 8 E E c c 06 2 2 D 0 4 o o c c o o O O 0 2 0 0 20 40 60 80 100 00 Depth cm Depth emn 0 08 E E c 0 06 c o 5 0 04 D o o c
29. 53 E mail djacques sckcen be Jirka Sim nek Tel 1 951 827 7854 Fax 1 951 827 3993 E mail jiri simunek ucr edu LIST OF TABLES Table 2 1 Description of variables that are used in both HYDRUS 1D and PHREEQC or newly defined in HYDRUS ID a iiie mide eie e i texere 20 Table 4 1 Transport and simulation parameters for Verification Problem 3 30 Table 4 2 Soil hydraulic transport and reaction parameters for Verification Problem 36 Table 4 3 Input simulation parameters and numerical results for Verification Problem d pastu ca Mte bbc aM Ld in o D e 37 Table 4 4 Initial and inflow concentrations for Verification Problem 6 sss 39 Table 4 5 Soil hydraulic properties and cation exchange capacities of five soil layers 45 Table 4 6 pH and solution concentrations used in the simulation umol T 45 Table 4 7 Overview of aqueous equilibrium reactions and corresponding equilibrium constants data from phreeqc dat database Parkhurst and Appelo 1999 45 Table 4 8 Log K parameters for multi site exchange complex sse 46 Table 4 9 Mineralogical reactions Log K molar volume volume percent and moles present in Opalinus Clay and secondary minerals considered 51 Table 4 10 Parameters for the kinetic dissolution reactions Eq 4 7
30. 8 in this example the upper boundary condition varies with time and consists of daily values of precipitation and evaporation no plants are assumed to be present and transpiration and root water uptake are not considered Soil hydraulic and physical parameters Table 5 1 of the dry Spodosol located at the Kattenbos site near Lommel Belgium were taken from Seuntjens 2000 Tables 3 1 and 7 1 The cation exchange complex was assumed to be associated solely with organic matter The cation exchange capacity hence is directly related to the amount of exchangeable protons on the organic matter taken to be 6 meq g of the organic matter proton dissociating groups on fulvic acids are 6 10 meq g and 4 6 meq g on humic acids Tipping 2002 The cation exchange complex was assumed to consist of six different exchange sites as defined in Table 4 8 The upper boundary condition precipitation and evaporation was defined using meteorological data from the Brogel station Belgium for a period of 10 years 1972 1981 Table 5 1 Soil hydraulic and other properties of six soil horizons Horizon Depth p Organic 0 0 a n K cm g cm Carbon m ms A 0 7 1 31 2 75 0 065 0 48 1 6 1 94 1 1 10 E 7 19 1 59 0 75 0 035 0 42 1 5 3 21 3 6 10 Bhl 19 24 1 3 4 92 0 042 0 47 1 6 1 52 4 5 10 Bh2 24 28 1 38 3 77 0 044 0 46 2 8 2 01 1 0 10 BC 28 50 1 41 0 89 0 039 0 46 2 3 2 99 1 0 10 CI 50 75 1 52 0 12 0 03 0 42 2
31. D of 8 cm an immobile water content of 0 05 cm cm and a first order exchange coefficient of 0 0125 d Precipitation and evaporation rates were typical for the Campine region in Belgium The soil profile was discretized into 100 elements of 1 cm each Chloride was applied during the first 53 days of the simulation with a concentration of 0 1 mmol 1 Results in Figure 4 2 show a perfect match between HYDRUS 1D and HP1 From these results we conclude that the transfer of water contents and concentrations between the transport and reaction modules was done correctly 1E 005 f Sa x oer oe a ee 7 NC 1E 006 Ao L E t E d M i g HYDRUS 1D 1E 007 e i e HP1 oO 0 100 200 300 Time d Figure 4 2 Outflow concentrations for Verification Problem 2 The full line was generated with HYDRUS 1D while dots were obtained with HP1 28 4 1 2 Transport of a nonlinearly adsorbing solute subject to first order decay This section considers numerical simulations of a nonlinearly sorbing solute undergoing first order decay during steady state water flow We assume only heterogeneous reactions of a contaminant Cont with a sorbing surface Sor Adsorption is assumed to be instantaneous and described with the Freundlich equation Sk CT 4 1 where S is the adsorbed concentration M M C is the aqueous concentration ML Ka is an empirical adsorption coefficient Ma0 L M and np is the empi
32. DISCRETIZATION INITIAL CONDITIONS and BOUNDARY CONDITIONS flow and transport properties FLOW and TRANSPORT and output OUTPUT The keyword POROSITY defines the porosity Since the porosity is different for each soil layer porosities are read from the porosityl dat input file Numerical issues for solving the equations are defined in RUNTIME The read saturation identifies the file containing the saturation degree for each cell Information for this file saturation1 dat was obtained from an initial HYDRUS 1D run to find the steady state water content profile 4 3 1 3 Comparison between CRUNCH GIMRT CRUNCH OS and HP1 Simulations were performed with three different models 1 the global implicit option in the CRUNCH model CRUNCH GIMRT 2 the sequential non iterative option in the CRUNCH model CRUNCH OS and 3 the coupled HP1 model which also runs in a sequential non iterative mode The maximum time step for the latter two simulations was chosen such that both were run at a similar Courant number 0 5 Figure 4 9 shows selected simulation results Infiltration of the low pH solution causes a decrease in the cation exchange capacity i e an increase in the number of protonated sites on the six cation exchange complexes see Figure 4 9e This relation between pH in the soil solution and the deprotonated sites on the multi site cation exchange complex is apparent when comparing Figure 4 9c and Figure 4 9e This leads to a desorpt
33. Print int Step No int IMobil float ThImob char Get_Names int species_number 2 7 2 HYDRUS 1D The original Hydrus 1D code contained the following files see section 12 of Sim nek et al 1998 Hydrus for Input for Material for Output for Sink for Solute for Temper for Time for and Watflow for We added one additional file Exports h that interfaces Fortran with c Most of the changes were done in the main HYDRUS file where the main program was converted into a HYDRUS_MAIN subroutine From this subroutine we call the following c functions that are located in the Hydr tr c file call Initialization NumNP IMobil TempN ThOld ThImob NSD NS Conc Sorb tInit call Get Concentrations NumNP NSD NS Conc Sorb IMobil ThImob call run HYDRUS reactions NumNP NSD NS Conc t TempN ThNew Sorb dt iPrint TLevel IMobil ThImob and one new Fortran subroutine call PhreeqcMB NSD NS NumNP Conc Conc1 Sorb Sorbl1 x ThNew ThImob PhrExch The run_HYDRUS_reactions function calls PHREEQC to carry out all chemical and biological reactions and transfers information from HYDRUS to PHREEQC the Get_Concentrations function transfers information back to HYDRUS from PHREEQC the Initialization brings initial information from the PHREEQC part and PhreeqcMB calculates mass balances for the main components 20 Table 2 1 Description of variables that are used in both HYDRUS 1D and PHREEQC or newly defined in HYDRUS 1D Fortran variab
34. SSA Inc Madison Wisconsin USA 693 1008 Jacques D J Sim nek D Mallants M Th van Genuchten 2002 Multicomponent transport model for variably saturated porous media Application to the transport of heavy metals in 72 soils In Computational Methods in Water Resources Eds Hassanizadeh S M RJ Schotting W G Gray and G F Pinder XIVth International Conference on Computational Methods in Water Resources June 23 28 2002 Delft The Netherlands Elsevier Developments in Water Science 47 555 562 Jacques D J Sim nek D Mallants and M Th van Genuchten 2003 The HYDRUS PHREEQC multicomponent transport model for variably saturated porous media Code verification and application MODFLOW and More 2003 Understanding through Modeling Conference Proceedings ed E Poeter Ch Zheng M Hill and J Doherty International Ground Water Modeling Center Colorado School of Mines 23 27 Jaynes D B A S Rogowski and H B Pionke 1984 Acid mine drainage from reclaimed coal strip mines 1 Model description Water Resources Research 20 233 242 K llvenius G and C Ekberg 2003 TACK a program coupling chemical kinetics with a two dimensional transport model in geochemical systems Comp and Geosci 29 511 521 Kirkham D and W L Powers 1972 Advanced Soil Physics John Wiley amp Sons New York NY Lichtner P C 1996 Continuum formulation of multicomponent multiphase reactive transport In Reactiv
35. Verification Problem 9 ALKALINE 35 ense er PEHER IP FERE EOM REEL sous onsotbaonpenssenpsveusionssniaanbe 50 5 Examplezidiiccitviiadi E iE Qoo ecd ae ee oio oL E coo edi FAL P A 59 5 1 Long term transient flow and transport of major cations and heavy metals in a soil profile HEAVMET 4isetiee see eene g epe tla Re C R epe ERES Pee dE Seda RAT ERR Eie da ees 59 5 1 1 Problem definition and governing chemical reactions 59 5 1 2 Analysis of the simulation results 4 eee ee ecce ee eere eee eerte eee n een 61 SE Concluding iari P ad CR 67 Referentes ooi olx e RA netaa ssas Ee eSEE SPEO ORE ERE EEDE EEEa OESS REA DR CREBRA RAN ROI DER 71 Acknowl deements 52 5222 ciae ie rii ormai ose Rr Eb Ed c opsel eeu Dose M edu eCk ossos sessio DRM op UNE OR UNS 69 Refereneesacohsaolianedtonr picta oU EOS GSP E COS o assins oossoo sse ossen SUPE 71 APPENDIX A LISTING OF NOTATION eere ee eere ette ene enne ta seta sete sete esses etas e tasa 75 APPENDIX B RELEASED VERSIONS AND BENCHMARKING eere 79 iii ABSTRACT Jacques D and J im nek 2004 User manual of the Multicomponent variably saturated transport model HP1 Version 1 0 Description Verification and Examples SCK CEN Mol Belgium BLG 998 79 p This report describes a new comprehensive simulation tool HP1 HYDRUSID PHREEQC that was obtained by coupling the HYDRUS 1D one dimensional vari
36. Wagenet R J and J L Hutson 1987 LEACHM Leaching Estimation And CHemistry Model A process based model of water and solute movement transformations plant uptake and chemical reactions in the unsaturated zone Continuum 2 Dept of Agronomy Cornell University Ithaca New York NY Walter A L E O Frind D W Blowes C J Ptacek and J W Molson 1994 Modeling of multicomponent reactive transport in groundwater 1 Model development and evaluation Water Res Res 30 3137 3148 Xu T 1996 Modeling nonisothermal multi component reactive solute transport through variably saturated porous media Civil Engineering School University of La Coru a Coru a Spain 310 p Yeh G T and H P Cheng 1999 3DHYDROGEOCHEM A 3 dimensional model of density dpeendent subsurface flow and thermal multispecies multicomponent HYDROGEOCHEMit cal transport EPA 600 R 98 159 150p Zheng C 1990 A modular three dimensional transport model for simulation of advection dispersion and chemical reaction of contaminants in groundwater systems Technical report U S Enviromental Protection Agency Zheng C and P P Wang 1998 MT3DMS A modular three dimensional transport model Technical report U S Army Corps of Engineers Waterways Experiment Station Vicksburg Miss APPENDIX A LISTING OF NOTATION Ci im Cim AANAAA ANAND ws Di Eo m s Je Ko m Ka activity of ith secondary exchange species amplitude of te
37. a particular cation is independent of the presence of other elements in the soil solution or on the soil solid phases Their coefficients are constant and independent of pH other cations complexing elements and ligands in the soil solution Unlike HYDRUS 1D the coupled HP1 code can include the effect of these factors on adsorption and consequently on the migration of multiple cations by using the ion exchange reactions of PHREEQC In this example we consider the transport of several major cations Na K Ca and Mg and three heavy metals Cd Zn Pb through a 50 cm deep multi layered soil profile during unsaturated steady state flow Each soil layer has different soil hydraulic properties and cation exchange capacities CEC Table 4 5 The top 28 cm of the soil is assumed to be contaminated with the three heavy metals initial pH 8 5 while an acid metal free solution pH 3 infiltrates into the soil Table 4 6 Assuming that the cation exchange complex is associated solely with organic matter CEC increases significantly with increasing pH due to the acid base properties of its functional groups The higher the pH the more functional groups of the organic matter are deprotonated and thus the higher the cation exchange capacity This behaviour is represented by a multi site cation exchange complex consisting of six sites each having a different selectivity coefficient for the exchange of protons see Appelo et al 1998 Finally chloride is
38. ably saturated water flow and solute transport model with the PHREEQC geochemical code The HP1 code incorporates modules simulating 1 transient water flow in variably saturated media 2 transport of multiple components and 3 mixed equilibrium kinetic geochemical reactions The program numerically solves the Richards equation for variably saturated water flow and advection dispersion type equations for heat and solute transport The flow equation incorporates a sink term to account for water uptake by plant roots The heat transport equation considers transport due to conduction and convection with flowing water The solute transport equations consider advective dispersive transport in the liquid phase The program can simulate a broad range of low temperature biogeochemical reactions in water soil and ground water systems including interactions with minerals gases exchangers and sorption surfaces based on thermodynamic equilibrium kinetics or mixed equilibrium kinetic reactions The program may be used to analyze water and solute movement in unsaturated partially saturated or fully saturated porous media The flow region may be composed of nonuniform soils or sediments Flow and transport can occur in the vertical horizontal or a generally inclined direction The water flow part of the model can deal with prescribed head and flux boundaries boundaries controlled by atmospheric conditions as well as free drainage boundary conditions The g
39. ally the solute transport equation The approach assumes that 1 the temperature effect on chemical reactions is important so that updated temperature information should be used for the geochemical equilibrium and kinetic calculations and that ii the effect of advection water flow in the solute and heat transport equations is larger than the effect of concentration and temperature on the water flow equation Therefore the water flow equation should be solved before the solute and heat transport equations The weak coupling method is also invoked in our modelling approach The same solutions sequence was used in the original HYDRUS 1D model im nek et al 1998 and is followed in HP1 as well Different approaches also exist to solve the multicomponent reactive transport equation 2 30 This equation contains terms describing the physical transport of the component the first three terms on the right hand side and a term describing the geochemical reactions the fourth term on the right hand side The physical transport part and the geochemical reactions can be solved either simultaneously a global implicit or one step approach or sequentially an operator splitting two step or sequential approach For a discussion of these two approaches the reader is referred to Steefel and MacQuarrie 1996 and Mayer 1999 In our model we use the sequential approach Following Walter et al 1994 the solution space has three degrees of freedom spat
40. amount of Cl mmol 1000 cm soil c pH and d amount of aqueous Cd mmol 1000 cm soil in the A horizon 0 7 cm depth sorption sites for the other cations The increased supply of monovalent cations due to upward water flow during summer further increased competition for the sorption sites leading to relatively more sorption of the monovalent cations and higher concentrations of divalent cations and heavy metals in the soil solution Another important process is the complexation of the heavy metals with Cl Because of the presence of relative stable aqueous complexes e g CdCl the amounts of Cd Pb and Zn will increase when Cl concentrations are high i e during summer periods Figure 5 4d shows the amount of Cd in the aqueous phase Differences in the amount of aqueous Cd between the summer and winter periods are more than two orders of magnitude Since water contents are lower during summer the differences in the concentrations will be even larger Similar but less distinct patterns can be observed also deeper in the soil profile Figure 5 5 This example shows that atmospheric boundary conditions can have a significant effect on the amounts of Cd and other heavy metals in the solution and hence also on their bioavailability since uptake processes by plants and soil micro organisms are often concentration dependent i e passive uptake of solutes together with uptake of water by roots as well as active uptake which can
41. ated hydraulic conductivity are nonlinear functions of the pressure head Three analytical models are available in HYDRUS 1D to describe these soil hydraulic properties Brooks and Corey 1994 van Genuchten 1980 and Vogel and Cislerova 1988 Only van Genuchten s functions will be used in the verification problems documented in this manual O h 6 2 2 K 4 K si S di 2 3 where 6 is the residual water content L L 0 is the saturated water content L L a L n and m 1 1 n are shape parameters is the pore connectivity parameter Ks is the saturated hydraulic conductivity LT and S 0 0 6 9 is effective saturation These parameters can be obtained by direct measurement of the 6 h or K h relationships e g chapters 3 3 and 3 4 in Dane and Topp 2002 by inverse optimization Hopmans et al 2002 or indirect estimation from basic soil textural properties using pedotransfer functions The latter approach in HYDRUS 1D uses neural network predictions from textural data as implemented in the ROSETTA program of Schaap et al 1998 The following additional features related to the soil hydraulic properties are also included in HYDRUS 1D i hysteresis in the water retention characteristic and the hydraulic conductivity function i1 description of small scale spatial variability in the soil hydraulic properties by means of scaling factors and iii temperature dependence of the soil hyd
42. ation of a coupled geochemical transport model Agreement No 58 5310 0 F105 between SCK CEN USDA ARS and Agreement No C0 90001412 01 between SCK CEN University of California Riverside The work of Dr D Jacques was supported by SCK CEN R amp D program on Environmental Remediation project E022031 and CO91002 The project coordinator at SCK CEN was Dr D Mallants 70 References Aagaard P and H C Helgeson 1982 Thermodynamic and kinetic constraints on reaction rates among minerals and aqueous solutions I Theoretical considerations Amer J Sci 282 237 285 Adler M 2001 Interaction of claystone and hyperalkaline solutions at 30 C A combined experimental and modeling study Ph Dissertation Universitat Bern 120 p Appelo C A J and D Postma 1993 Geochemistry groundwater and pollution A A Balkema Rotterdam 536 p Appelo C A J E Verweij and H Schafer 1998 A hydrogeochemical transport model for an oxidation experiment with pyrite calcite exchangers organic matter containing sand Applied geochemistry 13 257 268 Appelo C A J M J J Van der Weiden C Tournassat and L Charlet 2002 Surface complexation of ferrous iron and carbonate on ferrihydrite and the mobilisation of arsenic Environ Sci Technol 36 3096 3103 Brooks R H and A T Corey 1966 Properties of porous media affecting fluid flow J Irrig Drainage Div ASCE Proc 72 IR2 61 88 Carsel R F and R S Pari
43. ation or reclamation of sodic soils Sim nek and Suarez 1997 they can not be used to simulate the transport and reactions of other chemical species such as trace elements radionuclides and other chemicals PHREEQC 2 can simulate a large number of low temperature geochemical reactions in water soil and ground water systems including interactions with minerals gases solid solutions exchangers and sorption surfaces based on thermodynamic equilibrium kinetics or mixed equilibrium kinetic reactions e g Appelo et al 2002 PHREEQC 2 also allows one to simulate one dimensional reactive transport using a mixing cell solution approach see Appelo and Postma 1999 for details e g Appelo et al 1998 Postma and Appelo 2000 This model hence can be used to simulate reactive transport during steady state flow including a wide variety of geochemical reactions However the model can not deal with solute transport during transient water flow conditions In recent years the PHREEQC geochemical code has been coupled to various groundwater water flow and solute transport models For example the PHT 3D model Prommer 2002 couples PHREEQC 2 with MT3DMS Zheng and Wang 1998 the latter being an extension of the single species transport simulator by Zheng 1990 However this coupled program deals only with solute transport and reactions while the groundwater flow field needs to be computed using a separated simulation model Xu 1996
44. be simulated with HP1 using Monod or Michealis Menten kinetics Moreover the high concentrations of the heavy metals occur during periods with the highest 65 0 4 _ 0 36 c amp C 0 32 O o 9 0 28 oO 0 24 02 o o e E Q 1 o 1x10 e e E E O 1x10 4 4 pP I a 4 2 4 38 o eo E o eo e 9 410 E E o 1x10 1972 1974 1976 1978 1980 1982 Time year Figure 5 5 Time series of water content Cl concentration mmol 1000 cm soil pH and Cd concentration mmol 1000 cm soil in the third horizon Bh1 horizon at 22 cm depth micro biological activity during the year the summer months HP1 hence seems to be a very attractive tool for evaluating the effects of atmospheric boundary conditions on the long term mobility and leaching of the heavy metals in and from soils While not further addressed here such evaluations possibly could also account for the transport of colloids such as micro organisms humics or fulvics 66 67 6 Concluding Remarks The HP1 code presented in this report couples HYDRUS 1D a one dimensional variably saturated water flow and solute transport model with PHREEQC a general biogeochemical code The new code is able to handle 1 transient water flow in variably saturated media 2 transport of multiple components 3 heat flow and 4 mixed equilibrium aqueous speciation mineral dissolution precipitation cation
45. ccount for changes in the volume of minerals and corresponding changes in porosity hydraulic properties and solute transport parameters 2 Description of the model The HP1 model is the result of coupling the variably saturated water flow and solute transport model HYDRUS ID im nek et al 1998 with the geochemical model PHREEQC Parkhurst and Appelo 1999 Details about the governing equations initial and boundary conditions parameterization and adopted numerical methods are given in manuals of the original HYDRUS 1D and PHREEQC models In this chapter we give only a very concise overview of these topics mainly to provide a better understanding of the verification problems described in chapter 4 For more detail users are encouraged to examine the original manuals 2 1 Water flow in the vadose zone 2 1 4 The water flow equation Combination of the mass balance equation with the Darcy Buckingham law results in the Richards equation that describes water flow in variably saturated porous media The one dimensional form of the Richards equation can be written as a prol onl a K h ex cosa S h 2 1 where is the water pressure head L 0 is the volumetric water content L L time T x is the spatial coordinate L positive upward S is the sink term L L T a is the angle between the flow direction and the vertical axis and K is the unsaturated hydraulic conductivity LT Both the water content and the unsatur
46. cesocesoosesocessocesocescosssocssoecssocesocssoosesscessccesoose 21 31 Input datdie vci EU emn ap nad ak AAE EEE EEE EEOSE ia 21 3 2 Running the model and output cccessccsesscccescssssscssnsccssnsscssncccesssccssscasesseceesseeeesses 23 4 Verification problenis aoesetoee eese p Eee ER ERENRAS Ie E LREF TOT Y EE RENR PANE RINQRMSIM SEINS IKE AERE S d 25 4 4 Modelling the transport of a single component or decay chain 25 4 1 Physical equilibrium and nonequilibrium transport of chloride 25 4 1 3 Transport of a nonlinearly adsorbing solute subject to first order decay 28 4 1 3 Transport of first order decay chain of nonlinearly adsorbing contaminants during transient flow ccsssccsssscssssscssssscssssscssssccsssscsssssscssssscssssscees 35 4 0 Multicomponent transport during steady state flow eeeeeeeeeeee eere 39 4 2 1 Cation exchange reactions socec etse eeee et perte NR MEE N E PEERS I RN WHERE REN EU En EO STRIS rovs 39 4 2 2 Mineral dissolution cisco eor erret rene ES URR pe EepE Uo HI p Eno CO e ERR HEDLRNGE RAPERE AER CUR 42 4 3 More complicated verification problems of HPI model ecce 44 4 3 1 Heavy metal transport in a medium with a pH dependent cation exchange complex Verification Problem 8 MCATEXCH e c eeeee esee eese ee ense seen seens 44 4 3 2 Infiltration of a hyperalkaline solution in a clay sample
47. components in the incoming water and spatial and temporal discretization parameters of the numerical problem The link between HYDRUS 1D and PHREEQC is defined in the Species in input file that contains names of nine elements the names of which must be the same as in the Phreeqc dat database Finally the initial solutions and exchange complexes are defined in the Phreeqc in input file Since the transport problem involves variable water contents with depth a SOLUTION keyword needs to be defined for nearly each cell In addition the cation exchange complex for each soil layer must be defined using the keyword EXCHANGE with each layer containing six exchange sites Also included in the input file are the keywords TRANSPORT to indicate that HYDRUS 1D will be used for transport modelling and SELECTED OUTPUT to specify the desired output Details about the keywords used in Phreeqc in can be found in the PHREEQC 2 manual Parkhurst and Appelo 1999 47 CRUNCH The database contains chemical reactions and equilibrium constants for the aqueous speciation reactions between end of primary and end of secondary and the exchange reactions between begin of exchange and end of exchange The input file reactive in defines the initial conditions for the five layers and the inflowing water CONDITION primary secondary and exchange species PRIMARY SPECIES SECONDARY SPECIES and ION EXCHANGE the discretization and the initial and boundary conditions
48. creased Br concentration to 0 umol per kilogram of water Only a small increase in the Cd breakthrough was obtained Figure 4 10 compared with the previous case All three cases produced essentially the same results using both CRUNCH GIMRT and HP1 CRUNCH GIMRT HP1 Cr 0 5 0 8 Cd umol kgw pH 8 5 Cl m 0 4 pH 8 5 0 0 2 0 4 0 6 0 8 1 Time y Figure 4 10 Effect of pH and Cl concentration on Cd concentrations of leaching water at the 50 cm depth 50 4 3 2 Infiltration of a hyperalkaline solution in a clay sample Verification Problem 9 ALKALINE 4 3 2 1 Problem definition and governing chemical reactions Adler 2001 described a column experiment in which a high pH Na Ca OH solution infiltrated into a core containing Opalinus Clay This clay is investigated in Switzerland as a potential host formation for long term disposal of high level radioactive waste Thury and Bossart 1999 Adler 2001 also defined a conceptual geochemical model to describe the outflow concentrations of this experiment A slightly adapted version of this conceptual geochemical model is used here as a benchmark comparison between CRUNCH and HP1 Infiltration of a high pH plume into the clay core leads to different reactions First the cation exchange complex associated with the clay minerals will interact with the compositional change in the aqueous solution Due to the high pH primary minerals in the clay will become
49. developed a coupled model TRAN PHREEQE that links PHREEQE with a two dimensional finite element code for flow and transport in aquifers PHREEQE is a precursor of PHREEQC without cation exchange surface complexation or kinetic reactions Another example of coupling between PHREEQC and a two dimensional transport model for water saturated conditions and constant temperature was presented recently by K llvenius and Ekberg 2003 in the TACK model Transport And Chemical Kinetics However this model does not simulate water flow A recent detailed review of available numerical multicomponent transport models is given by Sim nek and Valocchi 2002 including an overview of the mathematical equations representing the major chemical reactions and governing transport processes a brief discussion how these equations can be implemented in reaction multispecies transport models and a description of several possible applications This report describes a new comprehensive simulation tool HP1 HYDRUSID PHREEQC that results from coupling the HYDRUS 1D one dimensional variably saturated water flow and solute transport model with the PHREEQC geochemical code The model incorporates modules simulating 1 transient water flow in variably saturated media 2 transport of multiple components and 3 mixed equilibrium kinetic geochemical reactions The accuracy of the coupled HP1 model was evaluated by comparing simulation results with HYDRUS 1D PHREEQC for s
50. e keyword TRANSPORT together with the option cells indicating the number of cells No additional information is required after the keyword e g dispersivity and molecular diffusion are all defined in the input files of HYDRUS 1D although punch cells and punch frequency can be used to control the output The different input files are most conveniently created using the interactive graphic based user interfaces of the original models These interfaces are public domain and can be downloaded from the World Wide Web for both the HYDRUS 1D and PHREEQC models The HYDRUS 1D model is located at http www pc progress cz Fr Hydrus htm The main part of the input files required for the water flow heat and solute transport parts of the problem can be constructed using the HYDRUS 1D interface section B of im nek et al 1998 A graphical interface for PHREEQC can be found at http www geo vu nl users posv phreeqc index html Or at http wwwbrr cr usgs gov project GWC coupled phreeqci This interface is helpful for constructing the phreeqc in input file for the geochemical reactions in the problem 3 2 Running the model and output An extra file path dat must be placed in the same directory as the executable of the HP1 model ie HPl exe and two dynamically linked libraries with HYDRUS ID i e Hydrus dll and PHREEQC i e Phreeqc dll subroutines and functions The path dat file specifies the path to the input and output file folder All
51. e many processes that are interrelated and contain parameters that are dependent upon the state of the system Without attempting to be complete the following interactions may occur in natural systems e Effect of concentration and temperature on flow properties by affecting the water density and the viscosity and the surface tension at the air water interface e Effect of temperature on diffusion coefficients e Effect of temperature on the equilibrium constants and rate coefficients e Effect of water flow on solute and heat transport e Effect of mineralogical changes on water flow and solute transport parameters Not all of these interactions are included in the present model The effect of concentration on the flow properties and the effect of mineralogical changes on the water flow solute and heat transport equations are neglected in the present version of HP1 16 Yeh and Cheng 1999 discriminate between strong coupling in which the governing water flow solute transport and heat transport equations as well as the geochemical reactions are solved simultaneously and weak coupling in which the governing equations are solved sequentially In the latter method state variables obtained after solving a given equation are used to calculate properties and state variables in the next equation Yeh and Cheng 1999 used the following sequence for this purpose first solve the water flow equation then the heat transport equation and fin
52. e transport in porous media Eds P C Lichtner C I Steefel and E H Oelkers Reviews in Mineralogy no 34 Mineralogical Society of America Washington DC p 1 81 Mallants D D Jacques and T Zeevaert 2003 Modelling 26Ra Rn and 7 Pb in a proposed surface repository of very low level long lived radioactive waste ICEM 03 Conference Oxford 21 25 Sept 2003 ASME Proc On CD ROM session 25 paper 6 Mayer U 1999 A numerical model for multicomponent reactive transport in variably saturated porous media Ph D thesis Department of Earth Sciences University of Waterloo 285p Millington R J and J M Quirk 1961 Permeability of porous solids Trans Faraday Soc 57 1200 1207 Morel F and J Hering 1993 Principles and applications of aquatic chemistry New York John Wiley amp Sons Inc Parkhurst D L and C A J Appelo 1999 User s guide to PHREEQC Version 2 A computer program for speciation batch reaction one dimensional transport and inverse geochemical calculations Water Resources Investigations Reprot 99 4259 Denver Co USA 312 pp Perrochet P and D Berod 1993 Stability and the standard Crank Nicolson Galerkin scheme applied to the diffusion convection equation some new insights Water Resour Res 29 3291 3297 Postma D and C A J Appelo 2000 Reduction of Mn oxides by ferrous iron in a flow system column experiments and reactive transport modeling Geochim Cosmochim Acta 64
53. ek J D L Suarez and M Sejna 1996 The UNSATCHEM software package for simulating one dimensional variably saturated water flow heat transport carbon dioxide production and transport and multicomponent solute transport with major ion equilibrium and kinetic chemistry Version 2 0 Research Report No 141 U S Salinity Laboratory USDA ARS Riverside California 186pp Sim nek J T Vogel and M Th van Genuchten 1992 The SWMS 2D code for simulating water flow and solute transport in two dimensional variably saturated media Version 1 1 Res Rep 126 U S Salinity Lab Agric Res Serv U S Dept of Agric Riverside Ca Singer P C and W Stumm 1970 Acid mine drainage the rate determining step Science 167 1121 1123 Steefel C L 2000 New directions in hydrogeochemical transport modeling Incorporating multiple kinetic and equilibrium reaction pathways n Computational Methods in Water Resources XIII Eds L R Bentley J F Sykes C A Brebbia W G Gray and G F Pinder A A Balkema Rotterdam 331 338 Steefel C I and K T B MacQuarrie 1996 Approaches to modeling of reactive transport in porous media In Reactive transport in porous media Eds P C Lichtner C I Steefel and E H Oelkers Reviews in Mineralogy no 34 Mineralogical Society of America Washington DC p 83 129 Steefel C L and S B Yabusaki 1996 OS3D GIMRT Software for multicomponent multidimensional reactive transport user ma
54. ent for Verification Problem 6 Full line HP1 dashed line PHREEQC 4 2 0 Mineral dissolution Verification Problem 7 Transport with mineral dissolution MINDIS A 100 cm long soil column consisting of amorphous SiO2 and gibbsite AI OH is leached with a solution containing 5 107 mol I Si 1 10 mol I Al and 1 10 mol I Na to obtain an inflow pH of 11 15 Initial concentrations were 1 76 10 mol T Si 8 87 10 mol I AI and 1 10 mol Na corresponding to a pH of 6 33 In each l cm cell 0 015 mol amorphous SiO and 0 005 mol gibbsite is present The flow velocity is 2 cm day and the dispersivity is 1 cm Simulations are again carried out with both PHREEQC and HP1 Distribution of pH Si Al amorphous SiO and gibbsite versus depth after 150 days of simulation are presented in Figure 4 8 No significant differences between PHREEQC and PHI were apparent for pH Si Al amorphous SiO and gibbsite The keyword EQUILIBRIUM PHASES hence was coupled correctly in HP1 43 12 11 HP1 e PHREEQC I 10 Qa 9 8 0 20 40 60 80 100 Depth cm 0 003 0 0002 E E 0 00016 c c 2 o 0 002 i 0 00012 c c o o o o c c 9 Q 8E 005 o 0001 E Z Z 4E 005 e L a e 0 m 0 0 20 40 60 80 100 Depth cm 0 016 0 006 0 012 m ed J 0 004 E E S 0 008 S o lt E 0 002 0 004 amorf SiOo Gibbsite 0 0 0 20 40 60 80 100 0 20 40 60 80 100 Depth cm Depth cm Figu
55. er HP 1 produced relatively high amount of numerical dispersion The first order rate constant was chosen equal to 0 2 d or 2 31 10 sec We again compare HYDRUS ID results with PHREEQC simulations using different spatial discretizations and HP1 simulations using different maximum time steps Figure 4 4 The latter code is the reference Results of simulations with mixed equilibrium Freundlich adsorption and kinetic first order decay reactions were very similar for all three models thus supporting the same conclusions as in the previous case The transport of a nonlinearly adsorbing first order decaying contaminant was accurately modelled with HP1 when Courant numbers were reasonably small Cr smaller than 0 2 i e 10 times smaller than Cr from HYDRUS 1D Furthermore results obtained with HP1 were somewhat better compared to results obtained with PHREEQC for the same Peclet and Courant numbers 34 Ax 2 0 5 cm Concentration mol l 0 20 40 60 80 100 Depth cm HP1 2 Atna 0 25 d HP1 3 Atna 0 08 d HP1 4 Atna 0 04 d Concentration mol 0 20 40 60 80 100 Depth cm Figure 4 4 Depth profiles of Cont after 250 500 and 750 d for different simulations defined in Table 4 1 Verification Problem 4 with a first order decay coefficient of 0 2 d 35 4 1 3 Transport of first order decay chain of nonlinearly adsorbing contaminants during transient flow Problem definition In this example we cons
56. eviously reported by Jacques et al 2002 2003 4 1 Modelling the transport of a single component or decay chain This first group of verification problems concerns the transport of one solute either inert or adsorbing according to the Freundlich isotherm or of a sequential first order decay chain Although these problems can be solved directly with the HYDRUS 1D program we use them here to verify the coupling of HYDRUS 1D and PHREEQC with HYDRUS ID solving the transport part and PHREEQC the reaction part of the problem First we will test if the transfer of information i e solute concentrations and water contents between two modules of the coupled model i e HYDRUS 1D and PHREEQC is done correctly section 4 1 1 Second we will compare the numerical accuracy of both the coupled HP1 model and PHREEQC against the stand alone HYDRUS 1D using a single component transport problem involving either equilibrium using the Freundlich adsorption isotherm or kinetic first order decay reactions These geochemical reactions are coupled with transport within HYDRUS 1D using a one step or global implicit method i e various reactions are directly included into the governing transport equations In contrast a non iterative sequential approach is used in both the PHREEQC and HP1 codes In section 4 1 2 we will compare results of these three models i e HYDRUS 1D PHREEQC and HP1 for different Peclet and Courant numbers Finally the transport
57. exchange and kinetic geochemical reactions The non iterative sequential approach used to solve the coupled transport reaction equations generally worked very well The accuracy of the coupling approach was tested using nine verification examples with different degrees of complexity Overall compared to other models codes HP1 was found to be accurate when time steps were not too large The easiest and most efficient way to control the time step is by defining a relatively small maximum time step and or a small performance index both in HYDRUS 1D Our analyses indicate that the performance index should be about five time smaller than the value previously suggested for HYDRUS 1D simulations i e 0 4 rather than 2 The nine verification problems and the application example serve as guidelines to create the input files using the graphical user interfaces of HYDRUS 1D and PHREEQC The user interfaces also display the output created with the HP1 simulations 68 69 ACKNOWLEDGEMENTS The work of Dr J Sim nek was supported in part by SAHRA Sustainability of semi Arid Hydrology and Riparian Areas under the STC Program of the National Science Foundation Agreement No EAR 9876800 and the Terrestrial Sciences Program of the Army Research Office Terrestrial Processes and Landscape Dynamics And Terrestrial System Modeling and Model Integration Additional support was obtained through the bilateral agreement project Development and evalu
58. ger moles Nm number of moles of mineral m moles Ty 0 initial number of moles of mineral m moles Ny number of kinetically controlled homogeneous reactions Nn number of master species N number of minerals N number of secondary species Ns number of secondary aqueous species Ns number of secondary exchange species Ny number of master exchange sites P precipitation L P Peclet number P period of time necessary to complete one temperature cycle one day T q volumetric flux density L T qlo solute flux at the bottom of the soil profile of the jth master species M L T gin solute flux at soil surface of the jth master species ML T Q cumulative inflow of the jth master species M L cumulative outflow of the jth master species M L OF ion activity product of mineral m R decay rate ML T Rj source sink term due to geochemical reactions for ith species M L T Riim source sink term due to geochemical reactions for ith species in immobile M L T region Rin source sink term due to geochemical reactions for ith species in mobile region M L T Rn dissolution rate of a mineral M T AUT homogeneous equilibrium reactions of source sink term ML T gue heterogeneous equilibrium reactions of source sink term M L T R homogeneous kinetic reactions of source sink term M L T gus heterogeneous kinetic reactions of source sink term ML T S water sink term L L T Se effecti
59. h plots are shown in Figure 5 4 for the top 7 cm the A horizon for the water content pH mass of Cl mmol 1000 cm soil and mass of aqueous Cd mmol 1000 cn soil The alternation between precipitation wet conditions and evaporation dry conditions as dictated by the atmospheric boundary conditions clearly affected the dynamics of the water content Figure 5 4a including upward flow of water during the dry periods not shown The flow dynamics in turn significantly influenced the g Na Cd r 0 3 0 2 0 1 Het eat tt TOT TETTE hid TE Wr ILU METERS ug mp infi MT Viu mp ya pa n FLU i A l I i 1 bo di y 1972 1974 1976 1978 1980 1982 Time year 0 5 0 3 0 2 I I boa r 1 s D E NIU Fd agn T h jm A LO DEREN a ta A mul NU i uhi wn Hd PDA vut vy l P i lh Ua va jv i 0 1 1972 1974 1976 1978 1980 1982 Time year 1974 1976 1978 Time year Figure 5 3 Relative mass balance errors amp as a function of time for Na Ca and Cd 1980 63 1982 geochemistry near the soil surface The most mobile elements such as Cl and the monovalent cations H Na moved upwards during the evaporation periods leading to higher concentrations near the soil surface e g Cl in Figure 5 4b The increased proton concentration induced by the lower water content and the supply from deeper horizons in the soil profile produced lower pH values near the soil
60. here to have an equivalent problem as in HYDRUS 1D and MOL Cont the molality of the solute The third statement line 6 integrates the rate over the subinterval with the special variable TIME Finally the moles of reaction during the time interval are saved with the last special statement SAVE Note the negative sign on line 6 that results in a negative amount of moles saved in the last statement In general a positive sign represents decreasing amounts of a phase i e dissolution whereas a negative sign results in precipitation of that phase Consequently elements will enter the solution in the former case dissolution and will be removed from the solution in the latter case i e precipitation degradation or decay In this example the imaginary phase degradation is precipitating This is done to prevent the cessation of the kinetic reaction i e when the phase degradation is completely removed from the system The second keyword in Box 4 2 KINETICS defines the names of the rate expressions related to a specific cell In this example we have one rate expression called degradation Since degradation is used here as an imaginary phase and not a phase defined in the database the option formula is used to define the elements produced i e when the Box 4 2 PHREEQC input for first order decay reactions 1 RATES 2 degradation 3 start 4 10 rem parm 1 first order degradation coefficient sec 1 5 20 rate
61. hese in the vector of strings cRows nRows Then we call the check hydrus species function that checks whether these components are defined in the geochemical database Then instead of the original Transport function we call HYDRUS MAIN with three parameters which refer to c functions called from the HYDRUS Fortran code At the end we call function free string array which deallocates memory from cRows void stdcall HYDRUS MAIN void run HYDRUS reactions void Get Concentrations void Initialization int read text file const char file path name int nRows char cRows void free string array int nRows char cRows void check hydrus species int nRows char cRows Hydr tr c This file contains four new functions Initialization Get Concentrations Run Hydrus reactions and Get Names that are called from HYDRUS They are described below void Initialization int node number int IMobil float Temperature float Theta float ThImob int max species number int species number float Concetrations float Im Concetrations float Timelnit void Get Concentrations int node number int max species number int species number float Concetrations float Im Concetrations int IMobil float ThImob 19 void run HYDRUS reactions int node number int max species number int species number float Concetrations float time float Temperature float Theta float Im Concetrations float TimeStep int i
62. hreeqc phreeqc dat Na Cl Ca Mg where phreeqc dat is the database file This file forms a link between HYDRUS 1D and PHREEQC The master species in the file must be written in the same way as they are defined in the SOLUTION MASTER SPECIES block of the phreegc dat input file see Parkhurst and Appelo 1999 The order of the master species in the species in file refers to solute numbers in the se ector in files for the HYDRUS 1D model Note that a check is made whether or not names from the species in file correspond to master species from the geochemical database When the graphical user interface of HYDRUS 1D discussed later is used the number of master species is limited to 10 After the HYDRUS ID input files are created with the graphical interface the number of master species can be increased manually by making several changes in the Selector in file described in Table 10 6 of Sim nek et al 1998 increase the variable NS to the number of master species add records 9 and 12 for each additional solute and expand record 16 The Profile dat and Atmosph in files Table 10 8 and 10 9 respectively must also be expanded for the specified number of species We suggest to set solute transport parameters describing exchange with the solid or gas phases or describing degradation equal to zero It is important to use the same molecular diffusion coefficient for all master species see section 2 5 Note that this is not checked in the curren
63. ial temporal and chemical Physical transport is connected in the spatial and temporal domains and the geochemical reactions are only connected in the chemical domain The physical part coupled in space uncoupled over the components is obtained by solving Eq 2 30 without the reaction term 30C ac aqc Fa oD j sot NG 2 31 Ot ox ax ax and the chemical part uncoupled in space ie no transport but coupled over the components by simultaneously solving the equilibrium and kinetic geochemical reactions An overview of the coupled multicomponent reactive transport calculations is shown in Figure 2 1 The symbols used in this figure are n the nth time step H variables related to water flow pressure heads fluxes T variables concerning heat transport temperature C variables dealing with components and species in the system 17 G variables concerning the solid phase mineralogical composition exchange site surface site reactive surfaces p vector of parameters needed to solve the water flow equation Eq 2 1 p vector of parameters needed to solve the heat transport equation Eq 2 9 p vector of parameters needed to solve the solute transport equation Eq 2 31 p vector of parameters needed to solve the geochemical reactions Initial conditions with state variables H T C G Start a new time step n with state variables H T C G
64. ider the transport of three non linearly adsorbing contaminants Conta Contb and Contc that are involved in a sequential first order decay chain defined as Conta LL kc conta Contb Lk contb Con tc LL kc conte ut i 4 5 Kg np Jon Ka np contb Kate contc SoraConta SoraContb SoraContc where uw p are the first order rate constants connecting two contaminants Uw is the first order rate constant for Contc K4 and ny are the Freundlich isotherm parameters for the three contaminants and SoraConta SorbContb and SorcContc are the three surface species related to Conta Contb and Contc on three surfaces Sora Sorb and Sorc respectively Reaction parameters for the three contaminants are given in Table 4 2 Model simulations were carried out for a 1 m deep homogeneous soil profile during 1000 d assuming transient flow Upper boundary conditions were taken as daily precipitation rates representative of the Campine Region in Belgium Evaporation was neglected during the simulations The lower boundary condition was defined as free drainage A uniform initial pressure head of 60 cm was assumed for the entire soil profile For solute transport the following initial and boundary conditions were considered 1 low initial concentrations 10 mole I for all three contaminants 2 third type solute fluxes as the top boundary conditions with 1 0 1 and 0 mol 1 for Conta Contb and Contc respectively and 3 zero gradient bottom boundary cond
65. ilization criterion Transport parameters are given in Table 4 2 and simulation parameters 1 e temporal and spatial steps and performance index in Table 4 3 Distributions of Conta Contb and Contc versus depth are shown in Figure 4 5 at three times 250 500 and 1000 d Simulations obtained with HP1 1 were identical to the HYDRUS 1D results for Conta and Contb except at the leading edge of the concentration profile for the first print time for Conta Results for HP1 3 and HP1 4 shown a better agreement between the simulations than those of HP1 1 and HP1 2 Concentration profiles for Contc were also in very good agreement between different runs except for the peak concentrations which were smaller than those obtained with HYDRUS 1D Note again that simulations for HP1 3 and HP1 4 runs shows a better agreement than those for HP1 1 and HP1 2 runs Table 4 2 Soil hydraulic transport and reaction parameters for Verification Problem 5 Parameter Value Conta Contb ContC Hydraulic parameters 0 cm cm 0 078 0 cm cm 0 43 o cm 0 036 n 1 56 K cm d 24 96 Transport parameters D cm 1 p g cm 1 5 1g Reaction parameters K 05 25 5 ny 1 0 9 0 8 Log k 100 12 99 43 99 12 Lus i d 0 005 0 06 gt doa 0 02 Input for HP1 assuming the amount of adsorption sites S 10 Table 4 3 Input simulation parameters and numerical results for Verification Problem 5 Parameter HYDRUSID HPI HPI 1 HP1
66. ins names of the 9 components The phreeqc in file contains chemical information for the reactive transport problem Initial mineralogical and chemical conditions for each cell are defined with the following keywords SOLUTION initial chemical composition of the aqueous phase EXCHANGE size and initial composition of the cation exchange complex EQUILIBRIUM PHASES identifies the minerals in equilibrium with the aqeous phase and their initial amount KINETICS identifies the kinetic reactions with their parameters and their initial amount TRANSPORT to indicate that HYDRUS ID is used for transport modelling and SELECTED OUTPUT to specify the desired output Details on the keywords used in phreeqc in can be found in the PHREEQC 2 manual Parkhurst and Appelo 1999 54 CRUNCH The database contains chemical reactions and equilibrium constants for the aqueous speciation reactions between end of primary and end of secondary the minerals between end of secondary and end of minerals the kinetic dissolution or precipitation reactions of the minerals between begin of mineral kinetics and end of mineral kinetics and the exchange reactions between begin of exchange and end of exchange The input file reactive in defines the initial conditions for the clay and the inflowing water CONDITION the primary secondary and or exchange species and the minerals PRIMARY SPECIES SECONDARY SPECIES ION EXCHANGE and MINERALS the discretization
67. ion of Cd and other heavy metals and their subsequent leaching from the soil profile Cd leaching peaks after about 0 3 y Figure 4 9b with most Cd leached from the profile after 1 y Results obtained with CRUNCH OS and HP1 using Cr 0 5 showed very good agreement especially for the outflow curves This indicates that the coupling between HYDRUS 1D and PHREEQC was done correctly since the coupled model was tested in this relatively complicated example with a completely independent reactive transport model in contrast to the verification problems discussed earlier in sections 4 1 and 4 2 However the pH outflow Figure 4 9a and the concentration distributions versus depth Figure 4 9c f showed a small 48 increase in numerical dispersion as compared to simulations using the global implicit approach CRUNCH GIMRT 0 0 2 0 4 0 6 0 8 1 0 0 2 0 4 0 6 0 8 1 Time year Time year 2 5x10 2 0x10 1 5x10 1 0x10 1 H sites total CEC Cd sites total CEC 5 0x10 0 0x10 0 0 1 0 2 0 3 0 4 0 5 Depth m CRUNCH GIMRT HP1 Cr 0 5 4 CRUNCH OS HP1 Cr 0 1 Figure 4 9 Selected results for cation and heavy metal transport with a pH dependent cation exchange complex a and b pH and Cd concentrations in outflow at 50 cm depth respectively c and d distributions of pH and Cd concentrations versus depth after 0 3 0 5 and 0 7 y respectively e and f distribu
68. ition Modelling the decay chain with HP1 The decay chain was modelled by defining the rates of three kinetic reactions degrad_conta degrad contb and degrad contc to model the degradation of solutes Conta Contb and Contc respectively The definition of degradation rates is similar as in Verification Problem 4 To model the degradation of Conta into Contb the option formula for keyword KINETICS line 11 is changed to formula Conta 1 Contb 1 which indicates that when 1 mole of the imaginary phase degrad conta precipitates a negative rate 1 mole of Conta disappears and 1 mole of Contb appears A similar kinetic rate reaction is written for the transformations of Contb into Contc 36 Verification Problem 5 First order decay chain of nonlinearly adsorbing contaminants during unsteady flow SEASONCHAIN We now compare HYDRUS 1D results with four HP1 runs using different time stepping schemes The first two simulations were carried out with maximum time steps of 30 000 s HP1 1 and 10 000 s HP1 2 in order to obtain Courant numbers smaller than 1 based on flow velocities at the top boundary For the other two simulations we used the stabilisation criterion defined by Perrochet and Berod 1993 and implemented in HYDRUS 1D PeeCr lt a 4 6 where o is the performance index Simulations were performed with 0 4 HP1 3 and 0 2 HP1 4 The HYDRUS 1D module of HP1 automatically adjust the time steppeing to fulfil this stab
69. les c variables Type Description Initialization NumNP node number int Number of nodes IMobil IMobil int Mobile immobile water model TempN Temperature float Temperatures ThOld Theta float Initial water contents ThImob ThImob float Immobile water content NSD max_species_number int Max number of components NS species number int Number of components Conc Concetrations float Concentrations Sorb Im Concetrations float Immobile concentration tInit Timelnit float Initial time Get Concentrations see definitions above Conc Concetrations float Sorb Im Concetrations float Concl Concetrations float Sorb1 Im_Concetrations float run HYDRUS reactions t time float ThNew Theta float dt TimeStep float iPrint iPrint int TLevel Step No int PhreeqcMB Conc real Concl real Sorb real Sorb1 real x real PhrExch real Concentrations before chemical reactions Immobile concentrations before chemical reactions Concentrations after chemical reactions Immobile concentrations after chemical reactions Time Water content Time step Print time flag Step Number Component concentrations in the mobile phase before chemical reactions Component concentrations in the mobile phase after chemical reactions Component concentrations in the immobile phase before chemical reactions Component concentrations in the immobile phase after chemical reactions Nodal coordinates Change in the mass of a component 21 3 Description
70. lude first type or Dirichlet type boundary conditions defining a prescribed boundary concentration and third type or Cauchy type boundary conditions defining a prescribed boundary solute flux At an impermeable boundary i e where q 0 or at a boundary where water flows out of the domain the third type boundary condition reduces to a second type Neumann type boundary condition 2 3 Heat transport in the vadose zone 2 3 1 The heat transport equation The one dimensional heat transport equation neglecting water vapour diffusion is given by aC 0 T at Sb oT c qt _ 2 20 ar c 4 c sr 2 9 where A 0 is the apparent thermal conductivity of the soil MLT K and C 0 and C are volumetric heat capacities of the porous medium and the liquid phase respectively ML T K The volumetric heat capacity of the porous medium is estimated based on its constituents de Vries 1963 as follows C 0 C 0 C 0 C 0 C 8 2 10 where Cna Co and C are the volumetric heat capacities of the solid phase the organic matter and the gas phase respectively ML T K and 6 0 and 0 are the volumetric fractions of the solid phase the organic matter and the gas phase respectively L L The apparent thermal conductivity is defined as de Marsily 1986 A 0 4 0 amp C Jd Q 11 10 where f is the thermal dispersivity L and 2o 0 is the thermal conductivity of the soil defined as Chung and Horton
71. mineralogical geological and biological processes Cycles of precipitation and evaporation largely determine if contaminants remain near the soil surface Changes in the chemical composition or pH of the soil solution may impact the retention of heavy metals on organic matter or iron oxides Dissolution and precipitation of minerals generally buffer the transport of a solution with a different pH through the soil profile Simulation of these and related processes requires a coupled reactive transport code that integrates the physical processes of water flow and advective dispersive solute transport with a range of biogeochemical processes In this report we present a new code that resulted from the coupling of two existing codes HYDRUS ID im nek et al 1998 and PHREEQC 2 Parkhurst and Appelo 1999 The new code should significantly expand applicability of the individual codes by preserving most of their original features and capabilities HYDRUS 1D simulates water heat and solute movement in one dimensional variably saturated porous media for various boundary conditions including precipitation and evaporation A sink term accounting for water uptake by roots is also included in the model Solutes can be exchanged between the water and gas phase and may interact linearly or nonlineary with the solid phase assuming either equilibrium or nonequilibrium reactions between the dissolved and adsorbed solutes The only possible interaction between
72. mperature sine wave at soil surface chemical formula of ith secondary species chemical formula of ith secondary species initial reactive surface area of a mineral chemical formula of ith secondary exchange species chemical formula of ith aqueous secondary species chemical formula of ith mineral chemical formula of jth master species normalized water uptake distribution parameters to calculate soil thermal conductivity 7 1 2 3 number of equivalents of the i th secondary exchange species on the j th master exchanger aqueous concentration of ith species aqueous concentration of ith species in the immobile region aqueous concentration of ith species in the mobile region aqueous concentration of ith aqueous secondary species aqueous concentration of jth master species concentration of the sink term concentration of ith secondary species in the sink term concentration of jth master species in the sink term aqueous concentration volumetric heat capacity of the gas phase total concentration of jth master species Courant number total concentration of jth master species in the sink term volumetric heat capacity of solid phase volumetric heat capacity of organic matter volumetric heat capacity of porous medium volumetric heat capacity of liquid phase dispersion coefficient of ith species in liquid phase dispersion coefficient of ith species in free water longitudinal dispersivity evaporation soil water pressure head osmotic head
73. n Ax is decreased in the HP1 simulations 1 H meee A HP1 2 di P 1 Ax 2 1 cm P 2 Ax 0 5 cm E 06 P 3 6 Ax 0 25 cm 8 0 4 c o O 0 2 0 0 20 40 60 80 100 Depth cm 1 C H HP1 1 Atmax 1 d ae HP1 2 Hu Atma 0 25 d HP1 3 E o6 At 0 08 d HP1 4 Ata 0 04 d c 8 0 4 c o O 0 2 0 20 40 60 80 100 Depth cm Figure 4 3 Depth profiles of Cont after 250 500 and 750 d for different simulations defined in Table 4 1 for Verification Problem 3 Tests consider effects of grid size for PHREEQC a and effects of maximum time steps for HP b 32 Modelling first order decay with PHREEQC First order decay Eq 4 2 in PHREEQC is modelled using the keywords RATES and KINETICS The implementation in PHREEQC is shown in Box 4 2 The kinetic reaction itself is defined under the keyword RATES lines 1 8 In this example the reaction is called degradation line 2 Between start and end a Basic program is written for the kinetic reaction of the phase degradation consisting of standard Basic statements e g here only rem and special Basic statements for PHREEQC e g MOL SAVE TOT and TIME The first statement of the Basic program line 4 is only a comment indicating the meaning of the first parameter The second statement evaluates the rate equation Eq 4 2 with parm 1 being the value of 44 TOT water the amount of water in the cell which is included
74. nt physical nonequilibrium transport of chloride Verification Problem 2 HYDRUS 1D TRANSCL Steady state transport of a nonlinearly adsorbing contaminant Verification Problem 3 HYDRUS 1D PHREEQC STADS Steady state transport of a nonlinearly adsorbing contaminant subject to first order decay Verification Problem 4 HYDRUS 1D PHREEQC STDECAY Transport of first order decay chain of nonlinearly adsorbing contaminants during unsteady flow Verification Problem 5 HYDRUS ID SEASONCHAIN Transport of heavy metals subject to multiple cation exchange Verification Problem 6 PHREEQC CATEXCH Transport with mineral dissolution Verification Problem 7 PHREEQC MINDIS Heavy metal transport with a pH dependent cation exchange complex Verification Problem 8 CRUNCH MCATEXCH Infiltration of a hyperalkaline solution in a clay sample Verification Problem 9 CRUNCH ALKALINE Overview of examples Long term transient flow and transport of major cations and heavy metals in a soil profile HEAVMET
75. nual and programmer s guide PNL 11166 Pacific North west Laboratory Richland Washington 74 Stolk A P 2001 Landellijk Meetnet Regenwatersamenstelling Meetresultaten 1999 RIVM The Netherlands RIVM Rapport 723101 056 61 p Suarez D L and J Sim nek 1997 UNSATCHEM Unsaturated water and solute transport model with equilibrium and kinetic chemistry Soil Sci Soc Am J 61 1633 1646 Tipping E 2002 Cation binding by humic substances Cambridge University Press Cambridge UK 434 p Truesdell A H and B F Jones 1974 WATEQ A computer program for calculating chemical equilibria of natural waters J Res U S Geol Surv 2 233 274 Thury M and Bossart P 1999 The Mont Terri rock laboratory a new international research project in a Mesozoic shale formation in Switzerland Eng Geol 52 347 359 van Genuchten M Th 1980 A closed form equation for predicting the hydraulic conductivity of unsaturated soils Soil Sci Soc Am J 44 892 898 van Genuchten M Th 1985 Convective dispersive transport of solutes involved in sequential first order decay reactions Comp Geosci 11 159 147 van Genuchten M Th and P J Wierenga 1976 Mass transfer studies in sorbing porous meida I Analytical solutions Soil Sci Soc Am J 40 473 480 Vogel T and M Cislerova 1988 On the reliability of unsaturated hydraulic conductivity calculated from the moisture retention curve Transport in Porous Media 3 1 15
76. o et al 1998 An overview of the considered aqueous equilibrium reactions is given in Table 4 7 The role of chloride as a complexing agent is described by reactions 10 through 13 Other geochemical reactions that are considered are the heterogeneous multi site 10n exchange reactions The exchange coefficients for major cations and heavy metals were assumed to be the same for all exchange sites Table 4 8 gives parameters for this multi site exchange complex 4 3 1 2 Database and model input HP1 Chemical reactions and their equilibrium constants are stored in the database phreeqc dat The following keywords are used for this problem SOLUTION MASTER SPECIES SOLUTION SPECIES EXCHANGE MASTER SPECIES and EXCHANGE SPECIES The input file for water flow and solute transport is relatively straightforward and can be easily created the standard way with the graphical interface of HYDRUS 1D The problem involves 9 components that are transported i e Na K Mg Ca Cl Br Cd Zn and Pb Only parameters related to solute transport i e the dispersivity and the aqueous diffusion coefficient are defined in HYDRUS 1D All other solute transport parameters except for the Freundlich exponent which is equal to one are set to zero This information is stored in the HYDRUS 1D input files Selector in and Profile dat Note that these files also contain information about the steady state flux 0 05 cm day the boundary concentrations of the nine
77. obtained with HP1 were found to be identical to those obtained with HY DRUS 1D as illustrated by the outflow concentrations in Figure 4 1 1 a ee oeo oo e e e o So e o ee t eo E ett 2 x o eo oe 2 M ge T v 2 n rd equilibrium ee HYDRUS 1D s EC nonequilibrium 0 6 ry HP1 Concentration mol l 0 10 20 30 40 50 Time d Figure 4 1 Outflow concentrations for Verification Problem 1 Full and dashed lines are results for physical equilibrium and nonequilibrium transport respectively obtained with HYDRUS 1D Dots are results obtained with HP1 27 Verification Problem 2 Transient physical nonequilibrium transport of chloride TRANSCL In this second problem we simulate the transport of chloride through a 1 m deep soil profile subject to a transient upper boundary condition given by daily values of precipitation and evaporation over a 300 d period Physical nonequilibrium i e the presence of immobile water in the soil profile was considered in this problem Note again that all transport calculations were done by HYDRUS 1D which means that the test again applies only to the transfer of concentrations and water contents between HYDRUS 1D and PHREEQC We used parameters of the soil hydraulic properties typical for a loamy soil 8 0 078 cm cm 0 0 43 cm cm 0 036 cm n 1 56 and K 24 96 cm d from Carsel and Parish 1988 Solute transport parameters were as follows a dispersivity
78. of three sequential first order decaying contaminants is simulated for transient flow with both HYDRUS 1D and HP1 section 4 1 3 4 1 1 Physical equilibrium and nonequilibrium transport of chloride In this first section we simply test if the transport of an inert solute Cl is correctly simulated using HP1 for the following conditions i steady state flow with no immobile water ii steady state flow with immobile water and iii transient water flow Since this example does 26 not consider any geochemical reactions all transport processes in the coupled model HP1 are simulated only with the HYDRUS 1D module However since solute concentration and water content values still pass between the two modules we will evaluate in this example if the transfer of information between the two components of the coupled model is done correctly Verification Problem 1 Steady state physical non equilibrium transport of chloride EQCL NEQCL This problem simulates the transport of chloride i e a geochemically inert tracer during saturated steady state flow in a 20 cm long soil core The saturated hydraulic conductivity K is 1 cm d and the saturated water content is 0 5 cm cm The following solute parameters were used a dispersivity Dz of 8 cm for both the equilibrium and nonequilibrium cases and for the latter case an immobile water content Om of 0 1 cm cm and a first order exchange coefficient of 0 01 d Simulation results
79. ogistic growth function with the assumption of an exponential root distribution with depth 2 1 2 Initial and boundary conditions for water flow To solve Eq 2 1 initial and boundary conditions must be specified Initial conditions can be defined in terms of pressure heads or water contents Possible system independent boundary conditions are time series of pressure heads or soil water fluxes at the soil surface and or the bottom of the soil profile and a zero gradient free drainage pressure head boundary condition at the bottom of the soil profile In addition system dependent boundary conditions that depend on the status of the system are also available When atmospheric conditions precipitation evaporation and transpiration defining the potential water flux across the top boundary are specified the actual water flux across this boundary depends also on the soil moisture conditions When the potential surface flux precipitation is larger than the infiltration capacity any excess water on the soil surface is either assumed to be immediately removed by surface runoff or is permitted to build up on the soil surface At the bottom of the soil profile the following boundary conditions can be implemented i a seepage face boundary condition that assumes a zero flux when the bottom of the soil profile is unsaturated and a zero pressure head when it is saturated ii a tile drain boundary condition that approximates flow to horizontal sub
80. om and surface of the soil profile moles T L respectively M Mj are total amounts of the jth master species in the flow region at time and at the beginning of the simulation moles L respectively and vj is the charge of the jth master species All other symbols were defined in section 2 M and Mj in Eq 5 2 consist of the total amount of the jth master species in the aqueous phase 0C see Eq 2 28 and the sum of the amounts adsorbed on six cation exchange sites Since an accurate simulation should conserve mass the difference between Eqs 5 1 and 5 2 defined as the absolute error im nek et al 1998 should be small Ql and QZ were calculated using information from the soluteN out output files generated by HYDRUSID while Mi Mj were calculated using values from the output file from PHREEQC defined by the keyword SELECTED OUTPUT every 15 days Figure 5 2 compares Eqs 5 1 and 5 2 for Cl Na Ca and Cd for the 10 year simulation Since differences between values obtained with Eqs 5 1 and 5 2 are very small absolute errors are small as well and no mass balance errors are expected The accuracy of the numerical solution can be evaluated also by using the relative mass balance error amp Eq 6 52 in Sim nek et al 1998 ACI mol m F 62 1 2x10 1 5x10 j l i 1 2x10 8 0x10 ir 90x10 0 i 3 i 4 NS i aor p E i dli oy a j T jJ D Al i fy 1
81. on or are equal to one Note that in PHREEQC it is also possible to express ion exchange reactions using the Gapon convention Gapon 1933 see also Appelo and Postma 1999 2 4 5 Heterogeneous mineral dissolution precipitation For equilibrium precipitation dissolution reactions of a mineral Eq 2 14 is written as Ds A 2 20 where i 1 N N is the number of minerals and is the formula of the mineral while the superscript p refers to pure phases minerals Note that in the database of PHREEQC dissolution precipitation can be written in terms of any of the aqueous species These reactions can always be transformed to a canonical form of Eq 2 20 For equilibrium conditions the mass action equation is Nin p K J ore y 2 21 ja since the activity of a pure phase mineral is assumed to be 1 13 2 4 4 Kinetic reactions Both homogeneous and heterogeneous reactions can be treated as kinetic reactions Homogeneous reactions define the production or consumption of a master species from other master species in the aqueous phase N m Svi 47 0 222 j l where i 1 Nix Nx is the number of kinetically controlled homogeneous reactions in the aqueous phase the superscript k indicates the homogeneous kinetic reaction and Vii are the stoichiometric coefficients involved in the ith homogeneous kinetic reaction The rate equation itself can be of any form and be dependent upon the total concentrations
82. otal concentration in the sink term C p DV Cn 2 29 and assuming that the diffusion coefficients D and D are all equal D i e species independent diffusion and that concentrations of the sink term c and c are equal to each other then the transport equation for the master species becomes 06C QUY oic i aid m Rr MT 2 30 Ot Ox Ox Ox Be aes Note that it is possible to introduce heterogeneous equilibrium reactions in Eq 2 30 in a similar way and thus to define transport equations for the total concentration of the master species using only kinetic reactions as source sink terms Lichtner 1996 Mayer 1999 However this is not done here since the solution method for solving the reactive transport equations is based on a sequential non iterative approach Steefel and MacQuarrie 1996 This means that the transport equation 2 30 is solved without the reaction term Ro whereas the mass action equations 2 16 2 19 and 2 21 are solved sequentially see section 2 6 Note that when c is not equal to zero a component will be taken up by the roots but not a species When c is zero no solute is removed from the soil solution due to the root uptake It is however still possible to define some active uptake mechanism in PHREEQC i e one that is independent of water uptake contained in the R term of Eq 2 30 2 6 Coupling procedures Reactive transport systems as defined in the above sections involv
83. ous reactions biogeochemical reactions 1 2 2 Limitations Specific limitations of the PHREEQC model for various geochemical calculations are discussed in detail by Parkhurst and Appelo 1999 p 4 6 and are not further discussed here Of these limitations those related to flow and transport modelling are no longer of concern here since HYDRUS 1D is used to simulate transient water flow and solute transport One possible limitation involves the invoked method of coupling HYDRUS 1ID and PHREEQC i e a non iterative sequential coupling method SNIA This method is still being discussed in the literature e g Steefel and MacQuarrie 1996 Mayer 1999 with some authors claiming that mass balance errors may occur when this coupling procedure is used However by using appropriate time steps accurate results can be obtained as we will show with examples later in this manual see also Mayer 1999 for some guidelines In addition we believe that uncertainty in the assumed processes and its parameters likely will contribute much more to uncertainty in the model simulations than possible limited numerical errors caused by the coupling procedure This manual documents Version 1 0 of the coupled HP1 model The following features of PHREEQC are not yet operational in the current version surface complexation solid solutions and redox reactions Diffusion and advection of components in the gas phase are also not considered We further do not a
84. overning flow and transport equations were solved numerically using Galerkin type linear finite element schemes To test the accuracy of the coupling procedures implemented in HP1 simulation results were compared with 1 HYDRUS 1D for transport problems of multiple components subject to sequential first order decay ii PHREEQC for steady state flow conditions and iii calculations obtained from an independent geochemical transport model CRUNCH for several relatively complex problems Nine verification examples of increasing complexity are described in this report This report serves as both a user manual and reference document Detailed instructions for input data preparation and interpretation of output data are given in the manuals of the original HYDRUS 1D and PHREEQC codes The graphical user interfaces of both HYDRUS 1D and PHREEQC can be used for easy input data preparation and output display in the MS Windows environment KEYWORDS Biogeochemical model variably saturated water flow multicomponent solute transport vadose zone dissolution precipitation cation exchange aqueous complexation benchmark calculations The software has been verified against selected test cases However no warranty is given that the program is completely error free If you do encounter problems with the code find errors or have suggestions for improvement please contact one of the authors at Diederik Jacques Tel 32 14 333209 Fax 32 14 3235
85. raulic functions These additional features are described in sections 2 6 2 4 and 2 5 of Sim nek et al 1998 respectively The sink term S in 2 1 is defined as the volume of water extracted from the soil by the roots The potential root water uptake rate S x is obtained by multiplying a normalized water uptake distribution b x L with the potential transpiration T gt LT The b x function is obtained from the root distribution with depth whereas T depends on climate conditions and vegetation Both b x and T are input to the HYDRUS 1D model The actual root water uptake rate S x is obtained by multiplying S x with a root water stress response function e g Feddes et al 1978 van Genuchten 1987 to account for a possible reduction in root water uptake due to water and salinity stress conditions in the soil profile Soil water uptake reduction due to the salinity stress can be included using an osmotic head reduction function that can be either additive or multiplicative to water stress The actual water uptake distribution is then of the form S h hg x a h h DT 2 4 where SS is the root water uptake as a function of the pressure head h related to water stress the osmotic head hg L related to salinity stress and depth x related to the root spatial distribution while a h hg defines the reduction in root water uptake due to the water and salinity stress Root growth can be included using the Verhulst Pearl l
86. re 4 6 Small differences between HP1 and PHREEQC exist when a spatial discretization of Ax 0 002 m is used Figure 4 7 If Ax is decreased to 0 0005 m no significant differences between the two simulations were present at the end of the column and in the outflow concentrations We conclude that the keyword EXCHANGE is correctly coupled in HP1 Pore water concentration mol Concentration on exchanger mol l 0 01 0 008 0 006 0 004 0 002 0 01 0 008 0 006 0 004 0 002 2 Depth cm Depth cm Pore water concentration mol l Concentration on exchanger mol l 41 HP1 Ax 0 002 m TM PHREEQC Ax 0 002 m e PHREEQC Ax 0 0005 m 0 0008 0 0006 0 0004 0 0002 0 002 0 0016 0 0012 0 0008 4 CdX 0 0004 Depth cm Figure 4 6 Distribution of pH dissolved Ca Cl Zn and Cd concentrations exchangeable Ca Na Zn and Cd concentrations versus depth after three days of infiltration for Verification Problem 6 42 6 HP1 Ax 0 002 m 5 3 HP1 Ax 0 0005 m disces PHREEQC Ax 0 002 m I MEL LL 21 0 PHREEQC Ax 0 0005 m 3 2 0 3 6 9 12 15 Time days 0 01 S 0 001 Zn 0 008 0 0008 o o E E Y E 0 006 5 0 0006 1 wo 0 004 Ca 0 0004 Cd o o 5 5 0 002 0 0002 0 0 0 3 6 29 12 15 0 3 6 9 12 15 Time days Time days Figure 4 7 pH Ca Cl Zn and Cd concentrations in the efflu
87. re 4 8 Distribution of pH Si Al amorphous SiO and Gibbsite versus depth after 150 days of infiltration for Verification Problem 7 Full line HP1 dashed line PHREEQC 44 4 5 More complicated verification problems of HPI model Examples presented in this chapter compare the combined HP1 model against a different computer program CRUNCH that also simulates multicomponent reactive transport in porous media with the limitation that only steady state flow can be invoked CRUNCH is based on the GIMRT OS3D package Steefel and Yabuski 1996 Steefel 2000 The geochemical reactions and transport in CRUNCH are coupled in one of two ways 1 a global implicit approach GIMRT that simultaneously solves the transport and reaction equations or 2 SNIA GIMRT generally leads to smaller numerical errors A comparison between HP1 and CRUNCH GIMRT allows one to assess numerical discretization errors of the SNIA coupling as a function of the maximum Courant number Cr 4 3 1 Heavy metal transport in a medium with a pH dependent cation exchange complex Verification Problem 8 MCATEXCH 4 3 1 1 Problem definition and governing chemical reactions Cation adsorption to negatively charged soil solid phases can greatly affect the migration of cations in soils In HYDRUS ID equilibrium isotherms such as the linear Freundlich or Langmuir isotherms describe the adsorption desorption of cations the use of such isotherms assumes that the adsorption of
88. re not considered The general solute transport equations as given by 3 1 and 3 2 in Sim nek et al 1998 for conditions described above reduce to 00c Q6 qe Ot Le i z or Sc R 2 5 where i 1 N is the aqueous species number is the total number of aqueous species c is the aqueous concentration of the ith species ML q is the volumetric flux density LT S the sink term in the water flow equation Eq 2 1 c the concentration of the sink term ML D is the dispersion coefficient in the liquid phase LT and R is the general source sink term due to geochemical reactions ML T This sink source term contains heterogeneous equilibrium reactions and homogeneous and heterogeneous kinetic reactions see section 2 4 Physical nonequilibrium solute transport is modeled using a two region model that assumes that the liquid phase can be divided into a mobile flowing region Om ET and an immobile stagnant region O m L L Solutes are exchanged between the mobile and immobile regions by means of a first order exchange process The mathematical formulation of this nonequilibrium model is given by van Genuchten and Wierenga 1976 m Oc Oqc am n SU 20 0 a um Sc Rs 2 6 OO mCi im in R t Q Cim C im i im where Cim and Ciim are concentrations of the ith aqueous species in the mobile and immobile regions ML respectively Rim and Rjim are the source sink term
89. rical Freundlich coefficient The contaminant Cont is additionally assumed to be subject to first order degradation R u C 4 2 w where R is the decay rate M L T and 44 is the first order rate constant for solutes in the liquid phase T Modelling nonlinear Freundlich adsorption with PHREEQC To model instantaneous adsorption using the Freundlich adsorption isotherm we rewrote Eq 4 1 in terms of the amount adsorbed per unit volume of water S K pk 4 3 where S is the adsorbed concentration per unit volume of water M L and p is the bulk density M L Equation 4 3 corresponds with the following mass action equation SorCont Lu er rena ee 44 Sor Cont l i s SorCont K Sor Cont 4 4b where K4 is the adsorption constant in mass per unit volume of water M C E T 721 In PHREEQC this sorption reaction is modelled as a surface complexation reaction The amount of adsorption sites is chosen very large so that Sor in Eq 4 4b does not change significantly when the amount of adsorbed species SorCont remains small An outline of the PHREEQC input file is given in Box 4 1 In this example the new solution master species 29 Box 4 1 PHREEQC input for the Freundlich adsorption isotherm 1 SOLUTION MASTER SPECIES 2 Con Core 10 0 Come 3 50 3 SOLUTION SPECIES 4 Conti cmo TY CT E Ort 5 URE ACH mMAS TER SEECITES 6 Sore Soe 7 SUREACH ROP HE LES
90. rms R and R include both rapid and slow reactions involving the given species However some reactions can be so fast that the rate of the reaction is controlled by the rate of transport of the species to or from the site of the reaction rather than by the reaction itself For these reactions equilibrium can be assumed Consequently the sink source term of the geochemical reactions can be divided as R Rerom 4 genter phinhom pkin het J J J E 4 2 25 R gu Ree qium 4 Ree where the superscripts eq kin hom and het refer to local equilibrium reactions kinetic reactions homogeneous reactions and heterogeneous reactions respectively A system involving both local equilibrium and kinetic reactions is in a state of local partial equilibrium Lichter 1996 A system of N N transport equations can be reduced to the number of primary species The R 4 can be expressed in terms of the reaction stoichiometry defined by Eq 2 15 J Ns Rowe Ly y Re 2 26 i l Substituting Eq 2 26 in Eq 2 23 and replacing R by Eq 2 24 allows the global transport equation for the jth master species to be written as obc A 00c 3 9c c at UAM cue mI cea Nu 5 Oc Vc Oqc x DL J 2 ox Sc LY WS R i l 2 27 where Ro is the term that includes all other heterogeneous equilibrium and kinetic reactions Defining C as the total concentration of a master species Ns Q eee Si vce 2 28 i l 15 and C the t
91. s due to geochemical reactions in the mobile and immobile regions ML T respectively and is the mass transfer coefficient for the ith aqueous species T Note that HYDRUS 1D can also consider chemical nonequilibrium transport kinetic adsorption desorption reactions However we strongly suggest not to use this option of HYDRUS 1D and to simulate chemical nonequilibrium reactions using options in the PHREEQC module The dispersion coefficient is given by 0D D la 0D T 2 7 where D wis the molecular diffusion of the ith aqueous species in free water L T D is the longitudinal dispersivity L and Tw is a tortuosity factor in the liquid phase that is related to the water content as follows Millington and Quirk 1961 eg w 0 S T 2 8 The dependence of the diffusion parameter D on temperature can be described using the the Arrhenius equation see section 3 4 of im nek et al 1998 2 2 2 Initial and boundary conditions for the solute transport equation To solve Eqs 2 5 or 2 6 the initial and boundary conditions must be specified Initial total aqueous concentrations of all aqueous species as a function of depth at time zero in both the mobile and immobile regions must be defined Concentrations of adsorbed secondary or precipitated species must also be specified at time zero when kinetic adsorption and precipitation dissolution reactions are considered Possible boundary conditions inc
92. s of 0 2 cm The flow velocity is 2 cm d and the dispersivity is 2 cm Simulations were performed for 15 days The maximum time step used in HP1 was 0 015 d Distribution versus depth after three days and outflow concentrations during 15 days for selected output variables ions and pH are shown in Figure 4 6 and Figure 4 7 respectively for simulations carried out with PHREEQC and HP1 Only small differences between the two simulations were present in the concentration profiles Figure 4 6 Deviations may be due to increased numerical dispersion in PHREEQC as noted in the Cl concentration profile As discussed in the previous section simulations with PHREEQC showed a larger dispersion compared to the simulations with HYDRUS 1D and HP1 if the same spatial discretization was used e g Figure 4 3a and Figure 4 4a If the spatial discretization in PHREEQC was Table 4 4 nitial and inflow concentrations for Verification Problem 6 mmol I Initial pore water Initial concentrations of Inflow concentrations composition exchangeable cations Al 0 5 0 92 0 1 Ca 1107 2 88 107 5 Cd 0 09 0 17 0 K 2 1 06 0 Mg 0 75 1 36 1 Na 6 0 62 0 Pb 0 1 0 34 0 Zn 0 25 0 76 0 Br 11 3 7 Cl 1107 10 pH 5 5 2 9 Calculated in equilibrium with the initial pore water composition Br is used to impose a charge balance at pH of 5 5 40 reduced to 0 0005 m ie 160 cells PHREEQC converged towards the numerical simulations obtained with HP1 Figu
93. scretization scheme displayed significant numerical dispersion as compared to the HYDRUS 1D solution Increasing the number of cells in PHREEQC and thus decreasing the Peclet number produced less numerical dispersion e g P 2 and P 3 However the simulation time became then very large because of the large number of nodal points in the discretization scheme Results of simulations using HP1 are compared with HYDRUS I1D in Figure 4 3b While the same spatial discretization was used for all simulations the maximum time step was decreased HP1 simulation with a maximum time step of 1 d corresponding to Courant numbers larger than one provided less accurate results as compared to HYDRUS 1D However the agreement gradually improved with lower maximum time steps and Courant numbers with the lowest two Cr values providing almost identical results HP1 results were found to be better than those of PHREEQC for the same Peclet and Courant numbers compare HP 2 and P 1 in Figure 4 3a To obtain a similar degree of accuracy for both HP1 and PHREEQC at least 200 cells were needed for the PHREEQC simulation Figure 4 3a compared to 100 for HPI Also the computational time for HP1 was significantly smaller as compared to PHREEQC 31 When the accuracy of the HP1 code is evaluated based on the simulation with HYDRUS 1D HP1 will need smaller time steps up to 25 times smaller with the same spatial discretization Note that a similar analysis can be done whe
94. sh 1988 Developing joint probability distributions of soil water retention characteristics Water Res Res 24 755 769 Casey F and J Sim nek 2001 Inverse analyses of transport of chlorinated hydrocarbons subject to sequential transformation reactions J Envir Qual 30 1354 1360 Casey F X M G L Larsen H Hakk and J Sim nek 2003 Fate and transport of 17p Estradiol in soil water systems Environmental Science and Technology 37 11 2400 2409 Casey F X M G L Larsen H Hakk and J Sim nek 2004 Fate and transport of testosterone in agriculturally significant soils Environ Sci Technol 38 3 790 798 Chung S O and R Horton 1987 Soil heat and water flow with a partial surface mulch Water Resour Res 23 2175 2186 Dane J H and G C Topp 2002 Methods of soil analysis Part 4 Physical systems SSSA Book Series no 5 SSSA Inc Madison Wisconsin USA 1692 pp de Marsily G 1986 Quantitative Hydrogeology Academic Press London de Vries D A The thermal properties of soils In Physics of plant environment Ed R W van Wijk North Holland Amsterdam pp 210 235 Gaines G L and H C Thomas 1953 Adsorption studies on clay minerals II A formulation of the thermodynamics of exchange adsorption J Chem Phys 21 714 718 Hopmans J W J Sim nek N Romano and W Durner 2002 Inverse methods In Methods of soil analysis Part 4 Co eds Dane and Topp Book Series no 5 S
95. species in ith secondary species stoichiometric coefficient of jth master species in ith secondary exchange species stoichiometric coefficient of j th master exchanger in i th secondary exchange species stoichiometric coefficient of jth master species in ith aqueous secondary species stoichiometric coefficient of jth master species in ith homogeneous kinetic reaction stoichiometric coefficient of jth master species in ith mineral bulk density tortuosity factor in liquid phase mass transfer coefficient for ith species 77 M M M L T T K equivalents LT LT L L L H LY L M L L L L L L L E L L EE E E EE MLT K MLT K T T Qs performance index 78 79 APPENDIX B RELEASED VERSIONS AND BENCHMARKING B 1 Version 1 0 Released date November 2004 Status beta version controlled distribution Versions used of the original codes HYDRUS 1D Version Version 2 0 July 1998 PHREEQC Version 2 4 e Manual Jacques D and J Sim nek 2005 User manual of the Multicomponent variably saturated transport model HP1 Version 1 0 Description Verification and Examples SCK CEN Mol Belgium BLG 998 79 Overview of verification problems Steady state physical non equilibrium transport of chloride Verification Problem 1 HYDRUS 1D EQCL and NEQCL Transie
96. species in solution The activity coefficients are defined with the Davies equation or the extended Debye H ckel equation Langmuir 1997 Parkhurst and Appelo 1999 2 4 2 Heterogeneous ion exchange reactions In PHREEQC the ion exchange sites are defined by exchange primary or master species X and ion exchange reactions are defined as half reactions For the Gaines and Thomas convention Gaines and Thomas 1953 the half reaction is written as N 2 EVIL AT ede 2 17 jal 12 where j 1 Ny Nx is the number of master exchangers i 1 Nse Nse is the number of the secondary exchange species and the superscript e refers to exchange reactions Let the activity of an exchange species be defined as a Y B 2 18 where a is the activity y is the activity coefficient of the 7 th exchange species and 5 is the equivalent fraction of the ith exchange species on the j th exchanger defined as b n T where b is the number of equivalents of exchanger je occupied by the i th exchange species n are the moles of the i th exchange species on exchanger je and T is the total number of exchange sites for the j th exchanger in equivalents Then the mass action equation for equilibrium conditions can be written as N Ki n amp MOT o8 2 19 j l The activity coefficients for the exchange species are calculated with the WATEQ Debye H ckel equation Treusdell and Jones 1974 the Davies equati
97. surface drains using selected analytical solutions and iii a deep drainage boundary condition that uses a functional relationship that relates the water table depth with the deep recharge from the soil profile Mathematical descriptions of these boundary conditions can be found in section 2 7 of Sim nek et al 1998 2 2 Solute transport in the vadose zone 2 2 4 The solute transport equation The HYDRUS 1D code allows simulations of the transport of multiple solutes involved in a sequential first order decay chain in three phases liquid solid and gaseous using the physical nonequilibrium advection dispersion equations However many solute transport features of HYDRUS 1D are not used in the coupled HP1 code since they are considered in the PHREEQC module These include interactions between the liquid and solid phases degradation production and the presence of sequential decay chain reactions These interactions are defined in PHREEQC using equilibrium or kinetic reactions In fact solute transport in the HYDRUS module is modelled as the transport of inert tracers i e no interaction with the solid phase and no solute sink terms since reactions are considered in the PHREEQC module Note however that it is still possible to simulate reactive transport with HYDRUS 1D in the coupled model when particular parameters are entered with non zero values In the present coupled model interactions with and thus also diffusion in the gas phase a
98. surface during the summer periods with their lower water contents Figure 5 4c The pH in actuality may have been affected also by soil carbon dioxide concentrations which usually change in response to the biological activity and moisture status of the soil Sim nek and Suarez 1993 This process however was not considered in our simulations here Although upward flow during the summer had almost no effect on the total amount of heavy metals near the soil surface results not shown due to the low mobility of these elements the concentrations of these elements did vary significantly during the various seasons Because of lower water contents concentrations of all aqueous species and elements increased significantly during the summer periods The changing aqueous concentrations in turn also caused changes in the cation exchange equilibrium by promoting monovalent cations to sorb on the cation exchange complex The lower cation exchange capacity defined here as the amount of deprotonated sites during summer lower pH implied the presence of fewer 64 a Water content b Log Cl mmol 1000 cm soil Depth cm i 5 2 1974 1976 1978 1980 1972 1974 1976 1978 1980 1982 Time year Time year c pH d Log Aqueous Cd mmol 1000 cm8 soil Depth cm 3 1 9 1976 1978 1980 1982 1972 1974 1976 1978 1980 1982 Time year Time year 1972 1974 Figure 5 4 Space time plots of a water content b
99. t and simulation parameters as well as Peclet numbers Pe q Ax D with Ax the characteristic lenght of a finite element and Courant numbers Cr q At 0 Ax Steefel and MacQuarrie 1996 are given in Table 4 1 HYDRUS 1D PHREEQC and HP1 simulations with several different combinations of the Peclet and Courant numbers are compared in 30 Table 4 1 Transport and simulation parameters for Verification Problem 3 Parameters HYDRUS PHREEQC HP1 H P 1 P 2 P 3 HPl 1 HPI 2 HP1 3 HP1 4 K cm d I 1 I 1 I 0 cm cm 0 5 0 5 0 5 0 5 0 5 0 5 0 5 0 5 q cem d 1 2 2 2 1 1 1 1 v cm d 2 1 1 1 2 2 2 2 Diem 1I jore jv povceqoceeegpeucpce M lel p g cm 1 5 z 1 5 1 5 1 5 i5 Ki O 5 5 5 i ny 0 8 0 8 0 8 0 8 0 8 0 8 0 8 0 8 Ki 0 15 2 3 7 5 7 5 7 5 7 5 7 5 Log Sor 100 100 100 100 100 100 100 Log_k 99 125 99 125 99 125 99 125 99 125 99 125 99 125 Ax em 1 1 05 025 a ccc Ata d 1 1 0 25 0 08 0 04 Pe 1 po i 05 095 NDS RO Do d Cr 2 1 33 0 66 0 33 2 0 5 0 16 0 08 Figure 4 3 Depth profiles of aqueous Cont concentrations are plotted at 3 different times 250 500 and 750 d Figure 4 3a compares three simulations carried out with PHREEQC the solution method is equivalent to SNIA sequential non iterative approach with an explicit time weighting scheme and HYDRUS 1D the global implicit method with implicit time weighing The concentration front for the P 1 space di
100. t prejudice to any other right in case of granting a patent or registration in the field of intellectual property SCK CEN Studiecentrum voor Kernenergie Centre d Etude de l Energie Nucl aire Stichting van Openbaar Nut Fondation d Utilit Publique Foundation of Public Utility Registered Office Avenue Herrmann Debroux 40 B 1160 Brussel Operational Office Boeretang 200 2400 Mol Belgium TABLE OF CONTENTS LIST OF TABLES k uietorskiecebhecostes vFetkobieP ed pred sieo isa eoe too Die Eob sS sopes esso b ssis S dV Rel FECE PER PERS PUE v I REESE BEG URES eec m vii T dRtFOdHCUHOB ood iie petet RO Dr HAGA DR MP ases cetex bi eo ed a tek b DROP VE E EUREN RE 1 LI HackproNf l ao deti e RE TUM I RENTUR toast NS pU DERE 1 1 2 Features and limitations of the coupled HPI model ecce 3 LXI PFeatur sqaeseNiU ieri ep Up rH Dos Pe M DURER PED ao I seu REEE ERES 3 1 22 2 Limitations mr 3 2 Description of the model dio en Srt e RESO ESO EEONSEPVINM SINON NINE INVE D UTE TRIS 5 2 1 Water flow tn the vadose z0R8 iieeos suse nts i I Eo Me eNS eu eN PE pMR UTR EVER E RM 5 2 1 1 The water flow equation etes nee NR ERE KENNEN RAE o LER enemas 5 2 1 2 Initial and boundary conditions for water flow e eeeeeeee 6 2 2 Solute transport in the vadose zone eese ee sees esee esee eene ette sten sten netus enn 7 2 2 1 The solute transport equation
101. t version of the model All boundary conditions are specified in the input files of HYDRUS 1D Initial chemical conditions are to be defined in the phreeqc in file This includes concentrations of all master species If the initial concentration of a given master species is zero a very small value e g 1E 20 mol l should be given otherwise it is possible that the master species will be neglected during the simulation There is no need to define the initial conditions of the total aqueous concentrations in the input files of HYDRUS 1D Initial conditions for water flow and heat transport are to be defined in profile dat For each node in the HYDRUS ID finite element mesh a SOLUTION must be defined in the pAreeqc in file The node number at the soil surface is 1 and numbering increases with depth Solutions at different depths with the same initial composition can be grouped in the phreegc in input file Since the phreeqc in data file is read prior to the profile dat data file the amount of initial water for each SOLUTION defined in the profile dat file is not known when the initial solutions are initialized Therefore the water content at each node and 23 hence for each solution should be specified for each SOLUTION key word in the phreeqc in file Alternatively the SOLUTION SPREAD key word can be used to define solutions with varying amounts of water and temperatures as a function of depth e The input file phreeqc in must contain th
102. teady state flow conditions and with calculations obtained with an independent geochemical transport model CRUNCH Steefel 2002 for several more complicated problems 1 2 Features and limitations of the coupled HP1 model 1 2 1 Features Any combination of the following features can be described with the HP1 model e One dimensional transient water flow for different boundary conditions including atmospheric conditions precipitation evaporation transpiration e Root water uptake as a sink for water e Root growth e One dimensional transient convective and conductive heat transport for time variable temperatures at the soil surface e One dimensional advective dispersive and diffusive transport of multiple solutes e Effect of temperature on transport parameters thermodynamic constants and rate parameters e Different functional forms for the soil hydraulic properties including hysteresis e Physical non equilibrium solute transport e Physical and chemical spatial heterogeneity e Equilibrium aqueous speciation reactions and kinetically controlled aqueous reactions such as radioactive decay e Multi site cation exchange related to type and amount of minerals present e Equilibrium and kinetic dissolution precipitation of primary and secondary minerals e User defined kinetic reactions e Simultaneous presence of different reactions sequential and parallel kinetic reactions equilibrium and kinetic reactions homogeneous and heterogene
103. the different solutes is a consecutive chain reaction in which the solutes are sequentially transformed along the chain by means of first order reactions which hence depend only on the concentration of the first solute van Genuchten 1985 These chain reactions can be used to describe the degradation of pesticides e g Wagenet and Hutson 1987 or chlorinated volatile organic compounds e g Schaerlaekens et al 1999 Casey and Sim nek 2001 and consecutive decay chains involving radionuclides e g Mallants et al 2003 endocrine disrupting chemicals Casey et al 2003 2004 and many other chemicals No other interactions between different species or components are currently considered in the HYDRUS 1D model The only attempts to include geochemistry into HYDRUS 1D related models were those by Sim nek and Suarez 1994 1996 and Suarez and Sim nek 1997 resulting in UNSATCHEM 2D which was based on the SWMS 2D code Sim nek et al 1992 a two dimensional precursor of HYDRUS 2D and UNSATCHEM These two models considered interactions only between major ions whose possible components Ca Mg Na K SO Cl alkalinity and CO2 and geochemical reactions speciation cation exchange and precipitation dissolution of minerals and their kinetics are predefined and hence can not be changed by the user Although these two codes can be applied to a wide range of important problems such as salinization of agricultural soils under irrig
104. tion of the fraction of deprotonated sites 1 H sites and sites with Cd versus depth after 0 3 0 5 and 0 7 y respectively 49 Furthermore at the boundary between the fourth and fifth layer the pH shows a small increase for both models using the sequential non iterative solution approach This increase seems to be an artefact of the implemented time step Reducing Cr to 0 1 in HP1 produced very good agreement with CRUNCH GIMRT with the artefact being no longer present Figure 4 9 Overall differences in results obtained with Courant numbers of 0 1 at least one order of magnitude smaller than typical values for Cr and 0 5 were quite small especially for breakthrough curves with both runs providing acceptable results To illustrate the effect of a decreasing cation exchange capacity with decreasing pH on Cd transport we changed concentrations of the infiltrating water to those found in the 28 50 cm horizon i e with a pH of 8 5 see Table 4 6 Simulations were carried out with CRUNCH GIMRT and HP1 using a Courant number of 0 5 this since no considerable differences were observed for the Cd breakthrough curve in previous calculations As shown in Figure 4 10 the Cd breakthrough is significantly retarded compared to that produced with the low pH infiltration water Finally to evaluate the effect of heavy metals complexation with CI ions we increased Cl concentration in the infiltrating water to 780 umol per kilogram of water and de
105. unstable and may dissolve Dissolution of primary minerals is described with kinetic dissolution rate equations The infiltration of Na Ca rich water and the increase of Al and Si from mineral dissolution will cause the precipitation of secondary minerals The infiltration of a hyperalkaline solution in a 7 4 cm long clay core was simulated for a period of 1 1 year The flow domain was discretized in 100 cells of 7 4 10 cm each The effective porosity of the Opalinus Clay was kept constant at a value of 0 13 during the simulations A constant flux of 2 403 10 m sec 7 58 cm y was applied resulting in a solute transport velocity of 1 85 10 m sec The diffusion coefficient was 0 5 10 m sec and the dispersivity zero The mineralogical composition of Opalinus Clay is given in Table 4 9 Due to the intrusion of the high pH plume primary minerals of the clay will become unstable and gradually dissolve The dissolution rate depends on the chemical composition of the aqueous phase and is described with a transition state based dissolution model Aagaard and Helgeson 1981 Parkhurst and Appelo 1999 2 3 Ra A uo a nr Kom 1 z Q K 4 7 m 0 where Rm is the dissolution rate of mineral m mol sec Am is the initial reactive surface area of mineral m m for HP1 m m for CRUNCH Mn Mmo is a factor accounting for The conceptual geochemical model used here is somewhat simplified compared to the model described by
106. ution processes of the primary minerals were treated as kinetic reactions in HP1 equilibrium was assumed for calcite Eq 4 7 was also used to describe the kinetic precipitation of two secondary minerals see Table 4 9 in CRUNCH These minerals were treated as being in equilibrium with the solution in HP1 Quartz kaolinite and illite were not allowed to precipitate in either model Table 4 10 Parameters for the kinetic dissolution reactions Eq 4 7 Mineral Log ko n Ao mol m sec m m Kaolinite 10 88 0 54 1 00 107 Illite 13 91 0 22 0 33 107 Quartz 10 20 0 30 2 00 10 Calcite 1 00 0 00 2 00 10 Dolomite 7 70 0 00 2 00 10 Gypsum 8 00 0 00 2 00 10 Hydrotalcite 1 00 0 00 1 00 10 Sepiolite 1 00 0 00 1 00 10 Data from Table 5 Part III Adler 2001 Minerals are assumed to be in equilibrium in HP1 High dissolution constants are used in CRUNCH 52 Table 4 11 Overview of aqueous equilibrium reactions and cation exchange half reactions with corresponding equilibrium constant Nr Reaction Log K Aqueous speciation reactions 1 H O OH H 14 0 2 HCO H CO H O 6 3447 3 HCO CO H 10 3288 4 Al H5O AIOH H 4 9571 5 AI 2 H O Al OH 2 H 10 5945 6 Al 2 HO A10 4 H 22 8833 7 Al 2H O HAIO 3 H 16 4329 8 AI SO AISO 3 01 9 AI Na 2 H O NaAlO AH 23 6266 10 AI 2 SOT AI
107. ve saturation Sp potential root water uptake 1 a z z ZU at 47 og adsorbed concentration adsorbed concentration per unit volume of water time temperature average temperature at soil surface total number of exchange sites for j th master exchanger potential transpiration average pore water velocity volume of mineral m initial volume of mineral m spatial coordinate chemical formula of jth master exchanger dimensionless water stress response function parameter in the water retention curve angle between flow direction and the vertical axis equivalent fraction of i th exchange species on jth exchanger thermal dispersivity absolute error in the solute mass balance relative error in the solute mass balance of the jth master species activity coefficient of ith secondary exchange species activity coefficient of ith aqueous secondary species activity coefficient of j th master exchanger activity coefficient of jth master species volumetric water content immobile water content mobile water content volumetric fraction of the solid phase volumetric fraction of the organic matter residual water content saturated water content volumetric fraction of the air phase apparent thermal conductivity of the soil thermal conductivity of the soil first order rate constant for solutes in the liquid phase first order rate constants for decay chain solutes in the liquid phase stoichiometric coefficient of jth master

User Manual of the Multicomponent Variably - PC

Contents

Download Pdf Manuals

Related Search

Related Contents