Acta Polytechnica

This work deals with a groundwater flow and solute transport model in the near-surface (predominantly unsaturated) zone. The model is implemented so that it allows simulations of contamination transport from a source located in a geological environment of a rock massif. The groundwater flow model is based on Richards’ equation. Evaporation is computed using the Hamon model. The transport model is able to simulate advection, diffusion, sorption and radioactive decay. Besides the basic model concept, the article also discusses potential cases that could lead to non-physical solutions. On three selected examples, which include, for example, rapidly changing boundary conditions, the article shows the solvability of such cases with the proposed model without unwanted effects, such as negative concentrations or oscillations of solution, that do not correspond to inputs.


Introduction
In the Czech Republic, it is planned to dispose of the spent nuclear fuel within the deep repository in hard crystalline rock.It is a similar concept to those adopted in Sweden (www.skb.se) and Finland (www.possiva.fi)[1][2][3] with one of the main differences being the near zero terrain gradient over their coastal deep repository sites.In the Czech Republic, the repository candidate sites have an average altitude of about 500 m with locally steep gradients, which means that it is necessary to focus closely on both saturated and unsaturated flow and solute transport.
The biosphere model is a part of the safety analyses for ensuring the safety of a planned deep repository of spent nuclear fuel and highly active wastes.For a purpose of computation of a possible radionuclide impact on biosphere it is necessary to know their distribution (concentration or activity) in its vicinity, i.e. in a groundwater (near its level) where there is a possibility of direct extraction via pumping and also in the unsaturated zone above the groundwater level.A prediction of radionuclide concentration distribution in groundwater is done by transport simulations along with groundwater flow simulations (both mainly in saturated zone).Such a model usually doesn't provide us with enough detailed concentrations in a zone near the surface where the solute transport has an entirely different character (due to the low saturation).This problem is not new.It is being studied both in the context of deep repositories and potential environmental contamination.Existing SW tools (Pandora [4], ResRad [5], Hydrus [6] used in this field do allow for transport simulation in the unsaturated zone, but they show numerical instabilities in the case of dynamically changing flow boundary condition.This article uses case studies to show how to deal with such instabilities. With a biosphere module in mind, it is necessary to evaluate the near-surface radionuclide distribution using a more detailed model of the unsaturated zone (based on a known concentration distribution in a saturated zone and expected precipitation amounts).For such simulations, there is a variety of existing tools.An overview of their selection along with their characteristic features may be found in Steefel, Appelo and Arora [7].The alternative is to implement one of the well-known concepts for transport in unsaturated and partially saturated environments, which is the case in a work presented in this article.The aim was to create a software tool for biosphere simulations, which uses Flow123d [8] and [9], verified against other codes in [10] for transport simulations in the saturated zone (geosphere).Flow123d provides us with pressure distribution, velocity vectors and radionuclide concentrations in both mobile and immobile phases of a rock massif saturated environment.For the calculation of a radionuclide concentration time-space distribution in an unsaturated zone of defined model subdomain, we use the precipitation amount data.So far, in the phase of model implementation and testing, we used data from meteorological station Praha Klementinum, which are (for daily precipitation amounts) available for past ca 200 years.The implementation includes the Hamon evapotranspiration model [11].The unsaturated zone flow simulation is based on Richards' equation and transport simulation is based on the advection-diffusion equation.For the approximation of both flow and transport, the numerical scheme shown in [6] and [12] is used (as is the case of, for example, Hydrus [6]).
This article deals with some of the aspects of the unsaturated flow and transport model implemented as an extension of a 3D hydrogeological and transport model of a site with a deep repository as a source of radionuclides that have the potential to migrate from the repository to a biosphere.The unsaturated zone model outputs serve as an input for radionuclide transport computations in a biosphere.The much needed mutual interconnection of all the modules along with the necessary customization of unsaturated zone model input and output form is the main reason for our own implementation, instead of using existing SW such as Hydrus.This also means that there is no need for installation or execution of any third party SW tools.We managed to ensure 1) same data structure for transport simulations in both saturated and unsaturated zones, which allows us to avoid difficult interconnection of inputs and outputs of various simulation SWs and 2) a support for sensitivity analysis of the whole path from a source to a biosphere.The downside is that we had to deal with problems linked to numerics, such as solution oscillations or negative concentrations.Besides the model itself, the article also shows three case studies to present the model's usability, including the sensitivity analysis performed in accordance with [13].

