Calibration of Hydrodynamic Models of Groundwater Flow in Aquifers Polluted by Liquid Radioactive Waste

: Pollution of aquifers by liquid radioactive waste is especially hazardous due to high permeability of aquifer rocks. Pollution can take place as a result of deep injection disposal of liquid waste or due to its leakage from a surface reservoir. The ecological hazard depends on distribution of aquifer permeability, which can be highly heterogeneous. The only source of data for determination of such distribution is measurements in exploratory boreholes at the polluted site. A technique is proposed for determination of spatial distribution of aquifer transmissibility on the basis of hydrological measurements in a limited number of boreholes. The problem is reduced to a minimization problem which was solved numerically.


Introduction
As it follows from some estimations [1,2], operation of all atomic power plants over the world leads to production of ∼ =11,300-11,500 t of spent nuclear fuel (SNF) per year. With a fuel cycle involving reprocessing and potential recycling of plutonium, SNF is dispatched after extraction from a reactor to radiochemical plants for reprocessing. Processing of 1 t of SNF leads to production of 13-34 m 3 of liquid radioactive waste (LRW) [3] with activity of more than 10 Ci/l [4].
In the USSR and later in Russia, deep injection of LRW was carried out on an industrial scale. LRW was injected into deep-seated, highly permeable aquifers confined from top and bottom by low permeable rocks. Such injection technology was implemented for disposal of LRW from the Mining and Chemical Plant in the Krashoyarsk region [5] (Figure 1).
Reprocessing of SNF is not the only source of LRW, which can be produced also from maintenance of operation of atomic power plants. This LRW is usually accumulated at territories of the power stations in special reservoirs. As a result of an emergence, LRW can gush out from the reservoir to the earth's surface and then move with groundwater in the upper free aquifer to the zone of groundwater discharge, to an open reservoir or a river. An example of such radioactive pollution was an accident at the atomic power plant (APP) Novovoronezhskaya in 1985 [6].
In both cases, a radioactive plume moves from a pollution source to the biosphere through a permeable sedimentary rock layer thickness which is much less than its areal dimensions.
Assessment of an ecological hazard of radionuclides released to the biosphere calls for calculation of groundwater stream velocities. The key parameter in this calculation is spatially distributed water transmissivity of the aquifer rocks [7]. The only source of input data for determination of this distribution is data of transmissivity measurements in observation boreholes which are drilled at the pollution site. Since drilling of such boreholes is rather expensive (especially in a case of a deep-seated aquifer), the amount of Assessment of an ecological hazard of radionuclides released to the biosphere ca for calculation of groundwater stream velocities. The key parameter in this calculation spatially distributed water transmissivity of the aquifer rocks [7]. The only source of inp data for determination of this distribution is data of transmissivity measurements in o servation boreholes which are drilled at the pollution site. Since drilling of such borehol is rather expensive (especially in a case of a deep-seated aquifer), the amount of measur ment data is limited. A calibration of a mathematical model of groundwater flow in th aquifer is a determination of such spatial distribution of water transmissivity of the aqu fer, which is consistent with hydrological measurements in the observation boreholes.
The deep injection site for disposal of LRW from MCP is located on the eastern ban of the river Yenisei in the Krasnoyarsk region ( Figure 2). The geological medium in th injection domain is represented by a small-size sedimentary basin. It is built up by a la ered clay-sandy Jurassic sequence with the maximum thickness of 550 m, which fills a erosional paleovalley in the Pre-Cambrian crystalline basement. Three (sandy to weak cemented sandstone) aquifers are relatively clearly defined in the geological section of th sedimentary sequence ( Figure 3) and are labeled from bottom to top as Horizons I, II, an III. The lowest aquifer (Horizon I) is confined from bottom by low-permeability baseme rocks.
The deep injection site for disposal of LRW from MCP is located on the eastern bank of the river Yenisei in the Krasnoyarsk region ( Figure 2). The geological medium in the injection domain is represented by a small-size sedimentary basin. It is built up by a layered clay-sandy Jurassic sequence with the maximum thickness of 550 m, which fills an erosional paleovalley in the Pre-Cambrian crystalline basement. Three (sandy to weakly cemented sandstone) aquifers are relatively clearly defined in the geological section of the sedimentary sequence ( Figure 3) and are labeled from bottom to top as Horizons I, II, and III. The lowest aquifer (Horizon I) is confined from bottom by low-permeability basement rocks.
The layered sequence is truncated both structurally and hydraulically by a subvertical tectonic fault zone (right-bank fault) (Figure 3a). This fault is traced for a distance of about 2 km to the north and south of the sanitary-protection zone. Uplift of the block of the basement rock to the west of the fault seals Horizon I from the west (Figure 3c) and separates the waste plume in Horizon I from the Yenisei river.
Intermediate-level LRW (with radioactivity up to 10 −2 Ci/liter) are injected into Horizon I. From 1972 to 1993, Horizon I has been used also for disposal of LRW with elevated radioactivity (up to 5 Ci/l) [5]. Processes 2022, 10, x FOR PEER REVIEW 3 of 13

