Modelling of the Fate of ¹³⁷Cs and ⁹⁰Sr in the Chornobyl Nuclear Power Plant Cooling Pond before and after the Water Level Drawdown

Abstract

During the accident in April 1986, the Cooling Pond (CP) of the Chornobyl Nuclear Power Plant (ChNPP) was heavily contaminated by fuel particles and radionuclides of cesium-137 (¹³⁷Cs) and strontium-90 (⁹⁰Sr). Starting from the end of 2014, a gradual decrease of the CP water level began leading to the transformation of the whole reservoir into eight separate sectors and raising the concern of the fate of ¹³⁷Cs and ⁹⁰Sr in the future. In this study, two mathematical models were applied to reproduce radioactive contamination of the CP from 1986 to 2021 and to provide a forecast of ¹³⁷Cs and ⁹⁰Sr concentrations in the CP water from 2022 to 2030. The hydrodynamic model THREETOX provided three-dimensional (3D) currents in the CP corresponding to hydrological conditions before and after water level drawdown, and these currents were used in the box model POSEIDON-F for the long-term simulations of the changes in ¹³⁷Cs and ⁹⁰Sr concentrations in water, bottom sediments, and biota. Seasonal changes in the distribution coefficient (K_d) describing the partition of ¹³⁷Cs between water and sediments were considered in the box model, which allowed us to reproduce the observed variations of concentration. Calculated concentrations of ¹³⁷Cs and ⁹⁰Sr in water and freshwater fish occupying different trophic levels agreed well with measurements for the entire post-accident period. After the water level drawdown, concentrations of ¹³⁷Cs in the CP water slightly increased in all eight sectors, while ⁹⁰Sr concentrations significantly increased in sectors close to ChNPP, which was explained by an additional ⁹⁰Sr source when comparing the simulation results and measurement data. Using the model forecast from 2022 to 2030, we predict that the concentration of both radionuclides will gradually decrease in new water bodies of the Cooling Pond except in the northern sectors, where the suggested additional source of ⁹⁰Sr will lead to a stabilization of ⁹⁰Sr concentrations.

1. Introduction

During the accident in April 1986, the Cooling Pond (CP) of the Chornobyl Nuclear Power Plant (ChNPP) was heavily contaminated by fuel particles and numerous radionuclides, among which cesium-137 (¹³⁷Cs) and strontium-90 (⁹⁰Sr) were a significant part, making it one of the most radioactively contaminated large bodies of water on the globe. The CP is an artificial reservoir that was initially built in 1976 and completed in 1981 to meet the technological needs of the operating nuclear power plant and is located on the right-bank floodplain of the Prypiat River. The banks of the CP are partly formed by a floodplain terrace and mainly by a protective dam, which is 25 km long, 70–100 m wide, and 5.7 m high. The area of the CP was 22.7 km², with the water depth mainly between 4 and 7 m, except in some areas of up to 18–20 m depth. The CP volume was 149 million m³ at the normal maximum level, and the water level of the CP was 6–7 m higher than the Prypiat River water levels. Since the CP water seeped intensively through the dam and flowed to the underlying aquifer, the continuous pumping of water from the Prypiat River was required to cover these water losses. During the period of the ChNPP operation, the heated water was release from the cooling system of the NPP, preventing ice formation in the CP. After the shutdown of the last working reactor unit in 2000, the entire surface of the CP was covered by ice under sub-zero air temperatures.

The ichthyofauna of the CP was formed due to (1) natural populations of aboriginal species of fish that inhabit the Prypiat River, its backwaters, and floodplain lakes where the CP was built; (2) targeted introduction of other species with subsequent acclimatization; (3) periodic replenishment of fish population of introduced species, which were not capable of reproduction in the conditions of the CP, by stocking with juveniles artificially obtained at fish farms; and (4) replenishment of the diversity of ichthyofauna by different fish species from the Prypiat River and nearby bodies of water. During the CP operation time, a highly productive semi-natural ecological system was established with high biodiversity of hydrobionts of various trophic levels and ecological groups including higher aquatic plants, phyto- and zooplankton, phyto- and zooperiphyton, zoobenthos, and fish species.

The ChNPP accident contaminated the ellipse-shaped CP with a length of about 11.5 km and width of about 2.2 km from the direct atmospheric fallout of aerosols and fuel particles on the water’s surface, the release of highly contaminated water from a system of accidental cooling, and water used to extinguish the fire [1]. In the years after the accident, the contamination was redistributed in the CP under the impact of three types of operational and hydrological conditions mentioned herein. For the 1986–2000 period, the currents were forced by the cooling system of the ChNPP, causing a circular water movement in the CP. Following the shutdown of all reactor units of the ChNPP, natural circulation took place in the CP driven by wind and seasonal temperature stratification. During this entire time, continuous pumping fed the CP with the Prypiat River to compensate for the CP water losses leading to a slow decrease in the content of radionuclides. From the end of 2014, the pumping system was stopped, and the CP drawdown began with a gradual decrease in the CP water level. Currently, the CP water level has dropped by about 6.5 m, transforming the CP reservoir into several small bodies of water.

Due to the drawdown of the water level, large shallow areas of the CP were drained, leading to the die-off of a significant part of the aquatic organisms and to signs of eutrophication in the small bodies of water that formed on the site of the CP, as evidenced by the disappearance of bivalve mollusks. Among the fish population, the number of species sensitive to water quality decreased, and the number of species more adapted to adverse conditions of the water environment increased.