Model concept
The unsaturated zone model was implemented with the aim to simulate the radionuclide transport from a geosphere to a biosphere.It assumes a known timespace distribution of radionuclide concentration in a part of geosphere that is close to groundwater level.This concentration distribution is acquired from the transport model of a site with a deep repository, which has concentration values in given times and discrete points of a simulation domain as an output.In the case of Flow123d simulation tool, the time-dependent elementwise constant values of concentration are available.It is possible to identify the near-surface zones of the model domain for which the more detailed unsaturated zone simulations will be performed with the aim of acquiring a higher resolution of the concentration distribution near a surface.
While in the saturated part of the geosphere, the flow has predominantly horizontal character, in the unsaturated zone, the flow of water and dissolved contaminants is mainly vertical.This is why the unsaturated zone model was limited to one dimension in the direction of z axis (with the direction of rising altitude).For each selected near-surface element, the separate unsaturated zone simulation is performed with the following basic inputs: • Time evolution of contamination concentration in water -acquired from an output of the 3D transport model; values are used as a transport boundary condition for the 1D model, • Vertical profile of soil composition in the unsaturated zone (alternatively of the van Genuchten parameters), • Depth where the saturation is equal to one, i.e. depth of the groundwater level, • Time evolution of precipitation amounts; used as a flow boundary condition, • Time evolution of evaporation values (implemented via Hamon model [11]), • Time evolution and location of sources in a model domain.
The definition of initial conditions requires the lower part of the model domain to be saturated and, in sake of numerical stability, requires the saturated part to be at least few meters thick.The initial flow conditions in the model domain as a whole are usually derived from the groundwater level (hydraulic/piezometric head) which could be determined by a direct measurement or an expert estimation.The flow boundary condition on the top of the model domain is given by the volume of the water entering the profile (could be time dependent) and on the bottom, by hydraulic head (could also be time dependent).Sources could be defined to cover the root withdrawal.
For the transport simulation, the top boundary condition is usually zero concentration (when we assume the source of contamination to be located deep in the geosphere).The bottom boundary condition is defined based on the concentration values acquired from a simulation of the saturated zone.This boundary condition interconnects the saturated and the unsaturated model.The transport model includes sorption, diffusion and radioactive decay.
The 1D model domain is discretized with a given step.The values of pressure head, water content and flux in each discretization node are the outputs of the flow simulation.The values of concentrations in each discretization node are the outputs of the transport simulation.

Model of flow
In the saturated zone, all pores are filled with water and the pressure head h ≥ 0. The groundwater flow could be both horizontal and vertical.And in the unsaturated zone, the pores could be filled with water or with air.The pressure head h < 0. The flow is predominately vertical.That is why we limited the model to one dimension.
The flow in the unsaturated zone is described by Richards' equation [14]: where h [m] is the pressure head, θ [1] is the water content, t [s] is the time, z [m] is the vertical axis, S [s −1 ] are the sources and K [m s −1 ] is the unsaturated hydraulic conductivity.
In equation ( 1), the principal variable is h.Variable θ(h) and parameter K(h) are linearly dependent on the current value of pressure head h.According to van Genuchten [15], this dependence could be described by following equations: where θ r [1] is the residual water content, θ s [1] is the saturated water content, α [m −1 ] is the inverse of the air-entry value, n [1] is the empirical shape-defining parameters and K s [m s −1 ] is the saturated hydraulic conductivity.
Van Genuchten parameters α and n depend on the soil composition: fraction of clay, silt and sand and bulk density of soil.There is no analytical formula linking the soil composition to van Genuchten parameters.Their values could be estimated, for example, by using the Rosetta SW [16] which is based on a quasi-empirical model.
For a particular simulation, it is necessary to also input (aside from soil parameters) the boundary and initial conditions.On each model boundary (there are only two since the model domain is 1D) it is possible to input either flux or pressure head.The boundary conditions could be changed in the course of simulation, both their type and their value.In each time step, it is necessary for the pressure head to be specified on at least one boundary.The initial conditions (in the form of pressure head values) are derived from the groundwater level in the model domain.It is possible to define sources S [s −1 ] which represent water withdraval through roots or a well.
The model domain is divided into a mesh of n − 1 line elements with n nodes.Richards' equation (1) has no analytical solution; in a model, it is solved using the Piccard numerical scheme [17].The solution is in a form of pressure head in each computational node and each simulation time step.

