Modelling of the Fate of 137 Cs and 90 Sr in the Chornobyl Nuclear Power Plant Cooling Pond before and after the Water Level Drawdown

: 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 ( 137 Cs) and strontium-90 ( 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 137 Cs and 90 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 137 Cs and 90 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 137 Cs and 90 Sr concentrations in water, bottom sediments, and biota. Seasonal changes in the distribution coefﬁcient ( K d ) describing the partition of 137 Cs between water and sediments were considered in the box model, which allowed us to reproduce the observed variations of concentration. Calculated concentrations of 137 Cs and 90 Sr in water and freshwater ﬁsh occupying different trophic levels agreed well with measurements for the entire post-accident period. After the water level drawdown, concentrations of 137 Cs in the CP water slightly increased in all eight sectors, while 90 Sr concentrations signiﬁcantly increased in sectors close to ChNPP, which was explained by an additional 90 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 90 Sr will lead to a stabilization of 90 Sr concentrations.


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 ( 137 Cs) and strontium-90 ( 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 The objectives of this study were to reproduce and forecast changes in 137 Cs and 90 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 137 Cs and 90 Sr concentrations in the newly formed water bodies on the site of CP from 2022 to 2030.

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 and described in detail later in the text. 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] 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 90 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 137 Cs and 35 TBq of 90 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 137 Cs in 2013 taken from the map [2] over eight sectors shown in Figure 2b and compared the total activities of 137 Cs in bottom sediments between measurement data and model simulation. 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 137 Cs and 35 TBq of 90 Sr [2]. In the POSEIDON-F model, atmospheric deposition was considered a dominant source term.
(a) (b) 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 137 Cs in 2013 taken from the map [2] over eight sectors shown in Figure 2b and compared the total activities of 137 Cs in bottom sediments between measurement data and model simulation.   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 137 Cs in 2013 taken from the map [2] over eight sectors shown in Figure 2b and compared the total activities of 137 Cs in bottom sediments between measurement data and model simulation.
(a) (b) Figure 2. Satellite image of the CP in 2020 after water level drawdown (a) and schematic representation of the box system describing the CP in eight sectors connected and disconnected before and after the CP drawdown, respectively (b).

Figure 2.
Satellite image of the CP in 2020 after water level drawdown (a) and schematic representation of the box system describing the CP in eight sectors connected and disconnected before and after the CP drawdown, respectively (b).
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 137 Cs and 90 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 137 Cs and 90 Sr concentrations in water and different types of fish.

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 3 /s when three units were in operation and 100 m 3 /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.

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: 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 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: 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 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 137 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 137 Cs in the freshwater environment depends on concentrations of ammonium [NH 4 ] and potassium [K] ions. Changes in [NH 4 ] 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 137 Cs because the most of 137 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: where measured values of [NH 4 ] and [K] available from July 2002 to September 2003, and these monthly values ( Table 2) 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 137 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 137 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 2 /3 of the whole volume of the CP was filtered through the dam for every year that was considered in the modelling. According to the results of the simulations, the highest concentration of 137 Cs in the CP water varies from 1.5 × 10 6 Bq/m 3 in the northern part (Figure 3a,d) to 1.0 × 10 6 Bq/m 3 in the southern part (Figure 3c) in 1986. In the first years after the accident, most of 137 Cs was deposited at the CP bottom, which resulted in reduction of 137 Cs concentrations in the water column by approximately 50 times to 2.0 × 10 4 Bq/m 3 . After that, the seasonal variations of 137 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 137 Cs activity in bottom sediments decreases due to release of activity to water. As a result, the highest concentration of 137 Cs in CP water occurs in early autumn each year, and this interesting phenomenon requires detailed investigations in future studies.

Results and Discussion
the CP was filtered through the dam for every year that was considered in the modelling.  (Figure 3c) in 1986. In the first years after the accident, most of 137 Cs was deposited at the CP bottom, which resulted in reduction of 137 Cs concentrations in the water column by approximately 50 times to 2.0 × 10 4 Bq/m 3 . After that, the seasonal variations of 137 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 Kd. According to Equation (7), the lowest Kd value is expected in the middle of summer that corresponds to high concentrations of ammonium and potassium ions. Low Kd means that the amount of 137 Cs activity in bottom sediments decreases due to release of activity to water. As a result, the highest concentration of 137 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 137 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 137 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 137 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.

Concentrations of 137 Cs in the CP Environment for the Period 1986-2021
In Figure 3, the modelling results show that seasonal fluctuations in 137 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 137 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 137 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 137 Cs in water do not differ by more than two times from the measurements.
Using inventories of 137 Cs in the bottom deposition from the map [2] integrated over the eight sectors (Figure 2b), we compare the total activities of 137 Cs in bottom sediments for 2013 with corresponding activities calculated by the model (Table 3). 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 137 Cs in fish is shown for Sectors 1 and 7, which have the largest amount of measurement data (Figures 4-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             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 137 Cs concentrations for fish species from higher trophic levels. In addition, the highest concentration of 137 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 137 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]. Figures 4 and 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 137 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 137 Cs radionuclide. In summary, model results of 137 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 137 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 137 Cs concentrations in fish follow the corresponding variations of 137 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 137 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. Figure 7 shows the results of simulated and measured 90 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 90 Sr in bottom sediments. After that, the decrease of 90 Sr concentration becomes much slower and occurs only due to the radioactive decay and filtration of 90 Sr radionuclide through the dam around the CP. After the drawdown of the water level, the concentration of 90 Sr in the CP has increased similarly to the concentration of 137 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.

Concentrations of 90 Sr in the CP Environment for the Period 1986-2021
One of the reasons for such a non-uniform increase is the presence of an additional source of 90 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 90 Sr due to the destruction of fuel particles on the drained areas of the CP. Results of experiments described in [2] show that 137 Cs radionuclide in the CP exists mostly in non-exchangeable forms, while a large part of 90 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 90 Sr is on the terrain surface and in contact with the air. As a result of rain events, some quantity of oxidized 90 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 90 Sr concentration after the water level drawdown.
Measurements show that the highest increase in 90 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 90 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 90 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 90 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 90 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 90 Sr. One of the reasons for such a non-uniform increase is the presence of an additional source of 90 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 90 Sr due to the destruction of fuel particles on the drained areas of the CP. Results of experiments described in [2] show that 137 Cs radionuclide in the CP exists mostly in non-exchangeable forms, while a large part of 90 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 90 Sr is on the terrain surface and in contact with the air. As a result of rain events, some quantity of oxidized 90 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 90 Sr concentration after the water level drawdown.
Measurements show that the highest increase in 90 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 90 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 90 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 A significant increase in 90 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 90 Sr concentration in northern sectors of the CP. Quantitative assessment of the role of new surface sources of 90 Sr fluxes to the CP and the increase in 90 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 90 Sr in bottom sediments in different sectors of the CP were compared (Table 4) with corresponding values obtained from the digital map reported in [2] based on measurement data in 2013. Similar to 137 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 90 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 90 Sr in fish compared with the 137 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

Near-Future Forecast of 137 Cs and 90 Sr Concentrations in the CP Water
As described in Section 3.1, some increase in 137 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 137 Cs concentrations is demonstrated in the water of Sector 1 without and with Kd changes (Figure 9). According to model results of 137 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 Kd has not changed due to the water level drawdown. However, after the CP drawdown, the concentration of

Near-Future Forecast of 137 Cs and 90 Sr Concentrations in the CP Water
As described in Section 3.1, some increase in 137 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 137 Cs concentrations is demonstrated in the water of Sector 1 without and with K d changes ( Figure 9). According to model results of 137 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 4 ] 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 137 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.
iations 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 137 Cs concentrations is demonstrated in the water of Sector 1 without and with Kd changes (Figure 9). According to model results of 137 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 Kd has not changed due to the water level drawdown. However, after the CP drawdown, the concentration of [K] and [NH4] 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 Kd values. Demonstrating uncertainties of Kd in simulations, the 30% change of Kd values in both directions shows the higher and lower declining trends of 137 Cs concentrations of in the CP (Figure 9b). In addition, the amplitude of seasonal variations of Kd 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 90 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 For the near-future forecast of 90 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 90 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 90 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 90 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 90 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 90 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 90 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 90 Sr leads to the stabilization of 90 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 90 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 90 Sr decay every year in sediments in these sectors. The inventories of 90 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 90 Sr in the context of its inventory in the bottom sediments of the CP, while it is considerable in view of the 90 Sr in water.
southern ones. In this study, this spatial difference is explained by the presence of an additional source of 90 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 90 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 90 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 90 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 90 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 90 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 90 Sr leads to the stabilization of 90 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 90 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 90 Sr decay every year in sediments in these sectors. The inventories of 90 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 90 Sr in the context of its inventory in the bottom sediments of the CP, while it is considerable in view of the 90 Sr in water.

