Acta Polytechnica https://doi.org/10.14311/AP.2022.62.0427 Acta Polytechnica 62(4):427–437, 2022 © 2022 The Author(s). Licensed under a CC-BY 4.0 licence Published by the Czech Technical University in Prague A MODEL OF ISOTOPE TRANSPORT IN THE UNSATURATED ZONE, CASE STUDY Josef Chudoba∗, Jiřina Královcová, Jiří Landa, Jiří Maryška, Jakub Říha Technical University of Liberec, Faculty of Mechatronics, Informatics and Interdisciplinary Studies, Studentská 2, Liberec, Czech Republic ∗ corresponding author: josef.chudoba@tul.cz Abstract. This work deals with a groundwater flow and solute transport model in the near-surface (predomi- nantly 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. Keywords: Richards’ equation, unsaturated zone, groundwater flow in unsaturated zone, solute transport. 1. 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–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 en- vironmental contamination. Existing SW tools (Pan- dora [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 concen- trations 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 un- saturated zone flow simulation is based on Richards’ 427 https://doi.org/10.14311/AP.2022.62.0427 https://creativecommons.org/licenses/by/4.0/ https://www.cvut.cz/en www.skb.se www.possiva.fi J. Chudoba, J. Královcová, J. Landa et al. Acta Polytechnica 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 struc- ture for transport simulations in both saturated and unsaturated zones, which allows us to avoid difficult interconnection of inputs and outputs of various sim- ulation 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]. 2. 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 time- space 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 avail- able. It is possible to identify the near-surface zones of the model domain for which the more detailed unsatu- rated 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 unsat- urated 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 mea- surement 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 con- dition is usually zero concentration (when we assume the source of contamination to be located deep in the geosphere). The bottom boundary condition is de- fined based on the concentration values acquired from a simulation of the saturated zone. This boundary condition interconnects the saturated and the unsatu- rated 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. 3. Mathematical and physical model 3.1. 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. 428 vol. 62 no. 4/2022 A model of isotope transport in the unsaturated zone, case study The flow in the unsaturated zone is described by Richards’ equation [14]: ∂θ(h) ∂t = ∂ ∂z ( K(θ, h) ( ∂h ∂z + 1 )) − S, (1) where h [m] is the pressure head, θ [1] is the water con- tent, 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 hy- draulic 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: θ(h) =  θr + θs − θr [1 + |αh|n]m , h < 0 θs, h ≥ 0 , (2) K(h) = KS · S0.5e · [ 1 − ( 1 − S1/me )m]2 , (3) Se = θ(h) − θr θs − θr , m = 1 − 1 n , 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 Ks [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 param- eters. 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 simula- tion, both their type and their value. In each time step, it is necessary for the pressure head to be speci- fied 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 pos- sible 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. 3.2. Model of transport In the unsaturated zone, the contamination may gen- erally exist dissolved in three phases – liquid, solid and gaseous. In our concept, we do not account for a con- tamination 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 contami- nants/isotopes is very low, which is why the model uses only the linear sorption isotherm. The transport model also includes diffusion and radioac- tive 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]. ∂θck ∂t + ∂ρsk ∂t = ∂ ∂z ( θDwk ∂ck ∂z ) − ∂qck ∂z − µw,kθck − µs,kρsk + n∑ m=1 m ̸=k µw,mθcm+ n∑ m=1 m̸=k µs,mρsm + γw,kθ + γs,kρ − rk, (4) where θ [1] is the water content, ρ [kg m−3] is the rock density, z [m] is the vertical axis, ck [kg m−3] is the concentration of isotope k in a liquid phase, sk [1] is the concentration of isotope k in a solid phase, Dwk [m 2 s−1] is the diffusion coefficient of isotope k, q [m s−1] is the flux, µw,k and µs,k [s−1] is the radioac- tive 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, rk [kg m−3 s−1] are sources. The relationship between concentrations in solid phase sk and in liquid phase ck is given by equation (5). sk = kD,k · ck, (5) where kD,k [m3 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. 4. 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 implemen- tation. Potential instabilities must be eliminated by 429 J. Chudoba, J. Královcová, J. Landa et al. Acta Polytechnica an appropriate choice of parameter values. This sec- tion 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 evapora- tion 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 ground- water 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 Cri (here, the lower index i stands for the element of dis- cretisation) 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 con- ditions 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 gen- erally 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 os- cillations 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] and D is the diffusion coefficient [m2 s−1]. Numerical oscilations are negligi- ble if P e < 5. The denominator in the equation for Péclet number computation includes the diffusion co- efficient, which consists of hydrodynamical dispersion and molecular diffusion. In the case of 1D simula- tion, it means DW = DL|q| θ + Dwτw, where DL is the longitudal dispersivity [m], Dw is the coefficient of molecular diffusion [m2 s−1] and τw = θ 7/3 θ2s is tortuos- ity. Péclet number may be lowered by including both dispersion and molecular diffusions in the simulation. Oscillations of the solution and/or negative concen- trations may also be caused by changing concen- tration 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 precip- itation 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 do- main so that it captures a greater part of the saturated zone. 5. Case studies This part of the article shows three case studies. The first one is a synthetic task that includes a simula- tion 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 radionu- clide migration in the case of a realistic precipitation boundary condition (data from Praha Klementinum; available since 1804). This case also includes the evap- oration 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, de- rived 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). 5.1. 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 soil properties are constant throughout the model domain. The soil is 40 % silt, 15 % clay and 45 % sand. Its dry density is 1500 kg m−3. This composition corresponds to the following parameters: KS = 1.912 · 10−6 m s−1, n = 1.469, θr = 0.0492, θs = 0.3687, a = 1.355 m−1. The flux, representing precipitation (after subtrac- tion of evaporation), is used as a flow boundary con- dition on the top boundary. The time evolution of 430 vol. 62 no. 4/2022 A model of isotope transport in the unsaturated zone, case study Figure 1. Pressure head time evolution at 5 m depth (left) and 1 m depth (right). Figure 2. Concentration time evolution at 4 m depth (left) and 1 m depth (right). the boundary condition has a sinus shape with a 365- day periodicity. Its values for individual simulation times were (for the expected influx of 365 mm y−1) calculated as okppF 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 dis- tribution 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 Fig- ure 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.17 m (difference of ca 1.44 m). Based on the flow simulation results, the transport of two isotopes (noted isotope 1 and isotope 2) was simulated. Isotope 1 has a half-life of 10 000 days, its product is the stable isotope 2. The following parame- ters were used for the simulation: coefficient of molecu- lar diffusion 1 · 10−3 m2 day−1 (i.e. 1.15 · 10−8 m2 s−1), longitudal dispersivity 0.1 m, distribution coefficient 0 m3 kg−1 for isotope 1 and 0.0001 m3 kg−1 [18] for isotope 2. 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 horizon- tal 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 concen- tration 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, pre- cipitation and evaporation amounts, roots, and wells. 431 J. Chudoba, J. Královcová, J. Landa et al. Acta Polytechnica Isotope Half-life [y] Diffusion coefficient [m2 s−1] Distribution coefficient [m3 kg−1] 14 C 5 700 1e-9 0.0005 36 Cl 301 000 1e-9 0.0 41 Ca 102 000 1e-9 0.0001 59 Ni 76 000 1e-9 0.01 79 Se 356 000 1e-9 0.0005 107 Pd 6 500 000 1e-9 0.0001 126 Sn 230 000 1e-9 0.0 129 I 16 100 000 1e-9 0.0 135 Cs 2 300 000 1e-9 0.01 238 U 4 468 000 000 1e-9 0.1 Table 1. Values of half-life, molecular diffusion coefficient and linear sorption distribution coefficient of selected isotopes used in the simulation. 5.2. Case 2 – radionuclide transport with realistic transport boundary condition, sensitivity on infiltration amount This case deals with the transport of 10 real isotopes (14 C, 36 Cl, 41 Ca, 59 Ni, 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 1. The regional satu- rated model transport simulation period was 500 000 years. As a consequence of a given radionuclide re- lease 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 con- centration 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, α = 7.5 m−1, n = 1.89, Ks = 1.228 · 10−5 ms−1. The saturated part of the model domain is 8 m thick (groundwater level 2 m below surface). This corre- sponds to bottom flow boundary condition of a con- stant 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 concentra- tion 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. 5.3. 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, α = 7.5 m−1, n = 1.89, Ks = 1.228 · 10−5 m s−1. 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 con- ditions 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 esti- mated based on mean temperatures using the Hamon model [11] and included in the simulation as a neg- ative source uniformly distributed in the top 1 m of the model domain. The transport simulation assumes that contamina- tion with a concentration of 1 µg m−3 appears (as a 432 vol. 62 no. 4/2022 A model of isotope transport in the unsaturated zone, case study C l3 6 50 00 0 y 10 0 00 0 y 20 0 00 0 y pr ec ip it at io n [m m y− 1 ] 0. 3 m 0. 5 m 1. 0 m sa tu ra te d zo ne 0. 3 m 0. 5 m 1. 0 m sa tu ra te d zo ne 0. 3 m 0. 5 m 1. 0 m sa tu ra te d zo ne 0 1. 5E -0 6 1. 6E -0 6 1. 6E -0 6 1. 6E -0 6 8. 3E -0 7 8. 4E -0 7 8. 6E -0 7 8. 8E -0 7 1. 6E -0 7 1. 6E -0 7 1. 6E -0 7 1. 7E -0 7 20 1. 1E -1 2 8. 2E -1 2 1. 1E -0 9 1. 6E -0 6 6. 1E -1 3 4. 5E -1 2 6. 1E -1 0 8. 8E -0 7 1. 2E -1 3 8. 4E -1 3 1. 2E -1 0 1. 7E -0 7 50 2. 0E -1 3 1. 6E -1 2 2. 8E -1 0 1. 6E -0 6 1. 1E -1 3 8. 4E -1 3 1. 5E -1 0 8. 8E -0 7 2. 0E -1 4 1. 6E -1 3 2. 8E -1 1 1. 7E -0 7 C a4 1 50 00 0 y 10 0 00 0 y 20 0 00 0 y pr ec ip it at io n [m m y− 1 ] 0. 3 m 0. 5 m 1. 0 m sa tu ra te d zo ne 0. 3 m 0. 5 m 1. 0 m sa tu ra te d zo ne 0. 3 m 0. 5 m 1. 0 m sa tu ra te d zo ne 0 8. 5E -1 9 8. 6E -1 9 9. 3E -1 9 1. 0E -1 8 8. 5E -1 9 8. 6E -1 9 9. 3E -1 9 1. 0E -1 8 8. 5E -1 9 8. 6E -1 9 9. 3E -1 9 1. 0E -1 8 20 7. 0E -2 5 5. 1E -2 4 7. 0E -2 2 1. 0E -1 8 7. 0E -2 5 5. 1E -2 4 7. 0E -2 2 1. 0E -1 8 7. 0E -2 5 5. 1E -2 4 7. 0E -2 2 1. 0E -1 8 50 1. 2E -2 5 9. 5E -2 5 1. 7E -2 2 1. 0E -1 8 1. 2E -2 5 9. 5E -2 5 1. 7E -2 2 1. 0E -1 8 1. 2E -2 5 9. 5E -2 5 1. 7E -2 2 1. 0E -1 8 S n 12 6 50 00 0 y 10 0 00 0 y 20 0 00 0 y pr ec ip it at io n [m m y− 1 ] 0. 3 m 0. 5 m 1. 0 m sa tu ra te d zo ne 0. 3 m 0. 5 m 1. 0 m sa tu ra te d zo ne 0. 3 m 0. 5 m 1. 0 m sa tu ra te d zo ne 0 9. 3E -1 9 9. 3E -1 9 9. 7E -1 9 1. 0E -1 8 6. 9E -1 8 7. 1E -1 8 7. 5E -1 8 7. 9E -1 8 1. 4E -1 4 1. 4E -1 4 1. 4E -1 4 1. 5E -1 4 20 7. 0E -2 5 5. 1E -2 4 7. 0E -2 2 1. 0E -1 8 5. 5E -2 4 4. 0E -2 3 5. 5E -2 1 7. 9E -1 8 1. 0E -2 0 7. 5E -2 0 1. 0E -1 4 1. 5E -1 4 50 1. 2E -2 5 9. 5E -2 5 1. 7E -2 2 1. 0E -1 8 9. 5E -2 5 7. 5E -2 4 1. 3E -2 1 7. 9E -1 8 1. 8E -2 1 1. 4E -2 0 2. 5E -1 8 1. 5E -1 4 I1 29 50 00 0 y 10 0 00 0 y 20 0 00 0 y pr ec ip it at io n [m m y− 1 ] 0. 3 m 0. 5 m 1. 0 m sa tu ra te d zo ne 0. 3 m 0. 5 m 1. 0 m sa tu ra te d zo ne 0. 3 m 0. 5 m 1. 0 m sa tu ra te d zo ne 0 1. 6E -0 5 1. 6E -0 5 1. 6E -0 5 1. 5E -0 5 8. 1E -0 6 8. 1E -0 6 8. 1E -0 6 8. 1E -0 6 1. 3E -0 6 1. 3E -0 6 1. 3E -0 6 1. 3E -0 6 20 1. 1E -1 1 7. 8E -1 1 1. 1E -0 8 1. 5E -0 5 5. 6E -1 2 4. 1E -1 1 5. 6E -0 9 8. 1E -0 6 9. 1E -1 3 6. 6E -1 2 9. 1E -1 0 1. 3E -0 6 50 1. 9E -1 2 1. 5E -1 1 2. 6E -0 9 1. 5E -0 5 9. 7E -1 3 7. 7E -1 2 1. 4E -0 9 8. 1E -0 6 1. 6e -1 3 1. 2E -1 2 2. 2E -1 0 1. 3E -0 6 T ab le 2. C on ce nt ra ti on s of 36 C l, 41 C a, 12 6 Sn , 12 9 I in se le ct ed si m ul at io n ti m es , fo r va ri ou s pr ec ip it at io n in fil tr at io n am ou nt s in se le ct ed de pt hs . 433 J. Chudoba, J. Královcová, J. Landa et al. Acta Polytechnica Figure 3. Concentration evolutions at the geosphere/biosphere interface used as a boundary condition in the unsaturated zone transport simulations. Figure 4. Dependence of 129 I concentration on depth and infiltration amount at simulation time 50 000 y. consequence of the flow) in the previously clean sat- urated portion of the model domain during the first simulation time step. This concentration remains con- stant 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 1st 1804 to December 31st 2018). The simula- tion 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 val- ues in the summer season are caused by rainstorms. Colder seasons with higher precipitation amounts (usu- ally 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 (pre- cipitation). This has led to instabilities and non- 434 vol. 62 no. 4/2022 A model of isotope transport in the unsaturated zone, case study Figure 5. Evolution of 129 I concentration at selected depths. 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. 435 J. Chudoba, J. Královcová, J. Landa et al. Acta Polytechnica 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 pre- cipitation/evaporation amounts over a longer time period. 6. Conclusion For a simulation of flow and transport in an unsat- urated zone, several concepts and their implementa- tions (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 sim- ulation 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 repos- itory 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 contamina- tion 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 re- sults of the model in cases with varying boundary conditions and sources. Such cases are very challeng- ing 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 precip- itation 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 sys- tem 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. Acknowledgements This article was prepared with a support from the Technol- ogy Agency of the Czech Republic within the EPSILON program. www.tacr.cz. References [1] P. O. Johansson. Description of surface hydrology and near-surface hydrogeology at Forsmark. Site descriptive modelling, SDM-Site Forsmark. Tech. Rep. SKB R-08-08, Swedish Nuclear Fuel and Waste Management Co, Stockholm, Sweden, 2008. https: //www.skb.com/publication/1853247/R-08-08.pdf. [2] K. Werner, M. Sassner, E. Johansson. Hydrology and near-surfacehydrogeology at Forsmark – synthesis for the SR-PSU project. SR-PSU Biosphere. Tech. Rep. SKB R-13-19, Swedish Nuclear Fuel and Waste Management Co, Stockholm, Sweden, 2013. https: //www.skb.com/publication/2477948/R-13-19.pdf. [3] E. Jutebring Sterte, E. Johansson, Y. Sjöberg, et al. Groundwater-surface water interactions across scales in a boreal landscape investigated using a numerical modelling approach. Journal of Hydrology 560:184–201, 2018. https://doi.org/10.1016/j.jhydrol.2018.03.011. [4] P. A. Ekström. Pandora – a simulation tool for safety assessments. Technical description and user’s guide. Tech. Rep. SKB R-11-01, Swedish Nuclear Fuel and Waste Management Co, Stockholm, Sweden, 2010. https: //www.skb.com/publication/2188873/R-11-01.pdf. [5] C. Yu, E. Gnanapragasam, J.-J. Cheng, et al. User’s Manual for RESRAD-OFFSITE Code Version 4. United States Nuclear Regulatory Commission, Argonne National Laboratory, USA, 2020, https://resrad.evs. anl.gov/docs/RESRAD-OFFSITE_UsersManual.pdf. [6] J. Šimůnek, M. Šejna, H. Saito, et al. The Hydrus-1D Software Package for Simulating the Movement of Water, Heat, and Multiple Solutes in Variably Saturated Media, Version 4.17, HYDRUS Software Series 3. Department of Environmental Sciences, University of California Riverside, Riverside, California, USA, 2013, https://www.pc-progress.com/ Downloads/Pgm_hydrus1D/HYDRUS1D-4.17.pdf. [7] C. I. Steefel, C. A. J. Appelo, B. Arora, et al. Reactive transport codes for subsurface environmental simulation. Computational Geosciences 19(3):445–478, 2015. https://doi.org/10.1007/s10596-014-9443-x. [8] J. Březina, M. Hokr. Mixed-hybrid formulation of multidimensional fracture flow. In I. Dimov, S. Dimova, N. Kolkovska (eds.), Numerical Methods and Applications. NMA 2010. Lecture Notes in Computer Science, vol. 6046. Springer, Berlin, Heidelberg, 2011. https://doi.org/10.1007/978-3-642-18466-6_14. [9] J. Březina, J. Stebel, D. Flanderka, et al. Flow123d, version 2.2.1. User Guide and Input Reference. Technical university of Liberec, Faculty of mechatronics, informatics and interdisciplinary studies, Liberec, 2018, https://flow.nti.tul.cz/packages/2.2.1_release/ flow123d_2.2.1_doc.pdf. [10] M. Hokr, H. Shao, W. P. Gardner, et al. Real-case benchmark for flow and tracer transport in the fractured rock. Environmental Earth Sciences 75(18):1273, 2016. https://doi.org/10.1007/s12665-016-6061-z. 436 www.tacr.cz https://www.skb.com/publication/1853247/R-08-08.pdf https://www.skb.com/publication/1853247/R-08-08.pdf https://www.skb.com/publication/2477948/R-13-19.pdf https://www.skb.com/publication/2477948/R-13-19.pdf https://doi.org/10.1016/j.jhydrol.2018.03.011 https://www.skb.com/publication/2188873/R-11-01.pdf https://www.skb.com/publication/2188873/R-11-01.pdf https://resrad.evs.anl.gov/docs/RESRAD-OFFSITE_UsersManual.pdf https://resrad.evs.anl.gov/docs/RESRAD-OFFSITE_UsersManual.pdf https://www.pc-progress.com/Downloads/Pgm_hydrus1D/HYDRUS1D-4.17.pdf https://www.pc-progress.com/Downloads/Pgm_hydrus1D/HYDRUS1D-4.17.pdf https://doi.org/10.1007/s10596-014-9443-x https://doi.org/10.1007/978-3-642-18466-6_14 https://flow.nti.tul.cz/packages/2.2.1_release/flow123d_2.2.1_doc.pdf https://flow.nti.tul.cz/packages/2.2.1_release/flow123d_2.2.1_doc.pdf https://doi.org/10.1007/s12665-016-6061-z vol. 62 no. 4/2022 A model of isotope transport in the unsaturated zone, case study [11] W. R. Hamon. Estimating potential evapotranspiration. Journal of the Hydraulics Division 87(3):107–120, 1961. https://doi.org/10.1061/JYCEAJ.0000599. [12] T. Vogel, K. Huang, R. Zhang, M. T. van Genuchten. The HYDRUS code for simulating one-dimensional water flow, solute transport, and heat movement in variably-saturated media, Version 5.0. Research Report No 140. U.S. Salinity Laboratory, USDA, ARS, Riverside, CA, 1996. [13] K. Zhang, Y.-S. Wu, J. E. Houseworth. Sensitivity analysis of hydrological parameters in modeling flow and transport in the unsaturated zone of Yucca Mountain, Nevada, USA. Hydrogeology Journal 14(8):1599–1619, 2006. https://doi.org/10.1007/s10040-006-0055-y. [14] L. A. Richards. Cappillary conduction of liquids through porous mediums. Physics 1(5):318–333, 1931. https://doi.org/10.1063/1.1745010. [15] M. T. van Genuchten. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal 44(5):892–898, 1980. https://doi.org/10.2136/ sssaj1980.03615995004400050002x. [16] M. G. Schaap, F. J. Leij, M. T. van Genuchten. rosetta: a computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions. Journal of Hydrology 251(3-4):163–176, 2001. https://doi.org/10.1016/S0022-1694(01)00466-8. [17] M. A. Celia, E. T. Bouloutas, R. L. Zarba. A general mass-conservative numerical solution for the unsaturated flow equation. Water Resources Research 26(7):1483–1496, 1990. https://doi.org/10.1029/wr026i007p01483. [18] E. Juranová, E. Hanslík, S. Dulanská, et al. Sorption of anthropogenic radionuclides onto river sediments and suspended solids: dependence on sediment composition. Journal of Radioanalytical and Nuclear Chemistry 324(3):983–991, 2020. https://doi.org/10.1007/s10967-020-07174-w. [19] Historical data: Weather: Prague Clementinum. CHMI portal. [2020-11-02], http://portal.chmi.cz/ historicka-data/pocasi/praha-klementinum?l=en#. 437 https://doi.org/10.1061/JYCEAJ.0000599 https://doi.org/10.1007/s10040-006-0055-y https://doi.org/10.1063/1.1745010 https://doi.org/10.2136/sssaj1980.03615995004400050002x https://doi.org/10.2136/sssaj1980.03615995004400050002x https://doi.org/10.1016/S0022-1694(01)00466-8 https://doi.org/10.1029/wr026i007p01483 https://doi.org/10.1007/s10967-020-07174-w http://portal.chmi.cz/historicka-data/pocasi/praha-klementinum?l=en# http://portal.chmi.cz/historicka-data/pocasi/praha-klementinum?l=en# Acta Polytechnica 62(4):427–437, 2022 1 Introduction 2 Model concept 3 Mathematical and physical model 3.1 Model of flow 3.2 Model of transport 4 Model instabilities 5 Case studies 5.1 Case 1 – periodic precipitation evolution 5.2 Case 2 – radionuclide transport with realistic transport boundary condition, sensitivity on infiltration amount 5.3 Case 3 – precipitation and evaporation based on meteorological data 6 Conclusion Acknowledgements References