Model of transport
In the unsaturated zone, the contamination may generally exist dissolved in three phases -liquid, solid and gaseous.In our concept, we do not account for a contamination transfer to gaseous phase.The reasoning behind this is the long simulation time in comparison to rock residence time of isotopes in the gaseous form.
The concentration of dissolved contaminants/isotopes is very low, which is why the model uses only the linear sorption isotherm.The transport model also includes diffusion and radioactive decay.In each simulation time, we also assume an equilibrium between concentrations in solid and liquid phases given by the distribution coefficient.Based on these assumptions, the dependence of concentration on time and space could be described by the advection-diffusion equation ( 4) [6]. where is the radioactive decay of isotope k in liquid and solid phases, γ w,k [kg m −3 s −1 ] is the zero order reaction in a liquid phase, γ s,k [s −1 ] is the zero order reaction in a solid phase, r k [kg m −3 s −1 ] are sources.The relationship between concentrations in solid phase s k and in liquid phase c k is given by equation (5).
where k D,k [m 3 kg −1 ] is the distribution coefficient of linear sorption for isotope k.Boundary conditions on both ends of the model domain are given in a form of concentrations in the liquid phase.The model allows for the boundary conditions to change over time.
The initial conditions are given as concentrations in the liquid phase in individual computational nodes; the initial concentrations in the solid phase are then computed using equation (5).The differential equation ( 4) is numerically solved using the finite elements method.The computation uses the same discretisation as the model of flow.The solution is in the form of concentration of each contaminant in each node and at each simulation time.

Model instabilities
While simulating the flow and transport using the model described above, some instable states may arise.They needed to be dealt with during the implementation.Potential instabilities must be eliminated by an appropriate choice of parameter values.This section provides an overview of known potential unstable states along with a method that was used to deal with them in the model implementation.
The evaporation from the model surface may be defined by outflux boundary condition.This may cause a non-physical state where the withdrawn water doesn't exist in the model domain.This would lead to results limitedly approaching h → −∞ and θ → θ r , which would require the simulation time step to be ∆ t → 0. A possible solution is to define the evaporation as a source (water is withdrawn from the defined depth).This approach allows for the simulation of 1) water withdrawal through roots, 2) evaporation and 3) near-surface water circulation.Another possible solution is to define the model domain extent so that its bottom part remains saturated throughout the simulation.If the defined initial conditions imply non-zero flux in the model domain, it results in a rapid stabilization in the first few simulation time steps, which may lead to significant computed groundwater fluxes.This may lead to a non-physical solution and a necessity to shorten the simulation time step.To control the time step, the Courant number defined as Cr = |q•∆t| θ•∆z [1] may be used.The solution is stable when Cr ≤ 1.Based on the knowledge of max i Cr i (here, the lower index i stands for the element of discretisation) the time step may be adjusted accordingly.When changing the flow boundary conditions during the course of simulation, similar problems as in the previous case may arise.Those problems may be generally solved by lowering the simulation time step and making the changes in boundary conditions gradual, if possible.Because the flux may be discontinuous in time, it is necessary to estimate the Courant number and adjust the time step accordingly.
Unstable and non-physical solutions may be generally caused by long flow-simulation time step.The problem may be fixed by an estimation of the time step using the Courant number.
When simulating the transport, there may be oscillations in the solution and negative concentrations in cases where the advective flux dominates over the diffusion (along with suboptimal time and space discretisation).Maximum transport-simulation time step may be set based on Péclet number P e = q∆x θD [1], where q [m s −1 ] is the flux, ∆x is the element size [m], θ is the water content [1]  Oscillations of the solution and/or negative concentrations may also be caused by changing concentration in transport boundary conditions.Such problems may be solved by lowering the transport simulation time step.
An unstable solution may occur in the case of a too small model domain with respect to the pressure head value given as the flow boundary condition.If the bottom boundary condition is a pressure head, which is not h ≫ 0, then, if the precipitation amounts are high, the saturated/unsaturated zone interface moves up but the boundary condition remains the same.This results in a significant pressure head gradient in nodes close to the bottom boundary and consequently, unrealistically high fluxes.Such problems may be solved by extending the model domain so that it captures a greater part of the saturated zone.