During the 35 years following the accident, radioactive contamination of the CP was investigated by numerous field and model studies. A comprehensive review of radioactivity measurements conducted in the water and bottom sediments of the CP was performed by [2]. Measurements of radionuclide concentrations in aquatic biota were reported by [3,4]. Three-dimensional (3D) modelling of the transport of ¹³⁷Cs in the water and bottom sediments was performed for one year following the accident [5]. Radioactive contamination of aquatic organisms living in the CP was modelled by [6,7,8] for several years following the accident. However, no studies have focused on the long-term modelling of the radionuclide concentrations in the main components of the CP aquatic environment, which is required to predict changes in these concentrations after the water level drawdown. Long-term modelling requires validation using measurement data for the whole post-accident period before the drawdown and new data collected during the drawdown period by the joint Ukrainian and Japanese investigations [9].

The objectives of this study were to reproduce and forecast changes in ¹³⁷Cs and ⁹⁰Sr concentrations in water, bottom sediments, and aquatic organisms occupying different trophic levels before and after the CP drawdown. We used two numerical models that were customized to the CP freshwater environment with old and new measurement data to reproduce the radioactive contamination of the CP from 1986 to 2021. From these modelling results, we provided a forecast of ¹³⁷Cs and ⁹⁰Sr concentrations in the newly formed water bodies on the site of CP from 2022 to 2030.

2. Materials and Methods

Two connected models were used for the simulations: the 3D model THREETOX and the box model POSEIDON-F, which are described below. The THREETOX model simulated 3D currents in the CP for each type of hydrological conditions, and the simulated currents were used to calculate fluxes of water between boxes in the box model POSEIDON-F. The POSEIDON-F model, which is developed in this study as a freshwater version of the POSEIDON-R model, was applied for the long-term simulations of the changes in activity concentration in the water, bottom sediments, and biota starting from 1986. Different types of input data were required for THREETOX and POSEIDON-F models. They are summarized in

Table 1. List of input data and their sources used in modelling.

Type of Input Data	Model	Source
Meteorological data	THREETOX	ERA-Interim reanalysis data
ChNPP cooling system release rate	THREETOX	Average value for such type of reactor
Water heating rate in the ChNPP cooling system	THREETOX	Average value for such type of reactor
Bathymetry data	THREETOX, POSEIDON-F	Digital map reported in [2]
Water fluxes between boxes	POSEIDON-F	Calculated from the THREETOX model output
Water losses and their replenishment from the Prypiat River	POSEIDON-F	Estimations from [1]
Deposition densities of radionuclides and their total activity	POSEIDON-F	Estimations from [2]
Seasonal changes in ammonium and potassium ions concentration	POSEIDON-F	Measurement data from authors’ own databases
Concentration Factors for phytoplankton and macro-algae	POSEIDON-F	Calculated as a ratio of radionuclide concentrations in corresponding organism and water from authors’ own databases
Biological half-lives of radionuclides in aquatic organisms	POSEIDON-F	Data for freshwater biota adopted from [10]

and described in detail later in the text.

The Chornobyl accident led to the atmospheric deposition of radionuclides on large territories, with the maximal density around the ChNPP, where the Cooling Pond is located (Figure 1a). A map with estimated deposition densities of ⁹⁰Sr close to the ChNPP is shown in Figure 1b. The main processes in the first months after the accident were the decay of short-lived radionuclides and the adsorption of other radionuclides by sediments that were subsequently deposited to the CP bottom. Due to the large size of the CP, radionuclides deposited on its surface non-uniformly with different deposition densities, which correlated with deposition on the surrounding land [5]. The total amount of radionuclides deposited on the CP in April–May 1986 was estimated as 180 TBq of ¹³⁷Cs and 35 TBq of ⁹⁰Sr [2]. In the POSEIDON-F model, atmospheric deposition was considered a dominant source term.

The transformation of the entire reservoir into several small bodies of water is illustrated by comparing Google Maps satellite images of the CP taken in 2014 (Figure 1a) and in 2020 (Figure 2a). For the model set-up, the CP was separated into eight sectors (Figure 2b) to describe small bodies of water that remained after the water level drawdown. The contamination of bottom sediments was analyzed by [2], taking into account multiple measurements over the whole CP in the 1999–2013 period. In this study, we integrated deposited inventories of ¹³⁷Cs in 2013 taken from the map [2] over eight sectors shown in Figure 2b and compared the total activities of ¹³⁷Cs in bottom sediments between measurement data and model simulation.

The long-term simulations were conducted for the following operational and hydrologic conditions: (1) the 1986–2000 period, when the ChNPP was in operation; (2) the 2001–2015 period, when the ChNPP operation was stopped, keeping the same water level in the CP before the drawdown; (3) the 2016–2021 period, when the CP water level decreased to the new natural state as of the end of 2021; and (4) the 2022–2030 period, which is used to forecast ¹³⁷Cs and ⁹⁰Sr concentrations. For the model validation, the simulation results were compared with measurement data for the whole post-accident period collected by the authors from their own databases and various publications and obtained during field trips in 2020. Statistical processing was used to evaluate the agreement of the model results with measured concentrations. To achieve this, we selected the geometric mean (GM) and the geometric standard deviations (GSD), which are preferable for data sets that have a range of values with several orders of magnitude and could be described by the log-normal distribution. To estimate GM and GSD values, we used simulated-to-observed ratios of ¹³⁷Cs and ⁹⁰Sr concentrations in water and different types of fish.

2.1. 3D Model THREETOX

The modelling system THREETOX was developed to simulate the dispersion of radionuclides and other contaminants on local and regional scales [5,11,12]. The system includes models of hydrodynamics, ice thermodynamics, and sediment and radionuclide transport. In the current study, only the hydrodynamics model was applied. The model describes the motion of an incompressible fluid using a system of Reynolds-averaged equations of continuity and horizontal momentum using the Boussinesq and hydrostatic approximations. The system also includes the heat transfer equation and the equation of state, describing the dependence of the density of water on its temperature and pressure. The system of equations completes with equations of the k-ε turbulence model, which describes the coefficients of vertical turbulent viscosity and diffusion. Prognostic variables of the hydrodynamic model are the three components of the velocity vector field, temperature, and surface level elevation. A numerical algorithm was implemented in the horizontal curvilinear-orthogonal coordinate system.

