Effect of Density on the Migration of Liquid Radioactive Waste in a Non-Horizontal Freshwater Aquifer

Radioactive Waste in a Non-Horizontal

Freshwater Aquifer:

Case Study for the Krasnoyarsk-26 Nuclear Waste Disposal Site

MSc Graduation Research Project Maartje Boon

Faculty of Geology Hydrogeology Group Moscow State University

Effect of Density on the Migration of Liquid Radioactive Waste in a Non-Horizontal Freshwater Aquifer:

Case Study for the Krasnoyarsk-26 Nuclear Waste Disposal Site

Maartje Boon 2005/2006 Abstract

Deep-well injection is an alternative method of liquid radioactive waste (LRW) disposal within the Russian Federation. Hydrogeological modeling is used to make predictive calculations of the long-term migration of the LRW to assure the future safety of this method. LRW can be represented by two components: A dense neutral component 1, representing most of the waste and a light sorptive component 2 representing the radionuclides. Dense component 1, often, causes a density contrast, resulting in a density-driven convection with the more dense solutions sinking down the dip ofthe aquifer. The effect of this density driven convection is different for both components.

In this study the effect of density-driven convection on the long-term migration of both components of the LRW was analyzed for the non-horizontal Horizon I of the Krasnoyarsk-26 disposal site. The following tasks were under investigation: 1) A comparative study of the driving forces behind the migration of the two component waste plume in a sloping 2-D aquifer. 2) Analysis of the influence of aquifer vertical anisotropy on dense plume migration. 3) Simulation of post-injection (2 component) waste spreading at Horizon I of the Krasnoyarsk-26 site. 4) Analysis of the main driving forces that control waste spreading for Krasnoyarsk-26. In addition, the possibility to use a similarity solution, developed by Barenblatt to describe water flow in an initially dry aquifer, to describe flow of dense-water in an initially freshwater aquifer was analysed. The Modflow2000+MT3DMS and Seawat2000 codes were used to simulate density-independent and density-dependent flow and transport, respectively.

The modelUng results showed that dense component 1 is transported more towards the deeper part of the aquifer while light sorptive component 2 is transported in the direction ofthe regional convection. Due to sorption the flow velocity of component 2 is much smaller than component 1. The isotropic and vertical anisotropic results were similar. The numerical and analytical resuhs of Barenblatt's problem were similar for the non-flowing/flowing isotropic environment. For anisotropic environments they were different.

Table of Contents

A B S T R A C T 2 T A B L E O F C O N T E N T S 3

1 I N T R O D U C T I O N 5 2 B A C K G R O U N D 7

2.1 LIQUID RADIOACTIVE WASTE ( L R W ) 7 2.2 DRIVING FORCES BEHIND L R W MIGRATION I N SLOPING AQUIFER 8

2.3 SITE DESCRIPTION KRASNOYARSK-26 9 2.4 GEOLOGY OF KRASNOYARSK-26 1 1 2.5 WASTE INJECTED AT KRASNOYARSK-26 13 2.6 PREVIOUS MODELING ATTEMPTS FOR HORIZON I OF KRASNOYARSK-26 13

3 ( N O N ) - C O N S T A N T D E N S I T Y - F L O W A N D T R A N S P O R T M O D E L

F O R M U L A T I O N 1 4 3.1 FLOW EQUATION 14

3.1.1 Constant Density Flow Equation (Modflow2000) 16 3.1.2 Non-Constant Density Flow Equation (Seawat2000) 17

3.2 TRANSPORT EQUATION 19 3.3 COUPLING OF FLOW AND TRANSPORT 20

4 D E N S E F L O W I N H O R I Z O N T A L A Q U I F E R 2 2 4.1 PROBLEM DESCRIPTION 22 4.2 MODELING DESCRIPTION 29 4.3 RESULTS 3 1 4.4 DISCUSSION 32 4.5 CONCLUSIONS 33 5 E F F E C T O F D E N S I T Y O N T H E M I G R A T I O N O F L R W I N A S L O P I N G A Q U I F E R 3 4 5.1 PROBLEM DESCRIPTION AND METHODS 34

5.1.1 Problem Description ^4

5.1.2 Model and Modeling Description 35

5.1.3 Importance of Grid. -^7 5.1.4 Schemes

5.1.5 Center of Mass

5.2 RESULTS AND DISCUSSION 39

5.2.1 Schemes 39 5.2.2 Importance of Grid. 42

5.2.3 Cojnparison ofthe Different Driving Forces behind LRW Migration 45

5.3 CONCLUSIONS 56 6 C A S E S T U D Y F O R K R A S N O Y A R S K - 2 6 S I T E 58

6.1 PROBLEM DESCRIPTION AND METHODS 58

6.1.1 Problem Description ^8 6.1.2 Model Description

6.1.3 Modeling 65 6.2 RESULTS 66

6.2.1 Sinmlation of Injection Period 66 6.2.2 Initial Condition Set 1 of Calibration Scenario 1 66

1 Introduction

Tlie disposal of liquid radioactive waste (LRW) by means of deep well injection is applied at thi'ee locations within the Russian Federation: Tomsk-7, Dimitrovgrad and Krasnoyarsk-26. This method is based on the existence of a permeable layer bounded by low permeable layers. The waste is injected into the permeable layer while the low permeable layers should prevent its connection with the environment. Upon injection the waste migrates through the layer to the discharge zone where it returns into the environment, while simultaneously the radionuclides are decaying. By the time the waste reaches the discharge zone, the radionuclides should have decayed towards harmless levels.

Up till now the use of deep well injection to dispose of LRW at those three sites seems to be successfiil. To assure the future safety of this method a clear picture of the processes taking place in the subsurface needs to be formed and predictive calculations ofthe long-term migration of the waste need to be carried out. For these purposes hydrogeological modeling is used. Accurate modelling can only be done when all key parameters influencing the plume migration are taken into account. The most important factor influencing the migration of the plume is the regional groundwater flow within the injection formation. Secondly, the natural heterogeneity ofthe formation plays an important role and, thirdly, characteristics of the waste itself can influence the migration of waste.