Case studies
This part of the article shows three case studies.The first one is a synthetic task that includes a simulation of flow and transport with regularly changing boundary conditions.The second one examines the sensitivity of the output on the precipitation amount.And finally, the third one shows the expected radionuclide migration in the case of a realistic precipitation boundary condition (data from Praha Klementinum; available since 1804).This case also includes the evaporation computed using the Hamon model, which is based on temperature measurements, latitude, and date (last two quantities are used for the computation of daylight length).Bottom boundary conditions as well as initial conditions were, in all three cases, derived from the saturated transport simulation with contamination source located 500 m below surface.This regional saturated model is not discussed here; parameters derived from it are shown where relevant.Case studies are evaluated based on the concentration evolution (depth dependent).

Case 1 -periodic precipitation evolution
The first case simulates a flow and solute transport from a 10 m depth towards a surface.The model domain is divided into 100 computational elements (101 nodes).The z axis is in the direction of rising altitude (the surface is at z = 10 m).The simulation period is 5000 days with a simulation time step of 1 day.
The flux, representing precipitation (after subtraction of evaporation), is used as a flow boundary condition on the top boundary.The time evolution of  the boundary condition has a sinus shape with a 365day periodicity.Its values for individual simulation times were (for the expected influx of 365 mm y −1 ) calculated as okp p F lux = 365−438•sin 2πt 365 [mm y −1 ].A constant pressure head of 5 m is specified on the bottom boundary throughout the simulation.This corresponds to a groundwater level 5 m below surface, which is also how the initial conditions were defined.
The flow simulation result is the time-space distribution of the pressure head and saturation.For example, at 5 m depth (where the groundwater level was situated at the simulation start), the pressure head oscillates between 0.009 and 0.061 m (see Figure 1, left).At one meter below the surface (see Figure 1, right), we can see an initial rise, followed by periodic oscillations between −2.61 and−1.17m (difference of ca 1.44 m).
The initial concentration for both isotopes is 1 • 10 −9 kg m −3 below the groundwater level (depths of 5-10 m) and zero above it.The bottom transport boundary condition is a constant concentration of 1 • 10 −9 kg m −3 .The model also assumes a horizontal flow in its saturated part, which means that the concentration there is kept at a constant value of 1 • 10 −9 kg m −3 (via prescription of nonzero sources in computational nodes) throughout the simulation.
Figure 2 shows examples of results at two depths: 4 m and 1 m.It is clear that after 2 000 days, the concentration evolutions adopt the periodicity of the precipitation evolution.Near the surface, the concentration is significantly lowered by an infiltration of clean water.In a realistic environment, the measure of dilution would be predetermined by many factors, such as soil type, unsaturated zone thickness, precipitation and evaporation amounts, roots, and wells.