The THREETOX model was previously used for modelling the radioactivity contamination of the Chornobyl CP after the preliminary validation of the hydrodynamic module by measurement data of the temperature distribution in the surface layer of the CP before the accident [5]. The modified hydrodynamic module has been validated for the Cooling Pond of the South Ukrainian NPP [13]. Since currents are closely related to temperature distribution in the body of water, the correct description of heat fluxes between water and the atmosphere is important for model accuracy. In the case of the operating NPP, excessive evaporation forms the special microclimate above the CP. It was shown in [13] that the formula describing heat fluxes between the atmosphere and water for natural conditions [11] does not work correctly for the CP when the NPP is in operation. Application of the formula described in [14] allows us to obtain better agreement between calculated and measured water temperatures in the pond in comparison with other parametrizations. This approach was used in the current study for the time period when the ChNPP was in operation. Parametrization for natural conditions was used when the NPP was out of operation.

The model was customized for the CP as follows. A new curvilinear orthogonal grid was created for the CP with horizontal resolution from 40 m to 100 m (Figure S1a in the Supplementary Materials). In the vertical direction, 15 layers were used in the σ–coordinate system. The depth in every node of the grid (Figure S1b) was calculated by the linear interpolation of bathymetry data reported in [2]. During the 1986–2000 period when the ChNPP was in operation, the cooling system pumped water, causing the circular movement of water in the CP directed from the outlet to the inlet of the cooling system through the entire CP. In the model, the water release rate was set up at 150 m³/s when three units were in operation and 100 m³/s when two units were in operation. The heating of water between the inlet and outlet was 10 degrees. The atmospheric forcing used ERA-Interim reanalysis data [15], including air temperature, wind speed and direction, relative humidity, cloudiness, and air pressure. Note that when the ChNPP was in operation, stable circulation existed in the CP with currents up to 0.2 m/s. After the ChNPP was decommissioned, the circulation in the CP changed according to changes in wind speed and direction (see examples of obtained surface currents in Figure S2). After the water level drawdown, the CP was transformed into several small bodies of water (Figure 2a) with a very limited exchange of water fluxes between some of them.

2.2. Box Model POSEIDON-F

The box model POSEIDON-PC was initially developed for the simulation of routine discharges from nuclear installations located on the coast [16]. To describe the accumulation of radionuclides in marine organisms as the result of accidental releases, the BURN extinction of the POSEIDON model was developed [17]. In the model, the dynamic equations for organisms occupying different trophic levels are solved, allowing the model to reproduce the time delay between the uptake of activity by biota and variations of radionuclide concentration in the seawater. In this model, the “target-tissue” approach was also introduced based on the assumption that each radionuclide is accumulating in a single specific organ. Finally, the benthic food web was added to the model [18] to describe the transfer of activity from contaminated bottom sediments to marine organisms through the food chain. Currently, there is a POSEIDON-R version of the model, which is integrated into the RODOS decision support system [19,20]. In the current study, a freshwater version of the model named POSEIDON-F was developed and applied to the CP of the ChNPP.

In all modifications of the POSEIDON model, the water environment is considered as a system of boxes of different sizes, which can be subdivided into several vertical layers in the water column (Figure S3a). Each box contains a constant concentration of suspended sediments, which are continuously settling down to the lower water layers and then to the bottom. The model assumes that the fraction of radionuclide dissolved in the water (liquid phase) is in instantaneous equilibrium with that adsorbed on suspended or bottom sediments fraction (solid phase), with the constant ratio between their concentrations equal to the distribution coefficient K_d. The temporal variations of radionuclide concentration in each compartment are calculated taking into account the exchange of activity with adjacent compartments, suspended and bottom sediment, as well as radioactive sources and decay. The transfer of activity from one compartment to another depends on the concentration of radionuclide in the initial compartment and is governed by fluxes of water due to advection and diffusion processes on the faces between boxes. A simple three-layer model is considered for describing the migration of radioactivity in bottom sediments. The exchange of activity between the top bottom layer and near bottom water layer is described by diffusion and bioturbation processes. Only the diffusion process is considered between the top and middle sediment layers. In addition, the settling of sediments to the lower layers determines the continuous flux of activity in the bottom sediments directed downward. The equations for water column layers (Equation (1)), and upper (Equation (2)) and middle (Equation (3)) sediment layers read as follows:

$$\frac{\partial C_{0ik}}{\partial t}={\displaystyle \sum _m{\displaystyle \sum _j\left[\frac{F_{jimk}}{V_{0ik}}{C}_{0jm}-\frac{F_{ijkm}}{V_{0jm}}{C}_{0ik}\right]}}+{\gamma }_{0ik}{C}_{0(i+1)k}-({\gamma }_{1ik}+\lambda )C_{0ik}+\frac{L_{t,k}}{h_{1k}}{\gamma }_{2k}{C}_{1k}+Q_{sik},$$

(1)

$$\frac{\partial C_{1k}}{\partial t}=-({\gamma }_{2k}+{\gamma }_{3k}+\lambda )C_{1k}+\frac{h_{1k}}{L_{t,}{}_k}{\gamma }_{1k}{C}_{01k}+\frac{L_{m,k}}{L_{t,k}}{\gamma }_{4k}{C}_{2k},$$

(2)

$$\frac{\partial C_{2k}}{\partial t}=-({\gamma }_{4k}+{\gamma }_{5k}+\lambda )C_{2k}+\frac{L_{t,k}}{L_{m,k}}{\gamma }_{3k}{C}_{1k}.$$

(3)