Pollution Site at Atomic Power Plant Novovoronezhskaya
Pollution of groundwater with 60 Co at the site of APP Novoronezhskaya is caused by a leakage of concentrated technological LRW from a temporary storage in 1985. Estimated volume of the leakage was 480 m 3 , and activity of 60 Co in LRW at the time of the leakage was 1.34 × 10 11 Bq/m 3 [6].
APP and its temporary storage of LWR are located on the left bank of the river Don ( Figure 4a). The upper part of the underground medium at the site of APP is represented by a layered sequence of clayey-sandy sedimentary deposits confined from the bottom by tight basement rocks (Figure 4b). The contaminant plume moves within two upper aquifers above the argillaceous layer which is denoted with letter A in Figure 4b. As it is seen in (Figure 4b) they are partially divided from each other by a clayey interlayer, however they have a structural and hydrologic connection, as is shown in Figure 4b. An important argument in favor of this statement is very smooth relief of the water table (Figure 4b).

Confined Aquifer Example of the Injection Site of LRW from Mining and Chemicals Combined
If thickness of an aquifer is much less than its areal dimensions (it should be taken into account that horizontal scales in Figures 3b and c differ approximately by an order of magnitude from the vertical scale), a 2D planar approximation is quite sufficient for determination of groundwater flow in the aquifer. In a case of a confined aquifer, measurement data include values of hydraulic head in the observation boreholes and values of aquifer transmissibility, obtained by pumping tests in a part of the boreholes ( Figure 5). The calibration problem reduces to determination of such spatial distribution of transmissibility that measured transmissibility coincides with the corresponding local value of the distribution, and values of hydraulic head, which are calculated on the basis of this distribution, and are in good agreement with the measured values.