Isotope
Half  1.The regional saturated model transport simulation period was 500 000 years.As a consequence of a given radionuclide release scenario and transport parameters, the following isotopes had a high-enough concentration near the surface at the end of the simulation: 36 Cl, 41 Ca, 79 Se, 129 I. Concentrations of other isotopes were less than 1E-15 g m −3 .For the four isotopes mentioned above, the unsaturated transport simulation was performed with an aim to determine the sensitivity of the concentration distribution on the infiltration amount.
The model domain is 10 m thick and discretized with a 0.1 m step.The simulation period was 500 000 y with a 0.1 y time step.The soil is sandy loam with the following parameters: θ r = 0.065, θ s = 0.41, The saturated part of the model domain is 8 m thick (groundwater level 2 m below surface).This corresponds to bottom flow boundary condition of a constant 8 m pressure head.The simulations were carried out with various top flow boundary conditions that correspond to annual infiltrations of 0, 10, 20, 30, 40, 50, 60 and 80 mm y −1 .
The transport initial condition was zero concentration throughout the domain.The bottom transport boundary condition had a form of concentration time evolution based on regional saturated model results (see Figure 3).These concentrations are also kept in the saturated part of the model domain throughout the simulation.
Table 2 shows the calculated values of 36 Cl, 41 Ca, 126 Sn and 129 I concentrations at simulation times of 50 000, 100 000 and 200 000 y and in depths of 0.3, 0.5 and 1.0 m as well as in the saturated zone.The results shown are for the infiltration amounts of 0, 20 and 50 mm y −1 .There is a clear significant influence of the infiltration amount.For each isotope, simulation time, and given infiltration amount, the ratio between the unsaturated zone (at selected level) and the saturated zone concentration is about the same.The measure of dilution in the unsaturated zone is influenced mainly by the infiltration amount; the influence of the simulation time and isotope is negligible.Figure 4 shows the dependence of 129 I concentration on the depth and infiltration amount at a simulation time of 50 000 y. Figure 5 shows the evolution of 129 I concentration at selected depths.

Case 3 -precipitation and evaporation based on meteorological data
The model domain is 10 m thick and discretized with a 0.1 m step.The soil is sandy loam with the following parameters: θ r = 0.065, θ s = 0.41, The saturated part of the model domain is 7 m thick (groundwater level 3 m below surface).This corresponds to a bottom flow boundary condition of a constant 7 m pressure head.Flow boundary conditions on the top of the model domain are derived from meteorological measurements from the past 200 years [19].The influx to the model domain is based on daily precipitation data.The evaporation is estimated based on mean temperatures using the Hamon model [11] and included in the simulation as a negative source uniformly distributed in the top 1 m of the model domain.
The transport simulation assumes that contamination with a concentration of 1 µg m −3 appears (as a        Sn, 129 I in selected simulation times, for various precipitation infiltration amounts in selected depths.consequence of the flow) in the previously clean saturated portion of the model domain during the first simulation time step.This concentration remains constant throughout the simulation.The simulated tracer is 129 I; it is non-sorbing and its decay is negligible (its half-life is much higher than the simulation period).The simulation period is dictated by the length of the available meteorological data, i.e. 78 965 days (May 1 st 1804 to December 31 st 2018).The simulation time step was 0.01 days.Figure 6 shows the time evolutions of selected quantities.They capture the first five years of the simulation (the time axis has a MM.YY format).From the results, we can see the initial rapid change to an equilibrium state followed by seasonal oscillations around a mean value.For example, in the case of the concentrations in the 2 m depth (Figure 6c), the oscillation magnitude is about 10 % of the mean, in the 0.5 m depth, it is 18 % and at surface level, (Figure 6d) it is up to 200 %.This significant near-surface concentration uptick happens during summer months when the evaporation rate is higher, which, combined with lower precipitation amounts, leads to higher fluxes from the groundwater level towards the surface.The significantly low values in the summer season are caused by rainstorms.Colder seasons with higher precipitation amounts (usually from October to March) are linked with lower concentrations and oscillations.The concentrations drop by 14-15 orders in the 3 m thick unsaturated zone; from 1 • 10 −9 in the saturated zone to around (depending on time) 6 • 10 −24 kg m −3 near the surface.
During the simulation of this case study, there were rapid changes of the flow boundary condition (precipitation).This has led to instabilities and non-  physical solutions.The stability was increased by shortening the simulation time step to the final value of 0.01 days.Another way to increase the stability would be to increase the portion of simulation domain in which the evaporation and root withdrawal takes place.The 1 m thickness was used.The third way (not used in this study) would be to average the precipitation/evaporation amounts over a longer time period.