where C_0ik is the spatially averaged concentration of radionuclide in the water column layer i of box k; i = 1 corresponds to the near bottom water column layer; C_1k and C_2k are the averaged concentration of radionuclide in the upper and middle sediment layers of box k, respectively; λ is the radionuclide decay constant; F_ijkm is the water flux from layer i of box k to layer j of box m; V_0ik is the volume of layer i of box k; h_ik is the thickness of the water column layer i of box k; L_t,k and L_m,k are the thicknesses of the top and middle bottom sediment layers of box k, respectively; Q_sik is the source of the activity in layer i of box k; γ_0ik…γ_5k are the radionuclide transfer rates in the system water–suspended sediments–bottom sediments, t is the time.

The principle of the biological uptake model is based on grouping aquatic organisms into a limited number of classes by taking into account the trophic level and types of species (Figure S3b). Two pathways for radioactivity uptake are considered in the model: directly from the water and through the food chain. The model includes pelagic and benthic food chains. Pelagic organisms are grouped into phytoplankton, zooplankton, non-predatory (prey) fish, and predatory fish. The benthic food chain includes the organic part of the bottom deposit, macro-algae, deposit-feeding invertebrates, demersal fish, and benthic predators. Crustaceans, mollusks, and coastal predators feed on both pelagic and benthic organisms in shallow coastal areas.

Due to the rapid uptake from water and the short retention time of radioactivity, the concentration of radionuclides in phytoplankton is calculated using the Biological Concentration Factor (BCF) approach by multiplying the concentration of radionuclide in water by the corresponding BCF value. For the macroalgae, a dynamic model is used to describe radionuclide concentrations due to the longer retention times

$$\frac{d{C}_{ma}}{dt}=\left(C{F}_{ma}{C}_w-C_{ma}\right)\frac{\ln 2}{T_{0.5,ma}},$$

(4)

where C_w and C_ma are the radionuclide concentrations in the water and macroalgae, respectively, CF_ma is the corresponding BCF, and T_0.5,ma is the biological half-life of the radionuclide in the macroalgae. The concentration of a given radionuclide in other considered aquatic organisms is described by the following differential equation:

$$\frac{d{C}_i}{dt}=a_i{K}_{f,i}{C}_{f,i}+b_i{K}_{w,i}{C}_w-\frac{\ln 2}{T_{0.5,i}}{C}_i,$$

(5)

where C_i and C_f,i are the radionuclide concentrations in the i-th organisms and their food, respectively, a_i is the food extraction coefficient (assimilation rate), b_i is the water extraction coefficient, K_f,i is the food uptake rate, K_w,i is the water uptake rate, and T_0.5,i is the biological half-life of the radionuclide in the organism.

The activity concentration in the food of a predator C_f is expressed by the following equation, summing for a total of n prey types

$${{\displaystyle C}}_f={\displaystyle \sum _{i=0}^n{{\displaystyle C}}_{prey,i}{{\displaystyle P}}_{prey,i}\frac{{{\displaystyle drw}}_{pred}}{{{\displaystyle drw}}_{prey,i}}},$$

(6)

where C_prey,i is the activity concentration in prey of type i, P_prey,i is the preference for prey of type i, drw_pred is the dry weight fraction of predator, and drw_prey is the dry weight fraction of prey of type i. The index “0” corresponds to the organic deposit in bottom sediments. Values of the model parameters are discussed in [21]. The generic parameters of the model were calibrated for different aquatic environments [18,21]. The sensitivity of parameters was investigated in [18].

Note that, in reality, the population of each type of aquatic organism consists of individuals of different age classes. Each age class usually has different metabolic rates (i.e., food and water uptake rates, biological half-life) and even feeding structures, which must be described by different model parameters. For example, juvenile organisms have higher metabolic processes and were contaminated more intensively compared to adult organisms. This was observed in the CP during the first year after the accident, when the concentration of ¹³⁷Cs in older fish was lower than that in young fish of the same species [6]. The POSEIDON-F model considers only adult individuals of aquatic organisms. This means that the model does not distinguish the variability of organisms in terms of their size, age, and lifespan. To account for this, further development and testing of the model is required for the CP or other contaminated body of water where sufficient measurement data are available.

According to [22], the distribution coefficient (K_d) for ¹³⁷Cs in the freshwater environment depends on concentrations of ammonium [NH₄] and potassium [K] ions. Changes in [NH₄] and [K] are due to eutrophication processes such as the growth and death of planktonic species during the year, which leads to changes in K_d values. These changes are important to include for ¹³⁷Cs because the most of ¹³⁷Cs radionuclide is deposited in bottom sediments, and even small changes in K_d values may lead to significant changes in activity concentration in the water. Therefore, the corresponding changes were introduced in the POSEIDON-F model to replace the constant K_d with the following formula:

$$K_d=\frac{8}{[K]+5[N{H}_4]},$$

(7)

where measured values of [NH₄] and [K] available from July 2002 to September 2003, and these monthly values (

Table 2. Concentrations of ammonium [NH₄] and potassium [K] ions measured in the deep part of Sector 7 (see Figure 2b) from 07/2002 to 09/2003.

Month	Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep	Oct	Nov	Dec
[K], mg L−1	1.16	2.19	1.49	1.20	2.25	3.12	2.86	1.61	2.15	1.29	1.34	1.42
[NH4], mg L−1	0.27	0.11	0.45	0.35	0.40	0.29	0.42	0.60	0.23	0.20	0.07	0.01

) were used in all sectors of the CP for the entire simulation period.