Conclusions
A modelling system that included two coupled models, the 3D model THREETOX and the box model POSEIDON-F, was applied to simulate 137 Cs and 90 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 137 Cs and 90 Sr radionuclides in the CP of the ChNPP. The model was validated on the

Conclusions
A modelling system that included two coupled models, the 3D model THREETOX and the box model POSEIDON-F, was applied to simulate 137 Cs and 90 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 137 Cs and 90 Sr radionuclides in the CP of the ChNPP. The model was validated on the available measurement data for the entire post-accidental period. For 137 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 137 Cs and 90 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 137 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 137 Cs and 90 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 137 Cs and 90 Sr radionuclides. Using the POSEIDON-F model developed to simulate the fate of 137 Cs and 90 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 90 Sr in the areas close to the NPP. As a result, we hypothesized the presence of a 90 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 90 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 90 Sr can stabilize 90 Sr concentrations for the period until 2030.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/w15081504/s1, Figure S1: numerical grid for the THREETOX model and bathymetry of the CP before water level drawdown; Figure S2: examples of the surface currents simulated by the THREETOX model; Figure S3: structure of a box system in the POSEIDON-R model and scheme of radionuclide transfers from the water and bottom sediment boxes to aquatic organisms; Figure S4: measurement data of the 90 Sr concentration in the top sediment layer (0-5 cm) on the dried areas of the ChNPP Cooling Pond; Figure

Data Availability Statement:
The measurement data presented in this study can be requested from the authors of the study.