Confined Aquifer Example of the Injection Site of LRW from Mining and Chemicals Combined
If thickness of an aquifer is much less than its areal dimensions (it should be taken into account that horizontal scales in Figure 3b,c differ approximately by an order of magnitude from the vertical scale), a 2D planar approximation is quite sufficient for determination of groundwater flow in the aquifer. In a case of a confined aquifer, measurement data include values of hydraulic head in the observation boreholes and values of aquifer transmissibility, obtained by pumping tests in a part of the boreholes ( Figure 5). The calibration problem reduces to determination of such spatial distribution of transmissibility that measured transmissibility coincides with the corresponding local value of the distribution, and values of hydraulic head, which are calculated on the basis of this distribution, and are in good agreement with the measured values. Let us consider an areal model of groundwater flow in Horizon I, taking into account a leakage between Horizons I and II. Pumping tests and hydraulic head measurements were carried before injection. Salinity of the groundwater was uniform at that time. Hence, buoyancy forces caused by heterogeneity of the groundwater density were absent at that time. Then velocity components satisfied Darcy's law in the form is volume of groundwater which flows through unit square of the upper boundary of Horizon I per unit time (positive at ascending leakage).
According to [7], the value of q can be represented as Let us consider an areal model of groundwater flow in Horizon I, taking into account a leakage between Horizons I and II. Pumping tests and hydraulic head measurements were carried before injection. Salinity of the groundwater was uniform at that time. Hence, buoyancy forces caused by heterogeneity of the groundwater density were absent at that time. Then velocity components satisfied Darcy's law in the form where f = f (x, y) is hydraulic conductivity of rock in Horizon I, and H 1 = H 1 (x, y) is hydraulic head in Horizon I. The equation of water balance in Horizon I can be written in the form where δ = δ(x, y) is thickness of Horizon I; q = q(x, y) is volume of groundwater which flows through unit square of the upper boundary of Horizon I per unit time (positive at ascending leakage). According to [7], the value of q can be represented as where H 2 = H 2 (x, y) is hydraulic head in Horizon II; σ = σ(x, y) is parameter of the leakage between Horizons I and II.
Here, f tr = f tr (x, y) is hydraulic conductivity of the low-permeability clayey layer between Horizons I and II; ∆ = ∆(x, y) is thickness of the layer. Substitution of (1) and (3) into (2) gives where W = δ f is transmissibility of Horizon I, and Ω is modeling domain. Let us specify Ω as it is shown with a yellow dashed line in Figure 6.  (1) and (3) into (2) gives  is transmissibility of Horizon I, and  is modeling domain. Let us specify  as it is shown with a yellow dashed line in Figure 6. As it follows from Figure 3b, the right-bank fault forms an impermeable boundary for the part of Horizon I to the east of the fault. The zone of wedging out of Horizon I is also impermeable. Hence, we can specify zero Neumann conditions at the eastern and western boundaries of  . As it follows from measurements of hydraulic head in observation boreholes, the gradient of hydraulic head at the northern boundary is low ( Figure  6). This allows consideration of two practically identical variants of boundary conditions at the northern boundary of  :  As it follows from Figure 3b, the right-bank fault forms an impermeable boundary for the part of Horizon I to the east of the fault. The zone of wedging out of Horizon I is also impermeable. Hence, we can specify zero Neumann conditions at the eastern and western boundaries of Ω. As it follows from measurements of hydraulic head in observation boreholes, the gradient of hydraulic head at the northern boundary is low ( Figure 6). This allows consideration of two practically identical variants of boundary conditions at the northern boundary of Ω:
Zero Neumann conditions.
Values of H 1 at the southern boundary of Ω are determined by piece-linear interpolation of values, which are measured in observation boreholes.
M values of transmissibility in M points of Ω: Let us denote {W} an aggregate of distributions, W(X, Y) which satisfies conditions (6), and corresponding solutions of (4) satisfy conditions (5) at specified H 2 and conditions at the boundary of Ω. It is possible that such a distribution is not unique, and {W} consists of several elements (or even of infinite number of them). Hence, a problem emerges to specify an additional condition for selection of a unique element from {W}.
A characteristic feature of aquifers from unconsolidated sediments in layered sequences is relatively smooth variation of their hydraulic transmissivity on macroscale in horizontal directions, in comparison with hydraulic conductivity variations on a microscale [8]. As a condition for selection of a unique element from {W}, this can be written in the form Therefore, the calibration problem is reduced to a search of a conditional minimum of the functional Φ(W). In the considered case, the problem can be solved only numerically.
Independently of the method of numerical solution of Equation (4), it is desirable from the viewpoint of accuracy to use a grid where points {(x H i , y H i ), i = 1, . . . , N} and {(x W j , y W j ), j = 1, . . . , M} are nodes. As it is seen in Figure 5, points of boreholes are distributed in an irregular manner, and use of finite-difference methods with a rectangular grid is complicated. That is why we used a finite-element method [9].
Let us denote the vector of nodal values of H 1 as H 1 . If we use a weak Galerkin finite-element method [9,10], the boundary problem for Equation (4) takes the form A is a symmetric positively definite matrix what allows the solving of Equation (8) with the fast Cholesky method [11].
With use of finite-element approximation, the distribution W is fully determined through a vector of nodal values of W, and the functional Φ takes the form where W is vector of nodal values, W i , and S is a symmetric positively defined matrix.
Conditions (5) and (6) where I n , J m are a numbers of nodes corresponding to boreholes with numbers n and m, respectively. Hence, the calibration of the model of groundwater flow reduces to the search of a maximum of the quadratic form (9) at conditions (8) and (10), i.e., to a search for a conditional minimum in classic formulation.
Direct methods of search for an extremum with implicitly defined conditions are not obvious. That is why the problem was solved by an approximate method with use of a penalty function [12], and we determined an unconditional minimum of the functional where G is a parameter. The calculations were carried out at increasing G, up to such a value that the variation of W at future growth of G does not exceed 5%. Since points of boreholes {(x W j , y W j ), j = 1, . . . , M} are nodes of finite elements, values of W(x W j , y W j ) are fixed and equal to measured values. Hence, the second sum in the square brackets in the right-hand side of the expression (11) is equal to zero.
A minimum of (11) was determined by a modified method of gradient minimization. Obtained distributions of W and σ are presented in Figure 7.
Direct methods of search for an extremum with implicitly defined conditions are not obvious. That is why the problem was solved by an approximate method with use of a penalty function [12], and we determined an unconditional minimum of the functional (11) where G is a parameter. The calculations were carried out at increasing G , up to such a value that the variation of W at future growth of G does not exceed 5%. are fixed and equal to measured values. Hence, the second sum in the square brackets in the right-hand side of the expression (11) is equal to zero. A minimum of (11) was determined by a modified method of gradient minimization.