Standard parameters for the simulation of radionuclide uptake by aquatic organisms in the POSEIDON-F model are adopted from the marine to freshwater environment. According to [23], decrease of water salinity leads to an increase in the organism’s uptake rate for cesium and strontium due to a decrease in the concentration of competing potassium and calcium ions, respectively. For ¹³⁷Cs radionuclide, it is described by changing BCF for phytoplankton (CF_ph) and macro-algae (CF_ma), which are at the beginning of the food web and ensure the radionuclide transfer through the whole food web. In the present study, these CFs were defined from average relations between measured concentrations of ¹³⁷Cs in the organism and water. For the CP, they are CF_ph = 150 L/kg and CF_ma = 500 L/kg. In addition, the biological half-lives for the freshwater environment were adopted from [10] for all considered types of organisms.

To represent the described processes, a new box system was created for the CP, dividing the reservoir into eight sectors (Figure 2b), where each sector includes shallow (outer) and deep (inner) boxes. The outer boxes have one vertical layer in the water column, whereas inner boxes have two vertical layers. After the water level drawdown, only inner boxes remain in the model (Figure 2), therefore their shape approximately repeats the shape of new bodies of water, which were formed at each CP sector. The volumes and depths of boxes were calculated from bathymetry data of the CP, making a structure that is able to reproduce the circulation in the reservoir before and after water level drawdown. Fluxes of water between boxes were calculated by the averaging of 3D currents obtained using the THREETOX model. They are different for three hydrological states that existed in the CP from the Chornobyl accident to the present. For long-term modelling, the filtration of water through the dam around the CP and the pumping of corresponding volumes of water from the Prypiat River is an important factor for the balance of radionuclides in the CP. According to [1], around ⅔ of the whole volume of the CP was filtered through the dam for every year that was considered in the modelling.

3. Results and Discussion

3.1. Concentrations of ¹³⁷Cs in the CP Environment for the Period 1986–2021

Figure 3 shows simulated and measured ¹³⁷Cs concentrations in the CP water from 1986 to 2021 for Sectors 1 (Figure 3a), 4 (Figure 3c), and 8 (Figure 3d). Figure 3b is a zoom-in of Figure 3a from 2000 to 2008 to demonstrate seasonal variations of ¹³⁷Cs concentration. According to the results of the simulations, the highest concentration of ¹³⁷Cs in the CP water varies from 1.5 × 10⁶ Bq/m³ in the northern part (Figure 3a,d) to 1.0 × 10⁶ Bq/m³ in the southern part (Figure 3c) in 1986. In the first years after the accident, most of ¹³⁷Cs was deposited at the CP bottom, which resulted in reduction of ¹³⁷Cs concentrations in the water column by approximately 50 times to 2.0 × 10⁴ Bq/m³. After that, the seasonal variations of ¹³⁷Cs concentration in the CP water appeared due to periodical radionuclide adsorption–desorption processes in the system of “water–bottom sediments” that followed changes in K_d. According to Equation (7), the lowest K_d value is expected in the middle of summer that corresponds to high concentrations of ammonium and potassium ions. Low K_d means that the amount of ¹³⁷Cs activity in bottom sediments decreases due to release of activity to water. As a result, the highest concentration of ¹³⁷Cs in CP water occurs in early autumn each year, and this interesting phenomenon requires detailed investigations in future studies.

Measurement data show that the concentration of ¹³⁷Cs in water can vary from 2.5 to 6 times during the year, and this variability agrees with model results (Figure 3). The model results indicate that the obtained amplitude of variations is less than measured in some years, and rare measurement data of ¹³⁷Cs concentrations are about an order of magnitude smaller between sectors. The nature of such variability could be explained by changes in the concentrations of potassium and ammonium ions, which have similar biochemical properties as ¹³⁷Cs radionuclide, and this may indicate local eutrophication with extreme variations of ammonium and potassium ions affecting K_d. Since most radionuclide activity accumulates in sediments, even small changes in K_d lead to large changes in radionuclide concentration in water. In addition, the unknown source of the measured outliers can exist. The identification of the outlier’s source can only be carried out after the resumption of the CP monitoring.

In Figure 3, the modelling results show that seasonal fluctuations in ¹³⁷Cs concentrations in the water increased after the drawdown of the water level in the CP. It is important to note that this increase of ¹³⁷Cs concentrations is obtained for a situation where the changes in potassium and ammonium ion concentrations remain the same as before the CP drawdown. In general, the agreement between simulated and observed concentrations of ¹³⁷Cs in water is confirmed by values of GM = 1.19 and GSD = 1.64 for N = 1296 records. The GM value of 1.19 states that the average ratio of simulated-to-observed concentration C_calc/C_obs = 1.19. The GSD value of 1.64 implies that most of the calculated concentrations (about 68%) are in the range between C_obs × 1.64 and C_obs/1.64. The comparison of all data measurements with corresponding model results (Figure S5) shows that the vast majority of the calculated concentrations of ¹³⁷Cs in water do not differ by more than two times from the measurements.

Using inventories of ¹³⁷Cs in the bottom deposition from the map [2] integrated over the eight sectors (Figure 2b), we compare the total activities of ¹³⁷Cs in bottom sediments for 2013 with corresponding activities calculated by the model (

Table 3. Calculated and measured total activities of ¹³⁷Cs in bottom sediments in eight sectors of the CP in 2013.

Sector	1	2	3	4	5	6	7	8
Activity based on model results, TBq	8.5	6.3	22.8	23.4	25.6	22.6	18.2	23.2
Activity based on measurement data, TBq	11.0	25.1	58.8	25.3	38.0	19.6	25.4	20.0

). It is important to note that a relatively small number of measurement points in such a large reservoir can lead to inaccuracies in the calculation of the radionuclide inventories in bottom sediments based on measurement data. Overall, modelling results slightly underestimate the total activities with an average relation of calculated-to-measured data around 0.76, which could be due to the simplified description of radionuclide distribution in the water-sediment system used in the box model. Since the processes of suspended sediment transfer and bottom sediment movement were not considered in the model, these could play an important role in the redistribution of contaminated sediments at the bottom and in their deposition in the deepest parts of CP. This could explain the fact that the largest disagreement between the calculated and measured activities is observed in Sectors 1–3, which are located in the area of the highest water flow from the ChNPP outlet with the greatest disturbance of bottom sediments during the ChNPP operation. In addition, the unknown direct releases of radioactively contaminated water from the ChNPP site to the CP took place at the time of the accident and can have a significant effect on the local contamination of bottom sediments close to the NPP outlet. These unknown direct releases are not considered in present modelling study due to their minor activity in relation to atmospheric deposition and the impossibility of quantification, though these could be the aim of future modelling investigations.