Conclusion
For a simulation of flow and transport in an unsaturated zone, several concepts and their implementations (both commercial and non-commercial) exist.Despite that (for reasons stated in the introduction), we have decided to implement our own tool simulating the flow using the Richards' equation and transport (advection, diffusion, sorption and radioactive decay) using the advection-diffusion equation.The result is a SW that allows for a simulation of these processes in the unsaturated zone or, to be exact, in the simulation domain that consists of both the saturated and unsaturated zones.This SW was designed and implemented as a part of a more complex system that tracks the radionuclides released from a deep repository through the saturated zone, the unsaturated zone and the biosphere all the way to an individual (in the form of calculated effective dose).We have created a tool for a safety assessment of contamination sources, such as deep repositories of spent nuclear fuel and highly active wastes in a geosphere.Our implementation allowed us to have a precise control over the data exchange between individual models covering portions of the transport path.However, it forced us to deal with potentially unstable states and to verify the model so that we rule out any significant conceptual or implementation errors.
Three case studies were described in this article.They were designed to show the behaviour and results of the model in cases with varying boundary conditions and sources.Such cases are very challenging for existing SW tools used in the field and our ability to solve them is one of the main achievements of the SW concept (and implementation) shown in this article.Simulations with daily changes in precipitation and evaporation proved challenging because they can cause instabilities manifesting as negative concentrations or oscillations (not corresponding to periodicity in inputs).The presented cases prove the solvability of such situations by the presented model implementation as well as its usability within the system for a safety assessment of risks linked to biosphere contamination by radionuclides released from a deep repository.The system could also be used for cases when the contamination source is in the atmosphere and contaminates the environment through fallout and rain water.

2 s
and D is the diffusion coefficient [m 2 s −1 ].Numerical oscilations are negligible if P e < 5.The denominator in the equation for Péclet number computation includes the diffusion coefficient, which consists of hydrodynamical dispersion and molecular diffusion.In the case of 1D simulation, it means D W = D L |q| θ + D w τ w , where D L is the longitudal dispersivity [m], D w is the coefficient of molecular diffusion [m 2 s −1 ] and τ w = θ 7/3 θ is tortuosity.Péclet number may be lowered by including both dispersion and molecular diffusions in the simulation.

Figure 1 .
Figure 1.Pressure head time evolution at 5 m depth (left) and 1 m depth (right).

Figure 2 .
Figure 2. Concentration time evolution at 4 m depth (left) and 1 m depth (right).

Figure 3 .
Figure 3. Concentration evolutions at the geosphere/biosphere interface used as a boundary condition in the unsaturated zone transport simulations.

Figure 4 .
Figure 4. Dependence of 129 I concentration on depth and infiltration amount at simulation time 50 000 y.

Figure 5 .
Figure 5. Evolution of 129 I concentration at selected depths.

Figure 6 .
Figure 6.Time evolution of (from top to bottom) a) daily precipitation, b) evaporation, c) concentration 2 m below terrain and d) concentration at terrain level.

Table 1 .
Values of half-life, molecular diffusion coefficient and linear sorption distribution coefficient of selected isotopes used in the simulation.This case deals with the transport of 10 real isotopes ( 14 C, 36 Cl, 41 Ca,59Ni, 79 Se, 107 Pd, 126 Sn, 129 I, 135 Cs and 238 U) which are expected to migrate from a deep repository.The isotope parameters used in the simulation are shown in Table