Since points of boreholes
Obtained distributions of W and  are presented in Figure 7.  Table 1.
One can see that calculated and measured values of 1 H are in a satisfactory agreement.  Calculated values of W at the points of boreholes are equal exactly to measured values. A comparison of calculated and measured values of H 1 is presented in Table 1. One can see that calculated and measured values of H 1 are in a satisfactory agreement.

Free Aquifer Example of the Pollution Site at APP Novovoronezhskaya
The aquifer is confined from its bottom by practically impermeable argillaceous deposits. Atmospheric precipitations enter the aquifer through its upper free boundary. The balance equation in this case takes the form where h = h(x, y) is thickness of the water-saturated layer in the aquifer, and w l and w u are volume of moisture which enters the aquifer through unit square per unit time from the lower confined aquifer and atmospheric precipitation, respectively. Substitution of (1) into (12) gives The water influx from above is caused by atmospheric precipitations. Estimation of them is based on standard data of meteorological service from the station 27,947 presented in [13]. A plot of the atmospheric precipitation rate against time is presented in Figure 8. According to recommended values from [14], parameters ered case were specified as cr w = 37 mm and I  = 0.88. F Figure 8. Atmospheric precipitations at pollution site after [13].
Only a part of precipitation moisture enters the groundwater [7,14]. A fraction of precipitation is lost by runoff and evapotranspiration. It depends on relief, type of soil and plant cover, and presence of asphalt covering and buildings. It can be seen in Figure 4a that the area covered by the asphaltic surface is not significant. Estimation of the fraction of atmospheric precipitations which enter the modeling domain is based on a simple model presented in [14]. According to this model, water infiltration rate (without evaporation) can be estimated as w I is volume of infiltration which enters the underground medium per day through the unit open surface; w A is volume of precipitation moisture which enters the unit surface per day; w cr is critical rate of daily precipitations; α I is a dimensionless coefficient.
According to recommended values from [14], parameters of the model in the considered case were specified as w cr = 37 mm and α I = 0.88. Finally, an average value of w I was determined as where N is number of days of some time interval and (w I ) n is w I calculated from expression (14) for a day with number n. We considered a time interval of 30 years and obtained w av I = 1.489 mm/day. Taking into account evapotranspiration, w u can be determined from w av I as where c et is a dimensionless coefficient. Climatic conditions of the Voronezh region correspond to c et = 0.3, and finally w u = 4.46 × 10 −4 m/day. w l characterizes groundwater inflow from the lower aquifer beneath the argillaceous layer A (Figure 4b). This value is equal to vertical component of Darcy's velocity in the lowpermeability layer between these aquifers. Hydraulic conductivity of upper Cretaceous clays at the considered site is approximately equal to 10 −4 m/day, according to data from the enterprise "Hydrospetsgeologiya". One can expect that argillaceous clay of the lowpermeability layer is even less. We denote this value as f arg . The difference between hydraulic head in the lower aquifer and level of the groundwater in the upper aquifer is about 12 m. Thickness of the waterproof layer at the site is about 20 m (Figure 4b). Hence, we can estimate w l as Since w l /w u << 1, w l can be neglected in (13). Subsequent procedure of calibration is absolutely the same, as in the case of the confined aquifer. Obtained distribution of water conductivity is shown in Figure 9. One can see that spatial heterogeneity of hydraulic conductivity is not significant to what corresponds to results of exploratory drilling of the aquifer.
Subsequent procedure of calibration is absolutely the same, as in the case of the confined aquifer. Obtained distribution of water conductivity is shown in Figure 9. One can see that spatial heterogeneity of hydraulic conductivity is not significant to what corresponds to results of exploratory drilling of the aquifer.