The comparison between calculated by the POSEIDON-F model and measured concentrations of ¹³⁷Cs in fish is shown for Sectors 1 and 7, which have the largest amount of measurement data (Figure 4, Figure 5 and Figure 6). Since the POSEIDON-F model gives the concentration of radionuclides in several types of fish occupying different trophic levels, measurements in all available fish species were grouped according to the POSEIDON-F classification. In Figure 4, calculated concentrations of ¹³⁷Cs in non-predatory (prey) fish are compared with corresponding measurements in silver carp and crucian carp. From the model results, the highest concentration of ¹³⁷Cs in non-predatory fish is simulated at the end of 1986 with a rapid decline from 1987 to 1990 and a half-life of 0.6 years. During the 1990–2000 period, the elimination rate of ¹³⁷Cs decreased leading to a half-life of 3.5 years. Starting from 2000, the elimination became very slow, with a slight increase in ¹³⁷Cs concentrations after the water level drawdown.

Similar changes in calculated concentrations are shown for coastal predators feeding on both pelagic and benthic organisms, which was compared with perch and pikeperch (Figure 5), and for demersal fish, which was compared with bream (Figure 6). The difference corresponds to more smooth changes in ¹³⁷Cs concentrations for fish species from higher trophic levels. In addition, the highest concentration of ¹³⁷Cs radionuclide in organisms with a benthic diet in first years after the accident is lower compared to pelagic organisms due to the time delay between peaks of the water and sediment concentrations. After the equilibration, the concentration of ¹³⁷Cs in organisms with a benthic diet become higher due to an additional pathway of activity from bottom sediments to fish through the food web [18,24,25]. Figure 4 and Figure 5 show that predatory fish have higher calculated and observed concentrations compared with non-predatory fish, indicating the importance of trophic levels. This finding is consistent with [26], which investigated the effect of mass and trophic level on the concentration of radionuclides in fish, and with [6], where the ¹³⁷Cs concentrations in non-predatory and predatory fish in the CP were simulated for the first five years following the accident. A similar conclusion was made by [27] based on a series of measurements indicating that aquatic organisms of greater mass, as well as higher trophic levels, have a higher concentration of ¹³⁷Cs radionuclide. In summary, model results of ¹³⁷Cs concentrations agree well with measurements for non-predatory fish (GM = 0.8, GSD = 1.68, N = 60) and demersal fish (GM = 0.94, GSD = 1.26, N = 18) and underestimate measurements for coastal predators (GM = 0.68, GSD = 1.87, N = 62).

Since fish remain in the contaminated environment, the concentration of ¹³⁷Cs radionuclide in all fish species followed the trends of the corresponding concentration in water because radionuclides are assimilated from water and then transferred through the food chain. The differences in these processes for different species depend on the trophic level and feeding ration. Therefore, after the initial maximum concentration, there was a rapid decline, which changed to a gradual decline caused mainly by radioactive decay. The drawdown of the water level in the CP starting from 2014 has led to some increase in the radionuclide concentration in fish, after which the gradual decline continued. According to model results, the seasonal variations of ¹³⁷Cs concentrations in fish follow the corresponding variations of ¹³⁷Cs concentrations in water. Due to the delay in the uptake of activity in fish, we observe an increase in radionuclide concentrations in the late autumn–early winter and decrease in the amplitude of seasonal variations. Predatory fish have a longer delay and smaller variations in comparison with prey fish. The simulated seasonal variations of the ¹³⁷Cs concentrations in fish are quite low, and it is lower than the accuracy of the measurement methods. A detailed study with very precise measurements is needed in order to obtain corresponding experimental data.

3.2. Concentrations of ⁹⁰Sr in the CP Environment for the Period 1986–2021

Figure 7 shows the results of simulated and measured ⁹⁰Sr concentrations in the CP water from 1986 to 2021. The concentration decreases approximately 100 times during the first 5 years after the accident mostly due to the deposition of ⁹⁰Sr in bottom sediments. After that, the decrease of ⁹⁰Sr concentration becomes much slower and occurs only due to the radioactive decay and filtration of ⁹⁰Sr radionuclide through the dam around the CP. After the drawdown of the water level, the concentration of ⁹⁰Sr in the CP has increased similarly to the concentration of ¹³⁷Cs, and measurements show that this increase is not uniform in the CP. For example, the increase in northern sectors was observed to be much larger than in southern sectors.

One of the reasons for such a non-uniform increase is the presence of an additional source of ⁹⁰Sr radionuclide in the areas close to the ChNPP. The origin of such a source can be different and requires further research studies for its identification. In particular, it could be associated with the emission of ⁹⁰Sr due to the destruction of fuel particles on the drained areas of the CP. Results of experiments described in [2] show that ¹³⁷Cs radionuclide in the CP exists mostly in non-exchangeable forms, while a large part of ⁹⁰Sr radionuclide could be remobilized from sediments into the water after oxidation. The drainage of large areas of the CP has led to the situation when a lot of ⁹⁰Sr is on the terrain surface and in contact with the air. As a result of rain events, some quantity of oxidized ⁹⁰Sr remobilizes from sediments and moves downhill with the surface run-off to the new water bodies formed on the site of the CP. This process could be a reason for the additional increase of ⁹⁰Sr concentration after the water level drawdown.