LRW consists of sorptive radionuclides, which are only present in trace amounts, salts, metal ions, acids, alkalies, detergents, organic compounds and many other substances. Due to its high salt content LRW has a relatively high density compared to freshwater. This resuhs, often, in a density contrast between the injected LRW and the naturally present aquifer-water which can be up to 350 g/L (Rybal'chenko, A . I . et a l , 1998). Density contrasts induce flow: more dense solutions sink down the dip of the aquifer due to buoyancy forces. Therefore, density contrasts together with a non-horizontal aquifer form additional driving forces behind the transport direction ofthe LRW. Depending on the geometry of the aquifer bottom, this density driven flow can cause the transport direction of the LRW to deviate from the direction of regional-flow present in the aquifer (Malkovsky, V . I . et al., 1999).

Another characteristic ofthe waste that can effect its migration is the sorptivity ofthe radionuclides. Upon penetration of the LRW into the aquifer radionuclides are sorbed onto the host rocks. This sorption process causes retardation, which implies that the radionuclides are transported with a lower velocity than its surrounding water (Rybal'chenko, A . I . et a l , 1996 (2)).

Based on density and sorption LRW can be represented by two components: A highly dense non-sorptive component 1, representing most part of the waste, and a light sorptive component 2, representing the radionuclides present in trace amounts.

Because of this binary nature the following process is expected to occur during the migration of the LRW in a non-horizontal aquifer: The high dense component 1 causes the flow-field to change. At first this changed flow-field exerts an influence on both components but sorption causes the component 2 to be retarded. Therefore, eventually, the two components will separate and will both be influenced by different flow-fields.

of Krasnoyarsk, Siberia and occupies approximately 6.5 km^. The site is situated within a paleo-erosional depression in crystalline rocks of Precambrian age (gneisses), filled with Jurassic sand-clay sediments reaching 550 m below the ground surface. Three relatively high permeable aquifers divided from each other by low permeable aquitards can be distinguished and from bottom to top they are defined as horizon I , I I and I I I . The bottom topography of the aquifers is non-horizontal with its deepest part north ofthe injection zone. Injection of LRW is carried out in Horizon I and I I . The aquifers used for injection comprise medium-grained sands, weakly cemented sandstones and clay lenses. The presence of clay lenses causes vertical anisotropy i.e. the horizontal hydraulic conductivity is much larger than the vertical hydraulic conductivity (Rybal'chenko, A . I . , et a l , 1998;Compton, K.L. et al., 2000).

In previous modeling tasks for the migration of LRW in Horizon I of the Krasnoyarsk-26 site carried out by the All-Russia Design and Research Institute of Production Engineering (VNIPIPT) and the Institute of Geology of Ore Deposits, Petrography, Mineralogy, and Geochemistry of the Russian Academy of Sciences (IGEM) waste was represented by a single dense component and vertical anisotropy was not taken into account. Their research showed that gravity forces tend to move the dense plume predominantly to the north dovra the dip of Horizon I and settle into the depression north ofthe injection zone rather than moving northeasterly, following the pattern of regional flow as suggested in previous modeling attempts without density taken into account (Malkovsky, V . I . et a l , 1999; Compton, K.L. et a l , 2000). To form a picture of the effect of density differences on the migration of each of the LRW components of the waste plume in Horizon I of the Krasnoyarsk-26 disposal site and to verify i f the former described process of separation of waste components is actually taking place the following tasks were under investigation:

• A comparative study of the driving forces behind the migration of the two component waste plume in a sloping 2-D aquifer.

• Analysis of the influence of aquifer vertical anisotropy on dense plume migration.

• Simulation of post-injection waste spreading in Horizon I of the Krasnoyarsk-26 site.

• Analysis of the main driving forces that control waste spreading in Horizon I of the Krasnoyarsk-26 site.

In addition, the possibility to use a similarity solution, developed by Barenblatt to describe water flow in an initially dry aquifer, to describe flow of dense-water in an initially freshwater aquifer was analysed.

In this study complex chemical processes of water, waste and rock interaction were not taken into account. Focus was placed only on density and sorption effects. This simplified approach can show how density of waste and geometry of aquifers effects on spreading of sorbed and dense component. To make a more realistic prediction the models for geochemical interaction and radioactive decay can be added to this model. For this study the Modflow2000 code in combination with the MT3DMS code were used to simulate density-independent flow and transport. Density-dependent

2 Background

2.1 Liquid Radioactive Waste (LRW)

Liquid radioactive waste consists of radionuclides, which are only present in trace amounts, sahs, metal ions, acids, alkalies, detergents, organic compounds and many other substances. Its specific activity ranges from 3.7x10^ - 3.7x10 Bq/1 (Rybal'chenko, A . I . , et al., 1998). The specific activity is a measure of the amount o f radioactivity per unit amount of substance.

The classification of wastes into low-, intermediate- and high-level waste is based on specific activity and chemical composition. Table 1 shows the characteristics of each type of waste (C&E et a l , 1999).

Type of waste Specific Activity Salt Content Chemical Composition pH Type of waste

(Bq/I) (g/1)

Low-level 3.7x10^-3.7x10^ 1-30 sodium nitrate 6-8

Intermediate-level

3.7xl0'-3.7xl0'" 30-350 nitrates and soluble

complexes of aluminium and silicon

alkaline

High-level 3.7x10^-3.7x10" 250-350 nitrates,

soluble complexes of iron, chromium, manganese,

and nickel, and complex forming reagents

1-3

Table 1 Type of waste and its characteristics (C&E et al., 1999; Compton, K . L . et al., 2000).

Nuclide Co/CoCb' Distribution coefficient,Ki Retardation Desorption

(cm%) factor (R) (%) High-level LRW *Sr 1.0x10^ 1.2-1.3 6-9 21-22 '""Ru 7.0x10' 1.3-1.6 6-7 0.5-12 '^•'Cs 2.0x10' 2.7-2.8 13 0.7-1.3 "^Ce 6.5x10" 1.0-1.6 4-8 12-19 100 1.8-2.3 9-11 8-40 "'-»Pu 1000 1.4-1.6 13-30 0.5-3.0 '''Km 10000 1.1-L2 5-6 18-20 Intermediate-level LR W ^"Sr 2.5x10' 5.5-7.0 26-33 16-31 '"'•Ru 6.7x10' 6.0-10.5 30-50 12-30 '"Cs 6.0x10'' 4.5-6.5 30 11-13 ^^¥u 140 40-95 200-460 10-15 Low-level LRW 5000 35-37 200-300 5.5-9.0 900 7-14 45-85 1.4-4.0 '"Cs 270 60-90 350-450 1.5-2.0 '-•^Ce 80 - 95-115 1.5-4.0 ^^^Pu 10 110-140 500-625 3-6

Table 2 Retentive properties of radionuchdes in typical Russian reprocessmg wastes (Compton, K . L . et a l , 2000).

LRW has a relatively high density compared to freshwater due to its high salt content. As can be seen from Table 1 the sah content can reach values up to 350 g/L.

The radionuclides present in LRW are isotopes of strontium, ruthenium, cesium, cerium, promethium, and a number of other elements (Table 2). One important characteristic of radionuclides is their sorptivity. Radionuchdes are sorbed to the rocks of the aquifer which lowers their flow velocity and, therefore, more time is available for radioactive decay. Table 2 shows the retentive properties of radionuclides in typical Russian reprocessing wastes. (Compton, K.L. et a l , 2000; Rybal'chenko, A . I . , et al., 1998).

Based on density and sorptive properties the waste can be divided into two components: A highly dense component 1, representing most part o f t h e waste and a light sorptive component 2, representing the radionuclides present in trace amounts. Due to this two component nature one, often, speaks of binary waste. The initial concentration of each of the components is much larger than allowed according to the Russian legislation. The allowed value for NaNOs is 0.045 g/1 and the allowed value for the different radionuclides can be determined from column 2 in Table 2.

2.2 Driving Forces behind LRW Migration in Sloping Aquifer

The driving forces behind migration of a LRW plume in a sloping aquifer are regional-, density-dependent-, gravitational-, and thermal convection.

Upon injection waste gets carried away with the naturally present groundwater flow, the driving force behind this type of migration is called regional-convection.

Often a density contrast exists between the injected LRW and the naturally present aquifer-water. The LRW injected in Horizon I of the Krasnoyarsk-26 disposal site has a density of up to 240 gr/1, while the density of the naturally present aquifer-water does not exceed 0.3 gr/1 (Rybal'chenko, A.L, et al., 1998). Due to this density contrast gravity starts to effect the plume migration. The role of gravity is two-fold: First, due to density- dependent dispersion and the force of gravhy the waste wants to go to a low position, the driving force behind this type of migration is called density-dependent convection. The intensity of this process depends on the density contrast and the initial thickness of the dense body present in the aquifer. This process occurs in both horizontal and non-horizontal aquifers. Second, the dipping geometry of the aquifer resuhs in a down-slope movement of the dense waste, also driven by the force of gravity. This driving force is called gravitational convection.

2.3 Site Description Krasnoyarsk-26

Tlie Krasnoyarsk-26 disposal site of the nuclear weapon facility Mining Chemical Combine (MCC) is one of the three locations within the Russian Federation where LRW injection is practiced. The repository is located 60 km north of the city of Krasnoyarsk, Siberia (Fig. la) and 2.5 km from the Yenisei River, one of the largest Siberian rivers that runs from south to north in Western Siberia. The injection site is situated 12 km away from the radiochemical plant of the MCC and occupies an area of approximately 6.5 km^ sun-ounded by an exclusion zone of 52 km^ (Fig. lb). The site was put into operation in 1967 and plans are to keep injecting waste till 2010 (C&E, e t a l , 1999).

t

. D n i l r u v o r a a

lunst.

Fault zone

Boundary of injection zone Boundary of exclusion zone Lines of geological cross section Observation wells

Figure 1 a) Map of major deep well injection facilities in the Russian Federation, b) Site plan Krasnoyarsk-26 disposal site (Compton, K . L . et al., 2000).

Before the disposal site was created, special geological prospecting works and explorations were carried out to determine the sites feasibility and safety for LRW injection and to provide information on the distribution of geophysical properties and hydraulic heads in the area (Compton, K.L. et al., 2000).

This research is concerned with Horizon I (section 2.4). The details about the different wells concerning Horizon I can be seen in Table 3.

Well type Amount Special characteristic Injection well 7 arranged in a linear fashion.

Relief well 8 located at a distance of 1 km from injection wells in direction opposite to natural GW flow.

Observation well

46 inside repository 42 outside repositoiy

Table 3 Details of different wells in Horizon I (C&E, et al., 1999; Rybal'chenko, A.I. et al., 1998).

2.4 Geology of Krasnoyarsk-26

The injection site is situated within a paleo-erosional depression in crystalline rocks of Precambrian age (gneisses), filled with Jurassic sand-clay sediments reaching 550 m below the ground surface. East-West and North-South cross-sections of the injection site can be seen in Figure 2.

Three aquifers divided by aquitards are present within the paleo-erosional depression, fi-om bottom to top Horizon I , I I and I I I respectively. The aquifers and aquitards both consist of sand and clay sediments but the percentage of sand and clay differs. Sand-sublayers are predominantly present in the aquifers while the aqmtards consist mainly out of clay-sublayers. Horizon I I I is not fully saturated and does not occur everywhere which makes this horizon not suhable for waste disposal. Horizon I and I I , situated at depths of 350-500 m and 180-280 m below the surface respectively, are fully saturated and occur over a wider extent but do wedge out on the east. Horizon I and I I are used for waste disposal. From the west the depression is bounded by a fauU zone which runs North-South and is composed of clay. This fauh zone truncates the sedimentary strata both structurally and hydraulically, separating the injection zones from the Yenisey River. The bottom and edges of the depression are formed with gneisses and many-colored overlapping clays (IGEM, 1997; Rybal'chenko, A . I . et al.,

1998; Compton, K.L. et a l , 2000).

This research project concentrates on Horizon I which consists of gravel sands, breccias, non-sorted brecciated rocks with limestone cement and clay lenses. The aquifer-water naturally present in Horizon I has a sodium-hydrocarbonate composition with mineralization up to 0.3 g/1. Recharge in Horizon I occurs 7 km south ofthe injection site and the discharge area is according to [Rybal'chenko, A . I . et a l , 1998] into the Kan River, 12-14 km north of the site.

The initial steady state hydraulic head distribution of Horizon I (Fig. 3) shows that the regional flow is North/North-Easterly. Furthermore, the two iso-head lines of 165 m in the North of Figure 3 imply that there is no driving force for the North- Northeast movement between those lines. Since the wedging out of Horizon I on the east can be treated as an impermeable boundary, the west side of the depression is bounded by the impermeable fauh and it is shown that groundwater flows into the area at the southern boundary, a zone of discharge must be present between those equi-isohead lines. This implies leakage between Horizon I and I I (Malkovsky, V . I . et a l , 1999).

The bottom topography of Horizon I is non-horizontal (Fig. 4). A depression can be seen North/North-West of the injection zone.

i : 1:40000 s3» - 4000- 3000-S-36 170 17» tea i 2000- i ^ J 2 ANift A ; M A V I Ü 1000-JilB A - 2 6 ~ 165 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 5500 Center of injection Fault zone

Observation wells with measured predevelopment head

Figure 3 The steady state hydraulic head distribution m Horizon I.

12000 6000 5000 1 2 0 100 8 0 6 0 - 4 0 - 2 0 0 -20 -40 - 6 0 - 8 0 - 1 0 0 - 1 2 0 - 1 4 0 - 1 6 0 - 1 8 0 - 2 0 0 - 2 2 0 - 2 4 0 - 2 6 0 Elevation (m) Center of injection 3000 4000 5000 6000

2.5 Waste Injected at Krasnoyarsk-26

Exact data about the amoiint and type of waste injected at the Krasnoyarsk-26 site is not available. The data is secret because knowledge of this data would make it possible to calculate the amount of weapons produced at the MCC. [Robinson and Volosov, 1996] of the Green Cross, however, made some estimates of the characteristics of the waste disposed. Their resuhs can be found in [Compton, K.L. et aL, 2000].

At present time Horizon I is used for the injection of intermediate-level LRW. In the period from 1972 to 1993 high level LRW were also injected into Horizon I . Horizon I I is used for injection of low-level LRW, and Horizon I I I is used as a buffer (IGEM,

1997).

As this research is concerned with Horizon I , the following description of wastes w i l l be focused on the waste injected into this horizon.

The high level LRW were injected into Horizon I once to twice a year in batches o f 1000-2000 m l The intermediate level LRW were regularly injected in Horizon I from spring to fall at rates of up to 300 m^ per day (Compton, K.L. et a l , 2000).

According to [Rybalchenko, A . I . et al., 1998] the waste injected into Horizon I has a total sah content up to 240 g/1, a specific activity not exceeding 3.7x10^ Bq/1 and a radionuclide composhion of strontium, cesium, ruthenium, cerium, and several other short-lived nuclides.

2.6 Previous Modeling Attempts for Horizon I of Krasnoyarsk-26

Modeling attempts for the Krasnoyarsk-26 site taking density differences into account have been carried out by the All-Russia Design and Research Institute of Production Engineering (VNIPIPT) and the Institute of Geology of Ore Deposhs, Petrography, Mineralogy, and Geochemistry of the Russian Academy of Sciences (IGEM).

Both, the VNIPIPT and IGEM, assumed an essentially impermeable fault zone and vertically homogeneous aquifers characterized by lateral heterogeneity of thickness and transmissivity. Furthermore, they used a non-reactive, non-decaying one component contaminant to obtain a conservative estimate of plume movement. The IGEM used a multi-layer two-dimensional model, with a term reflecting hydraulic interaction between Horizon I and I I . They assumed an impermeable boundary on the east, based on wedging out of the aquifer. A site specific finite difference algorithm was developed to solve the flow and transport problem (Malkovsky, V . I . et al., 1999). The VNIPIPT used a single-layer, two dimensional (areal) model with a no-flux boundary at every point in the interior of the model, and constant head boundaries on the edges of the modeling domain where appropriate. They used internal Russian codes to solve for flow and contaminant transport (Compton, K.L. et al. 2000).

3 (Non)-Constant Density-Flow and Transport Model Formulation

For this study the Modflow2000 + MT3DMS codes were used to simulate flow and transport respectively, without density taken into account. To simulate density-dependent flow and transport, the Seawat2000 code was used. A description of these three codes and the assumptions and justification for their use in the simulations of this study can be found in appendix A.

In this section the flow equations used in Modflow2000 and Seawat2000 will be derived. Furthermore, the transport equation, which is similar in both the MT3DMS and Seawat2000 codes and the coupling of the flow and transport equations, will be discussed.

3.1 Flow Equation

The flow equation is based on the law of mass conservation. The law of mass conservation states that the mass entering the representative elementary volume (REV) of a porous medium less the mass leaving the REV is equal to the change i n mass storage with time (Fetter, C.W. 1999) (Guo, W. et al., 2002). Mathematically this can be written as:

Eq. 1 where: 8 d d V : differential operator — + — + — da op dy p : fluid density [ML"^]

q : specific discharge vector [LT"^]

p : density of water entering from a source or leaving through a sink [ML' ] q^ : volumetric flow rate per unit volume of aquifer representing sources and

sinks [T-^]

6 : porosity [-] t : time [T]

The first term on the left hand side is the net flux of mass through the faces of the REV and the second term represents the mass entering through sources or leaving through sinks located in the REV. The right hand side describes the time rate of change in the mass stored in the REV over a given period of time and can be expanded with the chain rule (Guo, W. et a l , 2002):

d{p0)_ 80 ^^dp 8t dt dt

Eq. 2

The changes in porosity here considered are restricted to those associated with fluid pressure and can therefore be expressed as:

de__de_8P_

8t ~ 8P 8t

Under isothermal conditions the equation of state for fluid density is: P = / ( A C ) where: Eq.4 P C

fluid pore pressure [ML'^T'^j

-3n

: solute concentration [ML' ]

Differentiation of equation 4 with respect to time gives:

dp _dp_dP_ d^dC_ 'dt~'dPdt dC dt

Eq. 5

After substituting equation 3 and 5 into equation 2, the following expression is obtained for the right side of equation 1 (Guo, W. et al., 2002):

d(p0) d9dP ^dpdP^^dpdC

— M L — p h t ) \-U

dt dP dt dP dt dC dt _{Eq. 6}

The specific discharge on the left hand side of equation 1 is given by Darcy's law:

(dP dz\ \da da) qp = kp fdP dz]

A

: principle directions of the hydraulic conductivity tensor : i component of specific discharge vector [LT'^j

: permeability in i direction [L^] : dynamic viscosity [ML'^T'^J : gi-avitational constant [ L T ]

da f dP dz^^ \da da + - ' p k / dp dP dz — + pg— dp dp + • JJ dy PK dP dz dy dy ^ dP dt dP dt dC dt + pq. Eq. 8

3.1.1 Constant Density Flow Equation (Modflow2000)

For constant density flow, the general flow equation given in equation 8 can be simplified:

The terms dp/dC anddp/dP on the right hand side of equation 8 are zero when density is constant and, therefore, cancel out. Furthermore, density can be taken out of the space derivatives on the left side of equation 8.

The pressure P in equation 8 can be expressed as:

P = pgih-z)

where:

h : hydraulic head [L] z : elevation [L]

Eq.9

For a constant density case the only variables in equation 9 are the hydraulic head and the elevation. In this case the pressure derivatives are given by:

dP dh dz _{= pg pg} dl dl dl

Eq. 10

where a,P,y,t can be substituted for /.

Since the elevation does not change with time, the last term of equation 10 cancels for the time derivative.

Using the former and dividing both sides by density the following expression of the flow equation is obtained:

The term dO

d{h-z)

on the right hand side represents the volume of water released

fi-om storage in a umt volume of a confined elastic aquifer per umt change i n hydraulic head and is called the specific storage Ss in terms of head [L'^].

on the left hand side represents the ease with which the medium The term

1 • transmhs water and is called the hydraulic conductivity K [LT' ] . Using these definitions equation 11 can be written as:

da K dh_ " da + • dp dh + dy[^ ' dy dh dt Eq. 12

Equation 12 is the constant density flow equation which is used in the Modflow2000 code (Harbaugh, A.W. et a l , 2000).

3.1.2 Non-Constant Density Flow Equation (Seawat2000)

To specify equation 8 for a non-constant density case a suhable expression for pressure needs to be used. With the expression of pressure of equation 9 in the non-constant density case there is next to the variables elevation and hydraulic head, the variable density. Therefore, substitiition of equation 9 for pressure into equation 8 in a variable density system will lead to cumbersome expressions. This can be avoided by defining pressure in terms of freshwater-density and freshwater-head (Guo, W. et al., 2002): P = Pfg(hf -z) where: Eq. 13 Pf h : freshwater density [ML" ] ƒ : equivalent freshwaterhead [L]

In this expression the only variables are the equivalent freshwater-head and elevation. The pressure derivative can be expressed as follows:

dP_ di dh, Pfg-di PfS dz_ di Eq. 14

Where a,p,y,t can be substituted for /.

Since the elevation does not change with time, the last term of equation 14 cancels for the time derivative.

_5_ da PK dh + • dy Pfg ƒ dz dz w da dh, + </• dz dz PfS-^-PfS-^,^PS V dy dy dy J ) ) dp Pk/ dh, dz dz + pq s = d_hf_^Qd^dC dt dC dt Eq. 15

Using the definition of the freshwater hydraulic conductivity iC^ = - ^ ^ ^ and

rearranging, the following equation is obtained:

d_ da i + — dy f = P V p K f . dhj da - + dhj dy 1 P-Pf . Pf ( dz da + -- + p-pf Pf dz dy )) dp + p q . p K f p ^ dhj Tp - + p-pf . pf dz dp d[hf-z} pd{hj--z) dc dt dt Eq. 16

Assuming that the difference between the compressibility coefficients, which is given by— of fresh and saline water can be neglected, the

term 1 dp P^hf-zY de ^ 1 ^ 1 ^QP ' d[hf-z) pd[h^-z) V

can be replaced by the coefficient of specific storage

in terms of freshwaterhead •$'ƒ [L' ] .

The coefficient of specific storage in terms of freshwater head represents the volume of water released from storage in a unit volume of aquifer per unit decline in freshwater head (Guo, W. et a l , 2002).

Furthermore, with the assumption that the viscosity of dense-water is essentially the same as that of freshwater because changes in viscosity are primarily the resuh of

temperature differences, the term ^ can be considered to be 1.

da 'd^ •+ +-dy pK fr dhj dy - + p - p f . Pf . ' p - p f Pf dz_ da + • df3 pK fP\ dz dy + p q s = ps dhj 'dp dh - + p - P f . Pf dz_ dp ƒ JJ dt dC dt Eq. 17

This is the flow equation for a non-constant density system in terms of freshwater-head which is used in the Seawat2000 code (Guo, W. et a l , 2002). For the special case where the aquifer has a horizontal bedding the horizontal and vertical x-,y-, and z-coordinate system can be used and equation 17 can be written as:

dx _V pK _{V ^ J} + - _pKfy ^dh, V ^ J d pK _fi dhj P - P f Pf _JJ + p q s dt dC dt Eq. 18

3.2 Transport Equation

The following solute-transport equation describes the fate and transport of component yt in a three-dimensional, transient groundwater flow system and is used in both the density-independent (MT3DMS) and density-dependent (Seawat2000) case:

d[9C')^ d dt da, where: -(g,C'')+g,C!+XR„ da^ Eq. 19 0 : porosity [-]

C* : dissolved concentration of component k [ML" ]

a. . : distance along the respective coordinate axis [L]

''^ 2 - 1

D.J : hydrodynamic dispersion coefficient tensor [L T" ] q_ : i component of specific discharge vector [LT"*]

: volumetric flow rate per unit volume of aquifer representing sources and sinks [T"']

C ƒ : concentration of the source or sink flux for component k [ML"^]

Y R„ • chemical reaction term [ML'^T"*]

term represents fluid sources or sinks within the REV and the last term is the chemical reaction term (Zheng, C. et a l , 1998).

Taking only sorption into account, the reaction term of equation 19 is given by:

ds

Eq. 20 where:

: bulk density of the soil = (l - 6»)/?,,,^ [MU^]

s : mass of adsorbed solute/ mass of bulk soil [-]

Assuming linear equilibrium sorption s can be expressed as follows (Fetter, C.W., 1999):

s = K,C'

Eq.21 where:

: distribution coefficient [L^M"']

Substitution of equation 20 and 21 into equation 19 and with the assumption of a constant porosity system the following expression for transport is obtained (Zheng, C. etal., 1998): ^8C^ 8 8t 8a, 8C or R9 dC' ^ 8 dt 8a, 0D, 'A 8a J J 8a. Eq. 22 Where R is the retardation factor [-] which is defined as:

Eq. 23

3.3 Coupling of Flow and Transport

For the density-independent case (Modflow2000 + MT3DMS), the flow and transport equation are uncoupled. First the flow equation is solved for the whole simulation period and, subsequently, this flow-fleld is used to solve for transport (Guo, W. et al., 2002).

For the density-dependent case (Seawat2000) the flow equation and transport equation are solved sequentially (two-step scheme). The resulting actual hydraulic heads from the flow equation are used to calculate the specific discharge, which is further used in the transport equation to calculate the transport of solutes (Guo, W. et a l , 2002). Next, the resulting solute concentration distribution of one of the component is used to update the density distribution. For this purpose a linear relation between density of saltwater and concentration developed by Baxter and Wallace (1916) is used:

p = Pj-hEC

Eq.24 where:

E : . Approximately 0.7143 (Guo, W. et a l , 2002) for the density-dependent dC

case (Seawat2000), 0 for density-independent case (Modflow2000)

C : salt concentration [ML'^]

4 Dense Flow in Horizontal Aquifer

Barenblatt, G.I. found in the early '50s [Barenblatt, G.I. et a l , 1984] a similarity solution to the 1-D flow equation for water flow in an isotropic initially dry 2=D horizontal aquifer using a sharp interface model. This similarity solution describes the evolution of groundwater level from an instant source of freshwater. Based on the similarity of the flow equation for wetting front and front of dense-water, the possibility to use a modified similarity solution to describe the evolution of dense/fresh water interface from an instant source of dense-water in a (non)-flowing iso/anisotropic imtially freshwater aquifer was tested in this study. Such a similarity solution can be of use to estimate model extend and parameters for density-dependent flow problems.

4.1 Problem Description

The objective of this part of the study was to check the possibility to modify the similarity solution, describing the evolution of groundwater level from an instant source of freshwater in an initially dry aquifer, to describe the evolution of the dense/fresh water interface from an in instant source of dense-water into a freshwater aquifer. This was done for a non-flowing/flowing and iso/anisotropic environment.

/ / / / / / / / / / Xs 4 w _H q ho H / h(x,t) / / / / / / /

^

I 1 Instant source of water / / Aquitard

Figure 5 Illustration of Barenblatt's model.

dt Adx dx

Eq. 25 where:

h : head of groundwater level above the aquhard [L] 0 : effective porosity [-]

Kq : hydraulic conductivity of saturated part of the aquifer [LT" ]

With the following initial conditions: h{xfi) = \ for - K ^ s < ^ < > 2 ^ . and

/2(x,0)= 0 outside this area, Barenblatt, G.I. found a similarity solution for equation 25 describing the evolution of groundwaterlevel h (Barenblatt, G.I., et al. 1984):

Kx,t) = with: a'0't 1/3 18 2 / 3 12 1 -{2xf Eq. 26 Xo(0

=

: initial head of groundwaterlevel above aquitard [L] : distance [L]

: length of base of the water plume [L] : initial width of the water plume [L]

Flow of Dense-water in a Confined Freshwater Aquifer (Zhang. H. 1988)

where:

k : intrinsic permeability [L ] /J, : dynamic viscosity [ML'^T'*] q : specific discharge [LT"*]

1 2

P : fluid pore pressure [ML" T" ] y : specific weight [ML"'^T"'^]

: umt vector in the i direction [L]

Taking the curl of equation 27 results into:

p dy dy —curlq e„

K dx ' dy

Eq. 29 q can be expressed as:

diff dy/^ q = curli// =

dz dx

Eq. 30 where:

I// : stream function

Substituting equation 30 into equation 29 leads to:

K dy // dx

Eq.31

Suppose that at each time the interface between dense- and freshwater (Fig. 5) with respect to z=0 is given by

z = u(x,t)

In terms of M the specific weight becomes:

yix,z,t)^ {p,g - Pfg)H{u{x,t)- z)+ P f g where: Eq. 32 Eq.33 Pf Ps g z H freshwater density [ML"^] dense-water density [ML"^] gravhational constant [LT" ] elevation of the interface [L]

Inserting equation 33 into equation 31 results into the following jump condition at the interface: dif/ dy/ dn = 0 = sina \M J Eq. 34 Eq. 35 where:

s : tangential local orthogonal coordinate [L] n : normal local orthogonal coordinate [L]

a : angle of the tangent at the interface with the horizontal [°]

For the jump in the x component of the specific discharge at the interface the following expression can be obtained from equation 34 and 35:

= qj^ ( x , u , t ) - q,^ {x,u,t)= - - Pfg)

Using the expression:

K_ K M Pfg du/dx l + {du/dxf Eq. 36 Eq.37

where K is the hydraulic conductivity [LT*], resuhs in the following equation:

du/dx

Eq. 38 where:

Pf

The total discharge through the aquifer can be expressed by:

Q = {x, z,t)dz + jq J, (x, z, t)dz

Eq. 39 where:

2 1

q^^ : specific discharge of dense-water in x direction [LT'*] qj.^ : specific discharge of freshwater in x direction [ L T ]

Also the following fresh and sah water continuity equations hold:

dx dt and dx dt Eq. 40 Eq.41 where:

: total discharge of freshwater in x direction [L^T*]

Q^^ : total discharge of dense-water i n x direction [L^T"'] 0 : total porosity [-]

To simplify this model the Dupuh approximation is used. With the Dupuh approximation, one assumes that the horizontal component of the specific discharge (qx) is constant over the height in each fluid and jumps across the interface (Fig.5). With the use ofthe Dupuh approximation equation 39 can be written as:

Q = uq^, {x, z, t)+{h- u)q (x, z, t)

Eq. 42

Equations 38 and 42 can be solved together for g „ and q^^ :

( ^ Kpih-u) du/dx , Q

Eq. 43

/ X Kpu du/dx Q

Eq. 44

A differential equation for u can be obtained by combining equation 40 and 44:

^^ = ± l K p { h - u ) u ^ ^ + {h-u)Q dt dx[ ' \+{du/dxy

0-dt d x y ' • l + (ÖM/öxj' J

x = 2 t = 6 u to obtain: ^du 1 5 [ ^ / , Y du/dx Q— = {Kp\\-u)u—, , dt 4 5x1 1 + (5M/5X)' With two more assumptions:

dx

2) M « 1 which implies u«h

equation 46 can be simplified to :

dt A dx ^ du^ V ö x y

Eq. 46

Eq. 47

Equation 47 has a similar form as Equation 25. When K^=Kp QXiéh = u the two equations are equal.

Based on this similarity the possibility to modify the similarity solution (Eq. 26) of Equation 25 to obtain a solution of Equation 47 was tested. For this purpose in Equation 26 h was replaced with u and with Kp to obtain the following modified similarity solution: u{x,t) = with: Xo(0^ 1/3 18 2 / 3 f / • o „ \ 2 A 12 1 -(2x)^ X, 0 J Eq. 48 2,\ i m y t 1/3 e 19 Mo = eu^x^ where:

Uq : initial head of dense-fresh water interface above aquhard [L]

K : hydraulic conductivity [LT"*]

To make use ofthe similarity solution in a flowing aquifer the following modification was done: Mo" ' 1/3 18^"

r

1

_ k dH dx where: : hydraulic gradient [-] dx for X < 0 for X > 0 Eq. 49 Eq. 50

In the anisotropic cases the horizontal hydraulic conductivity is much larger than the

vertical hydraulic conductivity : —f-«1

To take this into account an effective hydraulic conductivity needs to be used.

The 2-D flow equation is given by:

" dxdx " dz dz dt

Eq.51 Dividing by J^, leads to:

d du_ d du _ 0 du ^'dx dz dz ~ dt

K

Because — « 1 the second term can be neglected resuhing in:

8 du _ 0 du dxdx~ dt

Eq. 53 Therefore in the anisotropic case the value of the horizontal hydraulic conductivity

K.^ can be used as the effective hydraulic conductivity K^jj-.

To check the validity of the modified similarity solutions to describe dense flow in an initially freshwater 2-D horizontal aquifer, the analytical resuhs were compared with the results of numerical modeling. Flow of the dense plume was calculated and modeled for a flowing isotropic system, a flowing isotropic system and two non-flowing anisotropic systems respectively. For the numerical modeling of dense-flow the Seawat2000 code was used.

100

4.2 Modeling Description

The model was 1000 meters long, 100 meters thick and 300 meters wide and horizontal (Fig. 6). The fmhe difference grid-cells are 10 m horizontal, 3 m vertical and 100 m wide. No flow boundaries surround the model at the top and bottom. On the left and right side constant head boundaries exist. For the non-flowing simulations the difference in head between the left and right sides of the model is 0 m. For the flowing simulations the difference in head between the left and right sides is 0.5 m. The width of the imtial plume is 100 m and exists over the whole thickness of the model. Other input parameters and values can be seen in Table 5.

100

1000

Distances are in meters _ Initial plume

Constant head boundary

No-flow boundary

Figure 6 Model used for numerical solution of modified Barenblatt problem.

Model input parameter Variable (if applicable)

Value Horizontal hydraulic conductivity (isotropic and anisotropic) HK 365 m/yr Vertical hydrauhc conductivity (anisotropic) V K ^ 36.5 m/yr Vertical hydraulic conductivity (anisotropic) 3.65 m/yr

Effective porosity He 0.1

Specific yield 0.1 _{; 1^}

Bulk density Ph 2000 kg/m'

Initial concentration plume component 25g/l Table 5 Model mput parameters for modified Barenblatt problem.

The total simulation period of 400 years was divided into 5 stress periods. Each stress period was further divided into time-steps (Table 6). Seawat2000 uses stability constraints and accuracy reqmrements (Guo, W. et a l , 2002) to further divide the time-steps into transport-steps.

Stress Period Period length #of time-steps

I JO 10_

2 90 180 3 iOO 200

4 1ÖÖ 2ÖÖ

5 100 200

Table 6 Partition of simulation period into sttess periods and tune-steps.

Note that the theory behind the derivation of the 1-D flow equation (Eq. 47) to describe dense/freshwater flow in a 2-D aquifer reqmres the height of the initial plume to be much smaller than the height of the model (Fig. 5). The model used m this research (Fig. 6) does not f u l f i l l this requirement. To check whether this model (Fig 5) would still give the right resuhs some additional modeling was carried out with a larger model similar to the model in Figure 5. The resuhs were basically the same as can be seen in Figure 7.

Model Fi

I Model Fi

Distance (m)

4.3 Results

Figure 8 to 10 show the graphs ofthe migration of the waste plume in a respectively, non-flowing isotropic, flowing isotropic and non-flowing anisotropic model. The graphs show a distance versus height plot in meters, where height is the height of C0.5 (half the maximum concentration). The dashed line represents the resuhs obtained by numerical modeling and the solid line the analytical results calculated using the modified similarity solution. The anisotropy ratio of the plot in Figure 10 is VK/HK=3.65/365=0.01.

Figure 11 shows the calculated value of IQff in the modified similarity solution^usmg the numerical resuhs for the height of C0.5 at x=0 for an anisotropic (VK/HK=0.01) non-flowing horizontal aquifer (Fig.6).

(D X —m-r X —— f / I / ? 0 \ \ _\\ \ s (/ 1 1/ \ ^ _\\ \ \ _\v — </ ''fi // // / / \k \\ _\\ \^ (7 1/ it V\ _\\ u 1 \ / / '/ 1/ ' / 5 \\ \\ _\\ ll 1; li ' / 0 \\ _\\ h ))0 .600 -400 -200 * -1 200 400 600 81 t = 10 yrs t = lOOyrs Numerical Analytical Distance (m)

Figure 8 Numerical and analytical heights of C0.5 versus distance for an isottopic horizontal non-flowing model (Fig. 6).

- - - - — 4 0 ^ \ —™ — —— f/ }/ 1/ \ \ ^ \ \ — — ^ — —3&¬ d- 9ft-\ 9ft-\ \ \ \ —i 2(h ™ _ „ +S-— r-V 1 4Ö-/ 4Ö-/ h c

V

J t=10yrs t = 100 yrs t = 400 yrs Numerical Analytical -800 -600 -400 -200 0 200 Distance (m) 400 600

Figure 10 Numerical and analytical heights of C0.5 versus distance for an anisotropic (VK/HK-0.01) horizontal non-flowing model (Fig. 6).

100 160

200 250 300 360 400 460

Time (yrs)

Figure 11 Calculated K,ff from numerical results of an anisotropic (VK/HK=0.01) horizontal non-flowmg model (Fig. 6).

4.4 Discussion

To find the right relation between the effective hydraulic conductivity Keff and the horizontal and vertical hydraulic conductivity for this model empirically, the Keff for the modified similarity solution was calculated using the numerical resuhs ofthe C0.5 height at x=0. The resuhs can be seen in Figure 11. From here h can be seen that the Keff is time dependent. In early times flow is mostly into the vertical direction resulting in a low Keff. After awhile the components reach the bottom of the aquifer and horizontal flow dominates resulting in a large Keff- This explains the better agreement for later times. This time-dependency of the effective hydraulic conductivity makes the modified similarity solution not very suhable for anisotropic cases since the time-dependent relationship is not a simple linear one.

4.5 Conclusions

The modified similarity solution can be used to describe dense flow in an imtially freshwater aquifer for a non-flowing and flowing isotropic case.

5 Effect of Density on the Migration of L R W in a Sloping Aquifer

Before modeling the migration of the two-component LRW for the complicated Krasnoyarsk-26 she, a simple test case was carried out to form a clear picture o f t h e different processes taking place. Two simple 2-D rectangular test-models were used: a horizontal model 1 and a dipping model 2. For this test case, the influence of density on the migration of both components of the LRW for an isotropic and anisotropic environment was analyzed by performing a comparative study ofthe different driving forces behind the migration. Furthermore, the importance of the right choice of grid and scheme was looked at.

5.1 Problem Description and Methods

5.1.1 Problem Description

The objective is a comparative study of the different driving-forces behmd the migration of LRW to analyze the influence of density on the migration of both components of the LRW. Based on the resuhs of the [IGEM, 1997] (Section 2.2) temperature effects were not taken into account and, therefore, this study was limited to the coupling and decoupling of regional-, density-dependent- and gravhational convection. An isotropic and two anisotropic modeling geometries were considered to analyse the influence of aquifer vertical anisotropy on dense plume migration.

Two models of the same extend were used. Model 1 consisted of a 2-D rectangular horizontal aquifer (Fig. 6). Model 2 consisted of a 2-D rectangular dipping aquifer with flow in the upslope direction (Fig. 12).

For both models flow and transport were modeled for a density-dependent and density-independent case. Table 7 shows which driving force exerts influence on the migration ofthe plume in the different simulations:

Regional Density-dependent Gravitational Model 1 X X Density-dependent ^ Model 1 X Density-independent . . Model 2 X X x Density-dependent . Model 2 X Density-independent . Table 7 Active driving forces during the different simulations.

The LRW plume was represented by two components: A dense non-sorptive component 1 representing most ofthe waste and a sorptive component 2 representing the radionuclides (Section 2.1). Linear equilibrium sorption was assumed.

To keep the model as shnple as possible only the advection transport mechanism and sorption as chemical reaction were taken into account.

Second, the density-dependent and density-independent migration of waste for model 1 and 2 were simulated for a flowing and non-flowing system for both the isotropic and anisotropic environments using the grid and schemes that showed to give the best results.

The resuhs ofthe different simulations were qualitatively compared by looking at the visual resuhs and quanthatively by looking at the center of mass (section 5.1.5) ofthe plume both at different time instances during the simulation.

To decouple the regional convection from gravitational and density-dependent convection, the independent runs (Modflow2000+MT3DMS) and the density-dependent runs (Seawat2000) for each model were compared. Differentiation between the contribution of density-dependent convection and gravhational convection was done by comparing the density-dependent modeling resuhs (Seawat2000) of the horizontal model 1 with the density-dependent resuhs of the dipping model 2.

5.1.2 Model and Modeling Description

The description of model 1 can be found in section 4.2 of this report.

Model 2 was 1000 meters long^ 100 meters thick and 300 meters wide and made an angle with the horizontal of 5 degrees. Two different model grids were used: Grid 1 with horizontal gridlayers and grid 2 with dipping gridlayers parallel to bedding ofthe aquifer (Fig. 12). Both grids have the same resolution; each model gridcell is 10 m horizontal^, 3 m vertical and 100 m wide.

No flow boundaries surround the model at the top and bottom. On the left and right side constant head boundaries exist. The difference in head between the two sides is 0 m for the non-flowing simulations and 0.5 m for the flowing simulations. Flow is i n the up-dip direction. The width ofthe initial plume is 100 m and exists over the whole thickness of the model.

The initial concentration of the plume is set to 25 g/1 for both component m both model 1 and 2. The actual density of component 1 of LRW is around 250 g/1. The use of this high value caused the non-sorbed component to migrate beyond the model boundaries within the first few years. Therefore, a lower value of 25 g/1 was used for the simulation. The density of component 2 is much smaller than this but setting the concentrations of component 1 and 2 equal makes h easier to compare. Other input parameters and values for both models can be seen in Table 8.

Model input parameter Variable (if applicable)

Value Horizontal hydraulic conductivity (isotropic and anisotropic) HK 365 m/yr Vertical hydraulic conductivity (anisotropic) V K ^ 36.5 m/yr Vertical hydraulic conductivity (anisotropic) VK,„ 3.65 m/yr

Effective porosity De 0.1

Specific yield Ss 0.1

Bulk density Pb 2000 kg/m'

Initial concentration non-sorbed component 1 Cjii 25g/l Initial concentration sorbed component 2 25g/l Distribution coefficient sorbed component 2 Kd 1 xlO-Wd Table 8 Input parameters for model 1 and 2.

^ For the later runs the length ofthe model was extended to 2000m.

10

A •

10

3

t

Distances are in meters | | Initial plume

Constant head boundary

_ _ No-flow boundary

The Layer Property Flow (LPF) package was used for formulating inter-cell hydrauhc conductance terms. To solve the flow equation, the PCG2 solver was used. For the advection part ofthe transport equation two different solutions schemes were used: 1) The Upstream Finite Difference method with an upstream weighting based on a courant number of 0.75. 2) The S'" order TVD scheme also based on a courant number of 0.75. The transport equation was solved by usmg the GCG solver. To couple the flow and transport equation in SEAWAT2000 both, explich and implich schemes were used (appendix A). . j . ^ The total simulation period of 400 years is divided into 5 stress periods. Each stress period is fiirther divided into time-steps (Table 9). Seawat2000 uses stability constraints and accuracy requirements (Guo, W. et a l , 2002) to fiirther divide the time-steps into transport-steps.

Stress Period Period length #of time-steps 1 10 10

90 270 100 300 100 300 5 100 300

Table 9 Partition of simulation period into stress periods and thne-steps.

5.1.3 Importance of Grid

Correct representation of the hydraulic conductivity is important for the accurate solution of hydraulic heads, flow and transport (Anderman, E.R. et al., 2002). the dkection o f t h e vertical and horizontal hydraulic conductivity should coincide with the principle axes of permeability. In both Modflow2000 and Seawat2000, the codes used in this study to solve for flow and transport, the direction of the vertical and horizontal conductivity are ahgned with the vertical and horizontal axes of the grid. Therefore the axes ofthe grid should coincide with the main axes of permeability, ^or an isotropic case, the case where permeability is independent of direction the correct representation of the hydraulic conductivity is not dependent on grid. For the anisotropic case, however, permeability is dependent on direction and a correct choice of grid is required to get good resuhs. Both grids for model 2 used m this study do no fiilfiU this requirement for the anisotropic case. For grid 1 with the horizontal gridlayers, both the representation of the horizontal hydraulic conductivity and the vertical hydraulic conductivity are incorrect. For grid 2 with the dippmg gridlayers the horizontal hydraulic conductivity is represented right but the vertical hydraulic conductivity contains the same error as grid 1.

used for this is bounded by constant head and no-flow boundaries. The boundaries force the flow to go into the right direction. The eiTor in grid 2 should be less than grid 1 since the main direction of flow has the right representation of hydraulic conductivity.

5.1.4 Schemes

For each grid, density-dependent and density-independent flow and transport were modeled for an isotropic and two anisotropic cases for the following schemes:

Density-independent:

1) Modflow2000 + MT3DMS run using TVD solver for advection.

2) Modflow2000 + MT3DMS run using Finite Difference solver for advection.

Density-dependent:

3) Seawat2000 explich mn using TVD solver for advection.

4) Seawat2000 explich run using Finite Difference solver for advection. 5) Seawat2000 implich run using the Finite Difference solver for advection.

The Modflow2000 + MT3DMS run, models the migration of the plume whhout taking density differences into account and, therefore, only regional convection is modeled. The flow and transport equation are uncoupled.

The Seawat2000 explich and implich schemes both take density difference into account by coupling the flow and transport equation. In this way all three driving forces are modeled. The difference between the explich and implich approach is the way of coupling. The explich scheme uses fluid densities from the previous transport step and, thus, is actually one time-step behind. The implich scheme solves the flow and transport equation repeatedly for each transport step until consecutive differences in the calculated fluid densities are less than a user-specified value (Guo, W. et a l , 2002; Langevin, C D . et al., 2003).

Two different solvers for the advection part of the transport equation are used: the Finite Difference solver and the TVD solver. The Finite Difference solver is required for the Seawat implich run.

More information about the different codes and solvers can be found in Appendix A.

5.1.5 Center of Mass

The different schemes were quanthatively compared by looking at the position ofthe center of mass at different instances during the simulation.

'{xxc(x,z,t)dxdz \\zxc(x,z,t)dxdz c(x,z,t)dxdz ^jc(x,z,t)dxdz '{(x-X"""'fxc{x,z,t)dxdz , {{(z-Z"""'yxc(x,z,t)dxdz a / (0 = ^ ,(t' (O = f7-\\c(x,z,t)dxdz j^c{x,z,t)dxdz Eq. 54 where:

jmms ^ 2>»«'' : X and z component of center of mass respectively c{x,z,t) : concentration

(T^, : Standard deviation in the x and y direction respectively

5.2 Results and Discussion

The visual resuhs of schemes 1 and 2 and schemes 3, 4 and 5 are not distinguishable from each other. Therefore, the determination of the best scheme was based on the poshion ofthe center of mass. Figure 13 shows the center of mass of both components for the two different grids and for the iso- and anisotropic case at different time instances during the simulation using the different schemes. The strange wiggle for the center of mass of component 1 for the later times in Figure 13a,b,d is caused by boundary effects i.e. the model extend used for the first simulations was not long enough for a simulation period of 400 years.

From Figure 13 h can be seen that the choice of solver for advection in the transport equation is crhical. Using the Finite Difference solver for advection, the implich and explich schemes give almost exactly the same resuhs. When the TVD solver is used the explich and implich schemes start to deviate. The direction of flow stays about the same but the magnitude of flow for the TVD is a lot higher (except for component 2 for grid 1 in the anisotropic case). This trend can also be seen when comparing component 1 for scheme 1 and 2, the direction for the Finhe Difference and T V D solver is the same but the magnitude using the TVD solver is larger (except for grid 1 in the anisotropic case).

For the isotropic case the use of the TVD solver causes oscillations in the concentrations of component 1. The running time for schemes 1 and 2 and schemes 3, 4 and 5 were comparable.

From this h can be concluded that scheme 2, 3 and 5 are good schemes to use to numerically solve the problem in this section. In addhion, the documentation of the Seawat2000 code (Langevin, C D . et a l , 2003) states that for simulations with rapidly changing solute concentrations, the implich coupling option may provide a more accurate solution than the explich option. Therefore, the resuhs of schemes 2 and 5 ai-e used in the remainder of this report.

Legend Figure 13:

I Implicit (Finite) Component 1 I Explicit (Finite) Component 1 J Explicit (TVD) Component 1

Modflow (Finite) Component 1

I Modflow (TVD) Component 1

t = l

Implicit (finite) Component 2 Explicit (Finite) Component 2 Explicit (TVD) Component 2

I Modflow (Finite) Component 2 I Modflow (TVD) Component 2

Each dot represents the center of mass at a certain time. From the first dot to the last, the times belonging to each dot are as follows: 1-10-20-30-40-50-60-70-80¬ 90-100-120-140-..-..-360-380-400 years.

Distance (m)

0 --20 --30 --40 - -60--60 -c -70 -c -80 --90 -+ J -90 -Q. -90 -0) -100 -Q -100 --110 -120 -130 -140 -160 -160 • -170 -180 420 1020

iB2

Figure 13b Distance versus depth plot ofthe center of mass at different time mstances during the simulation using grid 2 for an isottopic envkonment.

-10 -20 -30 ^0 -60 ^ -70 ^ -80 £ -90 Q. 0 -100 ^ -110 -120 -130 -140 -160 -160 -170 -180 720 Distance (m)

Distance (m)

Figure 13d Distance versus depth plot of the center of mass at different thne instances during the simulation using grid 2 for an anisotropic envhonment.

5.2.2 Importance of Grid

The resuhs for the different grids (Fig. 12) in an isotropic and anisotropic (VKTHK^O.Ol) enviromnent for a density-dependent and independent case can be seen in Figure 14. Each dot in Figure 14 represents the center of mass at a certain time. From the first dot to the last, the times belonging to each dot are as follows:

1-10-20-30-40-50-60-70-80-90-100-120-140-..-..-360-380-400 years.

1

Grid 1 Species 1 (scheme 5) Grid 1 Species 2 (scheme 5) Grid 2 Species 1 (scheme 5) ; Grid 2 Species 2 (scheme 5) -«-Grid 1 Species 1 (scheme 2) Grid 1 Species 1 (scheme 2) Grid 2 Species 1 (scheme 2) — Grid 2 Species 2 (scheme 2) Grid 1 Species 1 (scheme 5) Grid 1 Species 2 (scheme 5) Grid 2 Species 1 (scheme 5) ; Grid 2 Species 2 (scheme 5) -«-Grid 1 Species 1 (scheme 2) Grid 1 Species 1 (scheme 2) Grid 2 Species 1 (scheme 2) — Grid 2 Species 2 (scheme 2) Grid 1 Species 1 (scheme 5) Grid 1 Species 2 (scheme 5) Grid 2 Species 1 (scheme 5) ; Grid 2 Species 2 (scheme 5) -«-Grid 1 Species 1 (scheme 2) Grid 1 Species 1 (scheme 2) Grid 2 Species 1 (scheme 2) — Grid 2 Species 2 (scheme 2)

Distance (m)

Grid 1 Species 1 (sclneme 5) Grid 1 Species 2 (scheme 5) Grid 2 Species 1 (scheme 5) Grid 2 Species 2 (scheme 5)

- 5 « - Grid 1 Species 1 (scheme 2) Grid 1 Species 2 (scheme 2) Grid 2 Species 1 (scheme 2) Grid 2 Species 2 (scheme 2)

Distance (m)

Figure 14b Distance versus depth plot of the center of mass at different time mstances during the simulation for the density-dependent (scheme 5) and density-mdependent (scheme 2) run for grid 1 and 2 in an anisotropic environment.

It can be seen that for the isotropic case the resuhs for the two grids without density taken into account (scheme 2) are similar. This is as expected since the main axes of permeability are independent of direction in an isotropic environment. When density is taken into account (scheme 5) the results of the two grids start to deviate: The direction of flow is the same for both grids, the flow magnitude differs. The influence of density is larger in the dipping gridlayer grid 2: the non-sorbed component 1 is carried fiirther down-slope and the sorbed component 2 is carried less far upslope. This can be explained by the governing equation for flow used in Seawat2000 (equation 17). The density term in the flow equation depends on the gradient of the elevation in the main directions of the grid. The gridlayers for grid 1 are horizontal, therefore the elevation gradient is zero along this axes and the density-dependent term disappears. For the dipping grid the elevation gradient is non-zero in the gridlayer direction resulting in a density term along this direction. Thus, for grid 1 there only exists a density term for the vertical direction and for grid 2 there is a density term for the gridlayer and for the vertical direction causing a higher influence of density.

-70 -110 ^ -130 Q. 0 Q -160 -170 1200 1400 Isotropic Species 1 1 Isotropic Species 2 Anisotropic Species 1 Anisotropic Species 2 Distance (m)

Figure 15 Distance versus depth plot ofthe center of mass at different tune mstances during the simulation for the density- independent (scheme 2) run for grid 2 for an isotropic and anisotropic environment.

Each small yellow/blue arrow represents one time-step. For the isotropic case in both the vertical and horizontal directions the same distance is travelled during one time-step. In the anisotropic case the distance travelled in the vertical direction is smaller, because of hs lower vertical hydraulic conductivity compared to the horizontal hydraulic conductivity. After the specified time period, the blue line for the isotropic case has travelled a larger distance with a bigger slope than the yellow line for the anisotropic case.

Grid 2 iso- & anisotropic

Grid 1 isotropic

Grid 1 anisotropic