Conclusions
A characteristic feature of aquifers is elevated permeability in comparison with enclosing rocks. In a case of a pollution by liquid radioactive waste, aquifers represent a substantial hazard in view of elevated velocity of groundwater flow. Pollution can take place in both confined and free aquifers. An important part of input data for calculation of groundwater velocity in the polluted aquifers is spatial distribution of their permeability, which can be substantially heterogeneous. This can take place even at low permeability gradients due to significant dimensions of the domain. The only data for determination of such distribution consists of results of measurements in exploratory boreholes. However, drilling of the latter is a time consuming and expensive procedure. That is why the number of such measurements is usually limited. In addition to that, measurements are of different types (transmissibility and hydraulic head in confined aquifers and hydraulic conductivity and elevation of hydraulic head in free aquifers). We proposed a technique for determination of spatial distribution of transmissibility or hydraulic conductivity on the basis of measurements of all the types, simultaneously. The technique reduces the problem of the distribution determination to a search of a conditional minimum of a functional, with conditions in the form of equations in partial derivatives. The successful solution of the problem is demonstrated on examples of a confined aquifer (site of deepinjection disposal of liquid radioactive waste in the Krasnoyarsk region) and a free aquifer (site of radioactive pollution at the territory of atomic power plant Novovoronezhskaya) in Russia.

Conclusions
A characteristic feature of aquifers is elevated permeability in comparison with enclosing rocks. In a case of a pollution by liquid radioactive waste, aquifers represent a substantial hazard in view of elevated velocity of groundwater flow. Pollution can take place in both confined and free aquifers. An important part of input data for calculation of groundwater velocity in the polluted aquifers is spatial distribution of their permeability, which can be substantially heterogeneous. This can take place even at low permeability gradients due to significant dimensions of the domain. The only data for determination of such distribution consists of results of measurements in exploratory boreholes. However, drilling of the latter is a time consuming and expensive procedure. That is why the number of such measurements is usually limited. In addition to that, measurements are of different types (transmissibility and hydraulic head in confined aquifers and hydraulic conductivity and elevation of hydraulic head in free aquifers). We proposed a technique for determination of spatial distribution of transmissibility or hydraulic conductivity on the basis of measurements of all the types, simultaneously. The technique reduces the problem of the distribution determination to a search of a conditional minimum of a functional, with conditions in the form of equations in partial derivatives. The successful solution of the problem is demonstrated on examples of a confined aquifer (site of deep-injection disposal of liquid radioactive waste in the Krasnoyarsk region) and a free aquifer (site of radioactive pollution at the territory of atomic power plant Novovoronezhskaya) in Russia.