THREE-DIMENSIONAL HINDCAST OF NITROGEN AND PHOSPHORUS BIOGEOCHEMICAL DYNAMICS IN LAKE ONEGO ECOSYSTEM, 1985–2015. PART

Despite a wide-ranging research, there is almost no information regarding the major biogeochemical fluxes that could characterize the past and present state of the European Lake Onego ecosystem and be used for reliable prognostic estimates of its future. To enable such capacity, we adapted and implemented a three-dimensional coupled hydrodynamical biogeochemical model of the nutrient cycles in Lake Onego. The model was used to reconstruct three decades of Lake Onego ecosystem dynamics with daily resolution on a 2 × 2 km grid. A comparison with available information from Lake Onego and other large boreal lakes proves that this hindcast is plausible enough to be used as a form of reanalysis. This model will be used as a form of studies of Lake Onego eco-system, including long-term projections of ecosystem evolution under different scenarios of climate change and socio-economic development.


Introduction
Being situated within the Baltic Sea drainage basin, boreal Lake Onego is the second-largest lake in Europe. Lake Onego is strongly phosphorus limited and still largely oligotrophic due to minor anthropogenic influence [1][2][3]. Total phosphorus (TP) inputs to Lake Onego, expressed either as a yield from the catchment area (15 kg TP km -2 yr -1 ) or per unit of the lake's water volume (3 mg TP m -3 yr -1 ), are similar to those in boreal Lake Ladoga (14 kg TP km -2 yr -1 ; 4 mg TP m -3 yr -1 , respectively), the largest lake in Europe. For comparison, the net anthropogenic input of phosphorus generated at the Great American Lakes' watershed has decreased over time from 21 kg TP km -2 yr -1 in 1987 to 10 kg TP km -2 yr -1 in 2012 [4]. The future evolution of the Lake Onego ecosystem may be driven by both the existing registered warming in its catchment area [5][6][7] and the projected effects of climate change [8][9], as well as socio-economic development [10][11][12]. Current projections indicate that river discharge will increase in the northern Baltic Sea region, subsequently increasing waterborne nutrient inputs [13] and mobilizing phosphorus reserves accumulated in the drainage area [14]. The interaction of these regionally occurring changes with the anticipated augmentation of local industrial, agricultural, mining, and forestry activities [15], including aquaculture [16], could generate synergetic ecosystem effects, such as those already occurring in the nutrient-rich Lake Winnipeg [17]. Hence, we need the capacity and a tool to reliably describe the current state of Lake Onego and to make prognostic estimates of future scenarios.
The development of mathematical models of large lake ecosystems and their implementation for producing recommendations on environmental protection measures has been ongoing for several decades [18][19][20][21][22][23][24][25], including modelling of Ladoga and Onego lakes [1,[26][27][28][29][30][31][32]. The main attention in Lake Ladoga and Onego modelling has been focused on spatial-temporal changes in nutrient concentrations, while the mechanism of these changes, although reproduced by the models as biogeochemical transports and transformations, remained largely hidden and unanalysed. Another important deficiency in prevailing modelling was the neglect of the sediments with its nutrient dynamics, which serve as a "memory" of the lake ecosystem's evolution and an important link that closes the biogeochemical cycles through the remineralization of nutrients. In the present study, we adapted and implemented the St. Petersburg Model of the Baltic Sea Eutrophication (SPBEM) with a particular focus on the analysis of the biogeochemical fluxes in Lake Onego, including its bottom sediments. Importantly for the long-term forecasting, SPBEM describes the coupled cycles of nitrogen, phosphorus, and silica, which enables its use in water bodies with spatially and temporally changing limitations by any of these nutrients. A good illustration of the SPBEM performance in waters with variable nutrient limitations is the reliable reproduction of the spatial gradients of limiting nutrients in the Neva River estuary in the easternmost Gulf of Finland [33].
Despite a wide-ranging research of the Lake Onego ecosystem [1,3,6,7,30], the drivers and components of Onego's biogeochemical cycles are still understudied. There is neither quantitative knowledge nor field data, especially on major biogeochemical fluxes, that could characterize the past and present of the Lake Onego ecosystem and could be used for reliable prognostic estimates of its future. Therefore, we implemented our ecosystem model, validated as extensively as the available data permitted, for the reconstruction of three decades of Lake Onego ecosystem dynamics with daily resolution. We also verified the model's consistency by analysing and demonstrating how it reproduces the major mechanisms that determine seasonal dynamics. The simulated three-dimensional fields of both pelagic and sediment variables, as well as the fields of most important biogeochemical fluxes, can be considered as a form of "biogeochemical reanalysis", albeit without formal data assimilation. To the best of our knowledge, this is the first use of a three-dimensional coupled hydrodynamic biogeochemical model to reconstruct past long-term biogeochemical dynamics in a large boreal lake, which presents new knowledge about the Lake Onego ecosystem. Although our model reproduces a limited number of components at the lower levels of lake ecosystem, simulated variables and fluxes can be used for further analysis together with the existing knowledge on the relationships with other, not simulated ecosystem components, for instance, between the nutrient dynamics in sediments and benthic communities [34,35].
Reanalysis with data assimilation has recently been implemented for short-term to long-term reconstructions of oceanic and marine ecosystems, with differing degrees of success [36][37][38][39]. The most important weakness of the data assimilation algorithm used in biogeochemical models based on mass balance principles is its inherent non-conservativeness. Consequently, the biogeochemical model, being reasonably calibrated during hindcast with data assimilation, may display deviant behaviour in forecasting when left unsupported by the data. The avoidance of such biases is especially important for our intention of producing prognostic estimates. Thus, the main purpose of this paper is a presentation of the 3D ecosystem model capable to a certain extent fill the historical deficit in observations of nutrient variables and, especially in estimates of the biogeochemical fluxes. According to one of the major functions of simulation modelling, we intend implementing this model as a complementary form of studies of Lake Onego ecosystem [40] providing a unifying formal platform for testing and discussing consistency of both model parameterizations and results of hydrological, hydrochemical, hydrobiological, and geochemical research. Furthermore, the model will be implemented as a major tool for a wide range of projections, from applied tasks of localization of fish farms, water intakes, and wastewater outlets to long-term large-scale ecosystem evolution under different scenarios of climate change and socio-economic development.
The entire paper is structured as follows. The model and set-up of the numerical run, including explanations on reconstruction of the initial and boundary conditions, are presented in "Model and data section". "Results and discussion" section starts with comparison between simulated selected characteristics of ice and water temperature with some available scarce observations. Scrutiny of the simulated long-term seasonal dynamics and average spatial distribution of the ecosystem variables and primary production is also accompanied by a comparison to available data, including nutrient content in the bottom sediments. The analysis of seasonal dynamics [41], starting with phenomenon of the spring phytoplankton bloom that had been missed in the earlier hydrobiological studies of Lake Onego and was revealed only in our simulation, is followed by presentation and analysis of the biogeochemical nutrient fluxes and integral nutrient budgets, also made for Lake Onego for the first time.

Model presentation
SPBEM is a coupled 3D hydrodynamical-biogeochemical model that performed reasonably well in hindcasts and scenarios of climatic changes and nutrient load reductions for the Baltic Sea [9,42]. SPBEM consists of biogeochemical and hydrodynamical modules. The former simulates the biogeochemical cycles of nitrogen, phosphorus, and silicon in the water column and bottom sediments. During over twenty years of implementation [43,44], the biogeochemical module has been thoroughly calibrated and extensively validated within BALTSEM model, plausibly reproducing ecosystem dynamics in the entire Baltic Sea [45], from the cold, annually ice-covered, almost fresh, and severely P-limited Bothnian Bay (i. e. very much Onego-like), to the warmer, mesotrophic Gulf of Finland and the Kattegat, with a single set of parameterizations and constants in both basin-wise horizontally averaged and true 3D versions [33,42,46,47]. Besides, the formulation contains no specifically salinity-determined parameterization and have already been favorably tested at Lake Ladoga [31]. Such simultaneous coverage of a wide range of ecological conditions makes us confident in application of largely the same set of formulations to Lake Onego. One of the technical advantages of using this model for the Lake Onego ecosystem is its capability to process and display not only indicators of trophic state (concentrations), but also the biogeochemical fluxes that determine the temporal and spatial dynamics of nutrient concentrations. All equations, parameterizations, coefficients, and constants of the biogeochemical module of SPBEM are presented in full detail, enabling its independent reproduction and implementation for other water bodies, in [33]. To avoid the extensive self-citation, we are not repeating those our formulations in the present paper.
The adaptation of the biogeochemical module of the SPBEM model to Lake Onego conditions consisted in redefining the autotroph variables and processes introduced in the original model [33,43]. Assuming the negligibility of nitrogen fixation in severely phosphorus-limited lake, we excluded both the process and its performers, the diazotrophic cyanobacteria functional group, from the implemented formulation. Thus, autotrophs were presented by only two variables, diatoms and non-diatoms, the later comprised all the other (summer) phytoplankton species, for example, chlorophytes, chrysophytes, and cyanobacteria. Such reformulation requested recalibration of several phytoplankton parameters, necessary also to better separate dynamics of summer "non-diatoms" from cold-water diatoms (Table A1). As was also shown by several test runs, all the other formulations, being extensively calibrated and tested in similar temperature and trophic conditions [33,44,46,47 and references therein] have not requested further fine-tuning. Besides, we omitted the silicon cycle from the model because, according to the occasional observations [48][49][50], silicate concentration during summer and autumn development of diatoms ranged within 200-400 mg Si m -3 even in the deep-water areas, and thus never became limiting. The modified scheme of the biogeochemical interactions is shown in Fig. 1.
The hydrodynamical module was built on the University of Massachusetts MITgcm model [51][52]. The model was configured to Lake Onego's bathymetry. The TKE turbulent closure scheme (GGL90 package) [53] was used to parameterize the sub-grid vertical mixing processes. The horizontal turbulent diffusion coefficients were set to be constant. Because Lake Onego is located in the subarctic zone, the SeaIce package included in the MITgcm model Three-dimensional hindcast of nitrogen and phosphorus biogeochemical dynamics in Lake Onego ecosystem… Трехмерная ретроспективная оценка биогеохимической динамики азота и фосфора в экосистеме Онежского озера… complex was used to simulate a percent of the ice coverage of grid cells. To emulate the conditions of a freshwater lake, the salinity of the ice was set to zero. The PTtracer package was used to solve the tracer equations of advectiondiffusion required for biogeochemical variables.

Study area
The fully coupled hydrophysical-biogeochemical model was run on a spherical grid with a horizontal step of 1.079′ in latitude and 2.331′ in longitude, which was approximately 2 × 2 km at the latitude of Lake Onego (Fig. 2). The z-coordinate was used with a uniform vertical step of 2 m from the surface to the bottom. Figure 2 presents also main limnic regions of the lake covering open waters and the most significant bays.

Initial and boundary conditions
The initial conditions for the hindcast simulation of the biogeochemical dynamics of Lake Onego ecosystem between 1985 and 2015 were generated during a spin-up simulation with the boundary conditions, including external nutrient inputs, repeated for 1984. With such repeating forcing, the quasi-steady state of seasonal dynamics was reached in 40 years, mainly because of a slow evolution of sediment nutrients.
Atmospheric forcing was set based on the ERA-Interim reanalysis fields [54]. The model used fields of pressure, wind speed components, air temperature, humidity, short-wave and long-wave incoming radiation, and precipitation.
Monthly river runoff was reconstructed from a water budget based on field observations for a period of 60 years [30]. The total reconstructed water discharge (Fig. 3, a) was split among 13 rivers (see Fig. 2), with empirical coefficients of the contribution of each river. The only river that drains Lake Onego into Lake Ladoga is the Svir River; its discharge was used to balance the reconstructed total river water input.
In the absence of reliable time-series of nutrient inputs entering Lake Onego from external sources (river runoff, atmospheric deposition, anthropogenic sources, etc.), reconstruction was carried out from a compilation of available published estimates. According to such estimates for 1965-2008, the rivers, wastewaters, and atmosphere contributed  Three-dimensional hindcast of nitrogen and phosphorus biogeochemical dynamics in Lake Onego ecosystem… Трехмерная ретроспективная оценка биогеохимической динамики азота и фосфора в экосистеме Онежского озера… 66 %, 23 %, and 11 % to total phosphorus (TP) input, respectively, and 67 %, 5 %, and 28 % to total nitrogen (TN) input, respectively [48]. Sabylina [55] compiled information on the total content and inorganic fractions of nutrients in river waters that was obtained in surveys performed in 1985-1986, 2001-2002, and 2007-2008. Assuming the similarity of concentrations in rivers that drain catchments with similar landscapes, Sabylina [55] provided reasoning for the aggregation of many streams and rivers into larger units flowing into Lake Onego. Correspondingly, we multiplied these concentrations by the reconstructed monthly water discharge for the indicated 13 rivers, filling the gaps over years without observations by linear interpolation between available concentration values. Dissolved organic content was further separated into labile and refractory components, assuming a 30 % bioavailability of dissolved organic nitrogen (DON) and 90 % bioavailability of dissolved organic phosphorus (DOP) [33,56].
The main sources of direct anthropogenic pressure on Lake Onego are the three industrial areas associated with the coastal towns of Petrozavodsk, Kondopoga, and Medvezhyegorsk. In contrast to the other two, discharging directly into the top of the bay, the sewage output from Petrozavodsk industrial hub is located on the open coast south of Petrozavodsk Bay (see Fig. 2). According to estimates from various publications [48,57], the total annual input from these areas varies from 170 to 250 tons of TP and was approximately 2,650 tons of TN per year. Annual nutrient deposition from the atmosphere on the lake surface ranges from 65 to 80 tons of TP and was approximately 2,260 tons of TN annually [48,57]. Since the mid-2000s, commercial fish farming has been actively developing in the waters of Lake Onego, generating annual nutrient loads of approximately 40 tons and 228 tons of TP and TN, respectively [57]. These fish farm inputs are comparable to those from the atmosphere and have been accounted for in boundary conditions since 2005.
Despite the almost trendless water discharge ( Fig. 3, a), the reconstructed external TP input showed an upward tendency, while the TN input distinctly decreased by approximately one-third (Fig. 3, b). As shows comparison of river concentrations between surveys of 2007-2008 [48] and 2019-2020 [58], these tendencies are continuing. In result, the prescribed bioavailable fraction of total nutrient inputs, comprising inorganic and labile dissolved organic compounds, reached approximately 8,300 tons N/yr and 800 tons P/yr in the 2010s, thus, decreasing the bioavailable weight N: P ratio from 13.7 in 1985/89 to 10.4 in 2011/15.

Results and Discussion
The field data availability (its coverage, regularity, and frequency) for Lake Onego is not sufficient to allow extensive model-data comparisons similar to the ones made, for instance, for the Baltic Sea [33,44,46,47] or Great American lakes [24 and references therein]. Therefore, omitting the formal validation that usually precedes further modelling analysis, we immediately start next Section with the results of simulation, on the way involving in the analysis as much as possible and whatever data are available for Lake Onego. In addition to this highly insufficient information, we also used some relevant material from other boreal and even temperate lakes, situated in the differing hydrometeorological conditions and impacted by the differing land cover patterns and land use practices at their watersheds. Despite these unsurprisingly emerging differences, we justify such involvement by the following considerations.
Our mechanistic model is based on a mass balance approach, describes internal biogeochemical cycles, and accounts for external sources and sinks (imports and exports), either prescribed as forcing functions or computed according to formulations. Consequently, all the simulated fluxes and concentrations (cf. Fig. 1) are strongly deterministically coupled and thus, confined. Therefore, their reliability should be judged by a simultaneous fitting of many simulated fluxes and concentrations in the known ranges reported for both Onego and similar boreal oligotrophic lakes. For example (see Table 1 in [41]), the sedimentation of organic matter (OM) cannot be much higher than nutrient uptake during primary production of OM simulated with the prescribed phytoplankton specific growth rates (Table A1), which plausibility is justified by a reasonable PP validation (see Table 1 and Fig. 10 below). Similarly, the sediment release (and denitrification) of nutrients cannot be much higher or lower than the sedimentation flux, otherwise the difference would cause fast accumulation or depletion of sediment nutrients that so far has not been reported from observations. Furthermore, the plausibility of simulated rates is estimated by a comparison to sediment rates from similar environments. This can be said about all the other processes in Fig. 1 above and Table 1 in [41]. That's why we consider the quantitative information from other lakes combined with Onego data as quite relevant and fully justified.

Important hydrophysical characteristics
The implementation of MITgem model for simulation of transport processes in the framework of SPBEM ecosystem model was justified by its successful application to large lakes such as Ladoga [31], Michigan [59,60], and Superior [61]. At the same time, we presume that 30-year atmospheric forcing based on the ERA-Interim reanalysis with a basic spatial resolution of 80 km and assimilating information from a comparatively sparser network of meteorological stations in this region is hardly suitable for reliable sequential reproduction of transient synoptic situations of 5-7 days duration over 30 simulated years.
In this paper, aimed at the long-term and seasonal biogeochemical dynamics, we also omitted the analysis of the most hydrophysical characteristics and present model-data comparison only for the ice coverage and water temperature as important integral indicators of the hydro-thermodynamics that significantly affect the ecosystem dynamics [62].
The realistic simulation of the ice cover dynamics, especially of its melting phase is a prerequisite for a timely reproduction of the phytoplankton spring bloom commencement, which in the model is determined by the light availability increasing due to the ice melt (see Fig. 9 in Sect. 3.2). A comparison of simulated dynamics to estimates obtained within visual-instrumental approach from satellite data by a semi-supervised algorithm [7] shows that the model captures both the timing and duration of seasonal ice phases rather accurately (Fig. 4), with the coefficient of linear correlation R = 0.99 and RMSE = 6.3 %. Most closely coincide the onset of freezing, observed from mid-December to early January, and the final ice clearing occurring from mid to late May. Calculation of the mean and standard deviation for the seasonal variation from 0 to 100 % has hardly much of interest.
In addition to its significance in Lake Onego's hydrophysics, water temperature is a crucial factor in ecosystem dynamics, including its effects on the phenology and intensity of biogeochemical fluxes. Available water temperature measurements in Lake Onego were scarce and irregular. The most frequent observations were made in the summer months and occasionally in May and September. Bearing in mind a likely irregular asynchrony between the day-by-day simulated temperature fields and those that had actually been occurring in reality, we did not attempt the detailed pair-wise comparison of 532 measurements made in June-August over 15 years of 1992-2007, a priory considering such comparison confusing and misleading rather than validating. Instead, these measurements of the surface temperature were pooled together regardless of the sampling location. The model results for the same months and years comprised the daily average water temperatures from all of the computational grid nodes. Over 70 % of all measurements were collected from the coastal areas and bays; thus, the contribution of the expansive and colder open waters to the observational statistics was lower than their contribution to statistics calculated from simulation. Taking this expected bias into consideration, the simulated water temperature (mean ± standard deviation 11.85 ± 3.92 °C; median 13.15 °C) matches the observations (13.05 ± 4.82 °C, median 14.40 °C) rather well. Model validation for the underrepresented open areas was made for three stations (cf. Fig.  2), covered with regular observations from 1985 to 1989, while later only episodic measurements were made. The best model-data comparability was found during the periods of heating and cooling of the lake, while the larger discrepancies occur during summer maxima (Fig. 5). Apparently, with a grid cell height of 2 m and daily averaging we cannot expect precise reproduction of the diurnal extremes in the surface temperature measured in the upper 0.5 m water layer. Such causes of the discrepancy could also partly explain the summer bias in the total model-data comparison above.
As shows a comparison between average vertical profiles of the water temperature, simulated and reconstructed from measurements (Fig. 6, cf. Fig. 2, a and Fig. 3, a in [41]), the summer thermic vertical structure is also simulated rather reasonably. Unfortunately, the further statistical analysis of model-data comparability is prevented by both paucity and irregularity of observations scattered over two decades with some years missed entirely.
The concept of "biological summer" can be characterized by its duration (BSD) and average water temperature (BST) and is often used in hydrobiological studies. At Lake Onego, the day of the transition of the surface water layer Three-dimensional hindcast of nitrogen and phosphorus biogeochemical dynamics in Lake Onego ecosystem… Трехмерная ретроспективная оценка биогеохимической динамики азота и фосфора в экосистеме Онежского озера… temperature to over a threshold of 10 °C was set as the phenological indicator of both the onset and end of "biological summer" [1,63]. In the simulation, a long-term tendency of increasing BST was clearly seen through significant interannual variations (Fig. 7). This tendency agrees well with both simulations and observations in the European boreal zone [30,42,64,65].  The BSD in Lake Onego calculated from the simulation was 92 ± 9 days, within the range of 71-107 days. Based on sparse field data, BSD in the shallow areas (100-110 days) was markedly longer than that in open deep water (85-90 days) [66]. In simulation, BSD in the area shallower 30 m was 99 ± 10 days comparing to BSD in deeper areas with 93 ± 10 days. The prolongation of "biological summer" occurred chiefly because of its extension into autumn, from September 27 (mean date for the first 5-year period) to October 5 (mean date for the last 5-year period). A similar asymmetry between trends of spring warming and autumn cooling has also been found in the Baltic Sea [64] and Karelian lakes [66].

Ecosystem variables
Daily long-term dynamics of phytoplankton primary production and ecosystem variables averaged over the entire Lake Onego are presented in Fig. 8. Assuming the classical Redfield molar ratio of C: N: P=106:16:1 being valid for the freshwater phytoplankton, including boreal oligotrophic lakes [67][68][69][70][71], the net primary production of diatom and non-diatom phytoplankton simulated in nitrogen weight units were converted into carbon weight units with a factor of 5.7 mg C mg -1 N. Presented concentrations of inorganic and total nutrients were averaged over the whole water volume. Consequently, the dissolved inorganic phosphorus (DIP), comprising summer phosphorus accumulation in the hypolimnion (cf. Fig. 3, b in [41]), was never fully depleted. The lake-wide averaged winter maximum values (Fig. 8, c, e) multiplied by the Lake Onego model volume (297 km 3 ) could be used to conveniently estimate total nutrient stocks. Diatom and non-diatom phytoplankton biomass, simulated as nitrogen units, were recalculated into wet weight units assuming a nitrogen content of 0.5 % and 1 % of the biomass of diatoms and non-diatoms, respectively [68,72]. For zooplankton biomass, a nitrogen content of 1 % was assumed [73].
The simulated long-term dynamics steadily reproduce a distinctive seasonal pattern: the strong phytoplankton spring bloom transitions into the summer quasi-steady-state phase followed by minor autumn blooming. The shortterm interannual variations were more pronounced in recent years in terms of primary production and biotic variables with higher plankton biomass (Fig. 8, a, b, d). The abiotic variables showed some apparent tendencies rather than distinctive periods (Fig. 8, c, e, f). The simulated average concentrations of TP (7-13 mg P m -3 ) and DIP (2-7 mg P m -3 ) classify Lake Onego as bordering the oligotrophic and mesotrophic trophic states, according to some classifications [58,74,75]. Considerable decreases in TN and dissolved inorganic nitrogen (DIN) concentrations were mainly related to the reduction of external inputs (Fig. 3, b). Such simulated quasi-stability of TP concentrations and clearly decreasing DIN concentrations (Fig. 8, c, e) is validated in [58], who found at the surface of pelagic part of Lake Onego a statistically significant decreasing trend of the DIN: TP weight ratio (mg N mg -1 P) from 33.7 in 1992-1995 to 23.7 in 2019-2021, the later estimated from the field surveys made in September 2019, June and August 2020. The observed decrease is well comparable to a simulated decrease of DIN: TP ratio from 36.5 ± 1.9 in 1992-1995 to 21.2 ± 1.1 in 2011-2015; these mean ± S.D. values are computed for the surface layer in I-IV limnic areas (cf. Fig. 2) over biological summer (cf. Figs. 7 and 11, a). Such a rapid response to the changing external inputs can be explained by the short nutrient residence times for the bioavailable fractions that can be estimated for the water body from winter maxima and external inputs, approximately 8 and 2 years for nitrogen and phosphorus, respectively. Annual integrals of the simulated phytoplankton primary production of 17.0-20.6 g C m -2 yr -1 were almost invariable and matched the global mean estimate of 20 g C m -2 yr -1 for the large lakes situated to the north of 60° N [76] but constituted only one-fifth of 94 g C m -2 yr -1 that was estimated for Lake Superior [77]. A small rise in plankton biomass (cf. Fig. 8, b, d) and annual production from the start toward the end of the simulation (from 18.5 to 20.2 g C m -2 yr -1 ) could be explained by the increased spring production, integrated over the time interval from the beginning of the year until the onset of biological summer, from 10.0 to 12.0 g C m -2 yr -1 . In turn, this biotic growth occurred because of increased surface winter DIP accumulation (maximum values increased from 5.9 to 6.5 mg P m -3 ), which was caused by a combination of increased external loads (cf. Fig. 3, b) and internal total P recycling. The later comprised pelagic regeneration due to zooplankton excretion and organic phosphorus remineralization, as well as phosphate release from the sediments, from 457 to 492 mg P m -2 yr -1 . A 1-week prolongation of the "biological summer" (cf. Fig. 4) insignificantly affected the small summer PP increase from 6.4 to 6.7 g C m -2 yr -1 .

Spatial inhomogeneity
In contrast to the nearly invariable average seasonal dynamics in Fig. 8, the spatial distributions of variables and fluxes, which were largely determined by their proximity to external nutrient sources, were highly inhomogeneous (Fig. 9, note the logarithmic scale). Because spatial gradients in Fig. 9 are significantly larger than the long-term tendencies of lake-wide averages (Fig. 8), we can also compare simulated concentrations (Fig. 9, b, c) to those measured in the later surveys [49,58]. In 2016 and 2019-2021, average DIN and TN concentrations of about 200 mg DIN-N m -3 and 400-500 mg TN-N m -3 , respectively in the Petrozavodsk Bay were expectedly (cf . Fig 8, e) slightly lower than averaged from simulation for 1985-2015 (Fig. 9, c), but, being non-limiting nutrients, remained almost the same in the Central Onego. In contrast, average DIP and TP concentrations of 4-9 mg DIP-P m -3 and about 20 mg TP-P m -3 in the Petrozavodsk Bay dropped under the detection limit of phosphates and to 6-8 mg TP-P m -3 in the Central Onego. Thus, the observed distributions are reasonably comparable to the simulated ones.
In the zones of elevated winter DIP concentrations (Fig. 9, b), annual primary production rates exceeding 25 g C m -2 yr -1 put large areas on the border of the oligotrophic status, while eutrophic areas were observed near the Petrozavodsk and Kondopoga industrial centres as well as along the coast south of Petrozavodsk Bay, where PP increases up to 120-150 g C m -2 yr -1 [2,49,78]. The simulated inhomogeneous distribution (Fig. 9) matched the estimates obtained from the direct measurements (Table 1).
In general, the summer offshore values were expectedly lower than the lake-wide mean primary production rates estimated for more eutrophic lakes Huron, Michigan, and Superior for 2010-2013 (216, 259, and 228 mg C m -2 day -1 , respectively); although primary production in all depth zones (shallow, mid, and deep) were similar across the lakes [79]. Fig. 9. Simulated long-term (1985-2015) averages of (a) annual primary production (g C m -2 yr -1 ); April surface concentration of (b) dissolved inorganic phosphorous (DIP; mg P m -3 ) and (c) dissolved inorganic nitrogen (DIN; mg N m -3 ) (note the logarithmic scale); average seasonal dynamics of (d) annual primary production (g C m -2 yr -1 ) and (e) phytoplankton biomass (g ww m -3 ) in the different limnic regions of Lake Onego; and averages of (f) benthic phosphorus (BEP; mg P m -2 ) and (g) benthic nitrogen (BEN; mg N m -2 ) with location of sites chosen for comparison with observations Three-dimensional hindcast of nitrogen and phosphorus biogeochemical dynamics in Lake Onego ecosystem…

Трехмерная ретроспективная оценка биогеохимической динамики азота и фосфора в экосистеме Онежского озера…
The general levels and seasonal patterns of PP and phytoplankton biomass in the bays and open waters also differed ( Fig. 9, d, e). A pronounced spring maximum in the oligotrophic central and southern Onego was followed by a weak autumn blooming after a long summer minimum. In contrast, there were two equal peaks of primary production in Petrozavodsk Bay at the end of May and July. Similar seasonal dynamics were simulated in Kondopoga Bay, although both the spring and autumn peaks were somewhat delayed. Except for the "revealed" in this simulation spring initiation of the vegetation season, these spatial and temporal developments were in good agreement with both available summer measurements of primary production (Fig. 10) and the general phenology of the Lake Onego ecosystem deduced from field observations [1,49,63,74,80,81]. Table 1 Phytoplankton primary production in Lake Onego according to observations (from Table 1 Fig. 10. Model-data box-and-whisker plots of the phytoplankton primary production (mg C m -2 day -1 ) in different limnic areas of Lake Onego (cf. Fig. 2). Note paucity and irregularity of observations The simulated spatial and temporal dynamics of phytoplankton biomass could also be validated by remotely sensed chlorophyll distributions [82][83][84]. However, we would not do such comparison here for two reasons. First of all, however important are synchronous chlorophyll fields obtained by remote sensing as a relative measure of phytoplankton biomass and its dynamics [85] they cannot be reliably used for estimation of in situ biomass without special empirical conversions by the prescribed C: Chl ratio, which at the temperate and boreal latitudes vary within an order of magnitude, from 20 g C: g Chl-a in spring to 200 g C: g Chl-a in summer [45 and references therein, 86]. In principle, it could be possible to qualitatively compare typical spatial patterns, including succession of phenological phases and its interannual variations. However, in contrast to Lake Ladoga [87,88], the specific retrieval algorithm for Lake Onego Chlorophyll has not been developed yet. Neglecting of such algorithm [89] in these boreal, highly brownified waters [3] would lead to questionable, if not even unreliable results.
The spatial distribution of primary productivity was well reflected in the sediments (Fig. 8, f, g). In the model, the sediment variables of benthic phosphorus (BEP) and nitrogen (BEN) were described according to a vertically integrated dynamic approach in areal units of mg P (N) m -2 for the biogeochemically active surface layer of unspecified thickness [33,43,90,91]. In geochemistry, nutrient content is usually calculated as a percentage ratio of the mass of phosphorus and nitrogen contained in a gram of ash-free dry sediment [1,34]. For model-data comparison ( Table 2), we chose four representative sites (see Fig. 8, g) and converted available measurements characterizing surface sediment layer of 5 cm thickness [34,74] into model units, using the bulk density values obtained from similar seascapes [92]. Wide ranges of published sediment characteristics compiled in Table 2 from maps in [74, рр. 100-101] and a few numbers in [34, Table 1], are based on samples taken over many years from the mosaic distribution of real sediment types, thus hardly guaranteeing the exact positioning, sufficient homogeneity and reproducibility of measurements at indicated locations. On the other hand, our simplified and spatially invariable sediment parameterizations totally disregard the diversity of sediment types, while registering a dynamic balance between the source (sedimentation) and the sinks (nutrient release and burial), and closing nutrient cycles by regeneration (see Table 1 in [41]). Taking into account these interpretations, a comparability between simulated values and estimates, recalculated from measurements, could be considered plausible.
According to the simulation, the total stocks of N and P in the surface sediments of Lake Onego, which cover a bottom area of 9,266 km 2 and are actively involved in the contemporary biogeochemical cycling of the lake, amount up to 80,000 and 40,000 tons, respectively. The low weight N: P ratio of two can be attributed to the low nitrogen content (ca. 0.45 % N of dry mass) in the contemporary deep-water sediments of oligotrophic lake, which, on the other hand, are highly enriched with phosphorus (ca. 0.25 % P of dry mass) associated with iron-humic complexes [3,34]. These values are close to those of the contemporary Lake Winnipeg's sediment nutrient content of 0.4 % N [93] and 0.12-0.26 % P [94]. Regarding such integral characteristics and considering statements of the occurrence of many-fold increases in P sediment content during recent decades [1], one should be concerned as to where such increases occurred and which areas were enveloped. Less than 1,000 tons of annual TP load (cf. Fig. 3) would not be sufficient for any substantial alteration of the present-day TP integral sediment stocks. Table 2 Sediment characteristics and nutrient content measured at the sampling sites (see Fig. 8, g)  Three-dimensional hindcast of nitrogen and phosphorus biogeochemical dynamics in Lake Onego ecosystem… Трехмерная ретроспективная оценка биогеохимической динамики азота и фосфора в экосистеме Онежского озера…

Conclusions
1. Despite a long history of extensive research, the drivers and components of Lake Onego's biogeochemical cycles are still insufficiently studied. Especially damaging such a deficit of knowledge could be to our capability of forecasting possible changes in the lake ecosystems in response to natural variations and anthropogenic impacts. To enable a quantitative description of the past, present, and, eventually, the future state of Lake Onego, we adapted and successfully implemented a three-dimensional biogeochemical model, considering the obtained results as a form of reanalysis.
2. The model was used to reconstruct three decades of Lake Onego ecosystem dynamics with daily resolution. Although the paucity of observations did not allow either formal model validation or the statistical demonstration of its skills, the comparison to all available information from Lake Onego and to a range of published estimates for other large boreal lakes led us to believe in plausibility of simulation.
3. The basin-wide biogeochemical reanalysis also eliminated a historical disparity in attention and, hence, in observations and knowledge, between bays (especially Petrozavodsk and Kondopoga bays) and vast open-water areas. Bays occupy relatively small bottom areas and water volumes and, consequently, are responsible for only a small fraction of biogeochemical phenomena and fluxes.

Acknowledgments
We are grateful to N.M. Kalinkina and T.F. Tekanova for their helpful discussions on the purpose and results of this study. We thank A.F. Balaganskyi for the river water discharge data.

Funding
The author A.I. conducted the present study within the framework of state assignment (Theme No. FMWE-2021-0014). The author O.S. was supported by the Swedish Agency for Marine and Water Management through their grant 1:11-Measures for the marine and water environment. The author N.F. conducted the present study within the framework of state assignment (Theme No. FMEN-2021-0007).