Measurements show that the highest increase in ⁹⁰Sr concentration is observed in northern sectors of the CP (Figure 7a,d), which are close to the ChNPP, and had the maximum deposition density that took place in 1986. The increase of ⁹⁰Sr concentration in southern sectors is much smaller (Figure 7b,c). The investigation of drained areas shows that the highest concentrations of radionuclides were found in parts of the CP that belong to Sectors 1 and 8 (see Figure S4). These results support the hypothesis that the significant increase in ⁹⁰Sr concentration in some sectors of the CP could be explained by its remobilization from sediments in the drained areas, which were the most contaminated due to the deposition of fuel particles after the accident. Such remobilization became possible only starting in 2015, when the first drained areas appeared due to the drawdown of the water level in the CP. To quantify this source, the application of the model gives that the amount of washed ⁹⁰Sr radionuclide should be around 3 GBq per year to fit the simulation results to the measurement data as shown by the dashed line in Figure 7a,d. The increase of ⁹⁰Sr concentration in southern sectors (Figure 7b,c) occurs naturally due to the redistribution of activity between water and sediments after the water level drawdown and does not require the addition of any external source of ⁹⁰Sr.

A significant increase in ⁹⁰Sr concentration in monitoring wells in the northern part of the CP after the water level drawdown was shown by recent groundwater monitoring data [28,29]. One of the explanations for this phenomenon could be a change in the direction of groundwater flow in this area resulting from the lowering of the water level in the CP by more than 6 m. New directions of groundwater flow may generate new sources of groundwater contamination in this area due to filtration through previously undisturbed post-accident temporary radioactive waste storage facilities on the territory at the ChNPP technical site. This process could also contribute to the increase in the ⁹⁰Sr concentration in northern sectors of the CP. Quantitative assessment of the role of new surface sources of ⁹⁰Sr fluxes to the CP and the increase in ⁹⁰Sr fluxes to groundwater of the northern part of the CP can be performed within the framework of the forthcoming monitoring and modelling studies. Such studies were interrupted due to the military invasion of the Chornobyl Exclusion Zone (CEZ) by Russian Federation troops in February–March 2022 and cannot be continued after the liberation of the CEZ by the Ukrainian army because of the large number of mines near the CP as a result of the consequences of military activities at the nuclear site. Therefore, a detailed analysis of the interaction of the CP waters with groundwater can only be conducted after the end of the Russian Federation’s military invasion of Ukraine.

Calculated total activities of ⁹⁰Sr in bottom sediments in different sectors of the CP were compared (

Table 4. Calculated and measured total activities of ⁹⁰Sr in bottom sediments in eight sectors of the CP in 2013.

Sector	1	2	3	4	5	6	7	8
Activity based on model results, TBq	0.47	0.36	1.36	1.46	1.63	1.30	1.04	1.30
Activity based on measurement data, TBq	3.21	4.80	13.1	3.95	4.70	3.02	4.00	3.80

) with corresponding values obtained from the digital map reported in [2] based on measurement data in 2013. Similar to ¹³⁷Cs, the model underestimates the total activities having the worse average relation of calculated-to-measured data, which is equal to 0.26. This could be explained by the larger fraction of exchangeable forms of ⁹⁰Sr radionuclide, which contributes to its more active redistribution in the CP under the water currents’ action. The possible inaccuracies in the calculation of the radionuclide inventories in bottom sediments based on measurement data can also contribute to this underestimation.

Since there are much fewer data of measurements for ⁹⁰Sr in fish compared with the ¹³⁷Cs data, we compare only the calculated concentration in non-predatory (prey) fish with measurement data in silver carp, crucian carp, and roach (Figure 8). The agreement between the calculated and measured concentrations in fish confirms the supposed value of the additional source of radionuclide in Sector 1 (Figure 8-left) and the absence or much lower value of the source in Sector 7 (Figure 8-right). The agreement between simulated and observed concentrations of ⁹⁰Sr in water in all sectors, where data of measurements are available, is confirmed by values of GM = 0.81 and GSD = 1.53 for N = 414 records. The corresponding comparison is plotted in the Supplementary Materials (Figure S6). For non-predatory fish, GM = 1.1 and GSD = 1.55 for N = 48.

3.3. Near-Future Forecast of ¹³⁷Cs and ⁹⁰Sr Concentrations in the CP Water

As described in Section 3.1, some increase in ¹³⁷Cs concentration and its seasonal variations in separate water bodies on the site of the CP are expected after the water level drawdown (Figure 3). To estimate such an increase, the POSEIDON-R model was applied with model parameters before drawdown, and the near-future forecast of ¹³⁷Cs concentrations is demonstrated in the water of Sector 1 without and with K_d changes (Figure 9). According to model results of ¹³⁷Cs concentration in Sector 1 (Figure 9a), the increase in yearly averaged concentration and seasonal variations is estimated as 36% and 1.34 times, respectively. Similar values are obtained for other sectors of the CP (not shown). These results correspond to the situation when the distribution coefficient K_d has not changed due to the water level drawdown. However, after the CP drawdown, the concentration of [K] and [NH₄] ions could change as a result of forming new biochemical conditions in each separated body of water that leads to the spatial and temporal variability of K_d values. Demonstrating uncertainties of K_d in simulations, the 30% change of K_d values in both directions shows the higher and lower declining trends of ¹³⁷Cs concentrations of in the CP (Figure 9b). In addition, the amplitude of seasonal variations of K_d could also be changed after the water level drawdown. All these factors need further research investigations, collecting field data with spatial and temporal measurements to be used in future modelling studies.

For the near-future forecast of ⁹⁰Sr concentrations, the increase in the concentrations after the water level drawdown differs among the eight sectors of the CP (Figure 10). In Figure 10, the largest increase is observed in the northern sectors compared to the southern ones. In this study, this spatial difference is explained by the presence of an additional source of ⁹⁰Sr that is introduced by the remobilization of exchangeable forms of radionuclide and its washout from drained areas or other unexplored phenomena such as changes in the groundwater flow direction. From simulation results, the increase of ⁹⁰Sr concentration in southern areas of the CP (Figure 10a) occurs only due to the redistribution of activity between water and sediments after the water level drawdown (Sector 5 is an example of such area). The concentration of ⁹⁰Sr in water in Sector 5 in 2021 is about 2.3 times higher than in 2014 before the water level drawdown and gradually decreases similarly as it was decreasing before the drawdown.

In northern areas of the CP (at least Sector 1 and Sector 8), a much higher increase in ⁹⁰Sr concentration is observed, reaching up to 3.9 times in Sector 8 and can only be simulated by introducing an additional source of unknown origin equal to 3 GBq of ⁹⁰Sr per year (Figure 10b). In the modelling of the near-future forecast, it is considered as a “black box” source term without identification of its origin that gives a direction for future research. The application of the “black box” source term allows us to make the prediction for the dynamics of ⁹⁰Sr concentration in the new bodies of water in the northern part of the CP from 2022 to the following years. In Figure 10b, the results of modelling show that the constant source of ⁹⁰Sr leads to the stabilization of ⁹⁰Sr concentrations in northern areas of the CP without a significant change.

To evaluate the annual source of 3 GBq, it can be compared with the inventories of ⁹⁰Sr in the bottom sediments and water in Sectors 1 and 8 where this source was introduced. The comparison shows that 3 GBq is less than 0.1% of 3.2 and 3.8 TBq (see Table 4) deposited in Sectors 1 and 8, respectively, while about 80 GBq of ⁹⁰Sr decay every year in sediments in these sectors. The inventories of ⁹⁰Sr in water from the model results were about 8 GBq in Sector 1 and 18 GBq in Sector 8 before the water level drawdown. As a result, the annual source of 3 GBq is a rather small source of ⁹⁰Sr in the context of its inventory in the bottom sediments of the CP, while it is considerable in view of the ⁹⁰Sr in water.

4. Conclusions

A modelling system that included two coupled models, the 3D model THREETOX and the box model POSEIDON-F, was applied to simulate ¹³⁷Cs and ⁹⁰Sr concentration changes in water, bottom sediments, and biota within the Cooling Pond (CP) of the Chornobyl Nuclear Power Plant (ChNPP), starting from the accident in 1986 to 2021 and including the CP drawdown period. In this study, we developed the POSEIDON-F model, which is a freshwater version of the POSEIDON-R model, to simulate the fate of ¹³⁷Cs and ⁹⁰Sr radionuclides in the CP of the ChNPP. The model was validated on the available measurement data for the entire post-accidental period. For ¹³⁷Cs simulation, the POSEIDON-F model was modified to represent seasonal variations of K_d values by connecting K_d with the concentration of potassium and ammonium ions, which vary according to the life cycle of aquatic flora. In addition, freshwater parameters for the food web were considered in the model. The simulations were continued from 2022 to 2030 to predict future changes in the ¹³⁷Cs and ⁹⁰Sr concentrations.

The model results were compared with measurement data collected by the authors from their own databases and various publications, and with the measured data obtained in field trips in 2020. The simulated-to-observed ratios of calculated and measured concentrations had the geometric mean (GM) range between 0.8 and 1.19, with the geometric standard deviations (GSD) from 1.26 to 1.68 for water, non-predatory and demersal fish, while the calculated concentration of ¹³⁷Cs in predatory fish underestimated the measurements with a smaller GM of 0.68 and larger GSD of 1.87. This indicates that, in general, the calculated ¹³⁷Cs and ⁹⁰Sr concentrations in water and freshwater fish occupying different levels of the food chain are in good agreement with measurements for the entire modelling period. We also found that the model underestimated the total activity accumulated in the bottom sediments based on measurement data prior to the CP drawdown. This discrepancy could be due to two reasons: (i) suspended sediment transport and bottom sediment movement, which are not considered in the model, could play an important role in the redistribution of contaminated sediments at the bottom and their deposition in the deepest parts of the CP; and (ii) uncertain direct releases of radionuclides with the ChNPP technical waters, which were not considered in the modelling study, could have a local effect on the contamination of CP bottom sediments close to the NPP outlet.

Since the CP is one of the most radioactively contaminated large bodies of water in the world, it was important to conduct the first study of the new conditions formed in the CP after water level drawdown and corresponding concentration changes in ¹³⁷Cs and ⁹⁰Sr radionuclides. Using the POSEIDON-F model developed to simulate the fate of ¹³⁷Cs and ⁹⁰Sr radionuclides in the freshwater environment and validated with measurement data available in the CP for the entire post-accidental period, this study showed that the increase in concentrations of both radionuclides in water and aquatic organisms after the CP water level drawdown occurred naturally due to the redistribution of radionuclides between water and bottom sediments. Although bottom sediments were more contaminated in the deeper areas of the CP, one redistribution process cannot explain the significant increase of ⁹⁰Sr in the areas close to the NPP. As a result, we hypothesized the presence of a ⁹⁰Sr source of unknown origin in these areas. There are several possible causes of this source requiring further field and modelling studies. In particular, it could be associated with the emission of ⁹⁰Sr due to the destruction of fuel particles on the drained shallow areas of the CP and the changes in the direction of groundwater flow after the CP drawdown. The value of this source was estimated to be 3 GBq per year to match the simulation results with the measurement data. The model also predicted that in the near future, the concentration of both radionuclides in new bodies of water, which were formed at the site of the former CP, will gradually decrease, with the exception of the northern sectors, where the revealed source of ⁹⁰Sr can stabilize ⁹⁰Sr concentrations for the period until 2030.