Dynamic modelling provides new insights into development and maintenance of Lake Kivu's density stratification

Lake Kivu is a 485 m deep, Central-East African rift lake with huge amounts of carbon dioxide and methane dissolved in its stably stratified deep waters. In view of future large-scale methane extraction, one-dimensional numerical modelling is an important and computationally inexpensive tool to analyze the evolution of stratification and the content of gases in Lake Kivu. For this purpose, we coupled the physical lake model Simstrat to the biogeochemical library AED2. Compared to an earlier modelling approach, this coupled approach offers several key improvements, most importantly the dynamic evaluation of mixing processes over the whole water column, including a parameterization for double-diffusive transport, and the density-dependent stratification of groundwater inflows. The coupled model successfully reproduces today’s near steady-state of Lake Kivu, and we demonstrate that a complete mixing event ~2000 years ago is compatible with today’s physical and biogeo-


Introduction
Lake Kivu is a large (2386 km 2 ) and deep (485 m) tropical rift lake, situated on the boundary between Rwanda and the Democratic Republic of the Congo (DRC), directly south of the Virunga volcano chain. It is fed by numerous small streams, and discharges, via the Ruzizi River, into Lake Tanganyika (Fig. 1a). In addition to surface streams, ~45% of the inflow into Lake Kivu is provided by subaqueous groundwater sources (Schmid and Wüest, 2012), some of which were identified in the northern part of the lake using temperature and salinity profiles (Ross et al., 2015a and Fig. 1a). Part of this intruding groundwater is hydrothermal, meaning that it is warm, salty and rich in carbon dioxide (CO 2 ). Due to their high density, these hydrothermal sources have a tendency to plunge before eventually stratifying (Ross et al., 2015a). In addition to the hydrothermal groundwater, two cooler and less salty sources were proposed based on numerical modelling (Schmid et al., 2005), and subsequently identified by conductivity and temperature profiling (Ross et al., 2015a). The model of Schmid et al. (2005) suggested that the discharge of the cooler sources is one magnitude larger than the hydrothermal discharge, and that they explain the strong thermo-and chemoclines at around 190 and 250 m (see Fig. 1b for an overview).
The groundwater intrusion into Lake Kivu has two major consequences for the whole water column: i) a continuous upwelling of up to 1 m yr − 1 (Pasche et al., 2009), and ii) a strong density stratification and therefore meromixis below 60 m. Below the main chemocline at 250 m, the lake water probably has not been in contact with the atmosphere for up to 1000 years and hence, the inflowing CO 2 accumulated to concentrations of up to almost 100 mmol L − 1 (Tietze, 1978;Schmid et al., 2005;Bärenbold et al., 2020a). Furthermore, the decomposition of settling organic matter and the reduction of CO 2 (Pasche et al., 2011) also led to the accumulation of methane (CH 4 ) in the deep water. The current total content of CO 2 and CH 4 is estimated to ~285 and ~62 km 3 , respectively (Bärenbold et al., 2020a). The CH 4 reservoir in Lake Kivu is a valuable resource for the neighboring countries Rwanda and DRC, but also a looming danger for the surrounding population. In fact, the smaller gas-rich lakes Nyos and Monoun in Cameroon have experienced limnic gas eruptions in the past (Kling et al., 1987;Sigurdsson et al., 1987), and 1746 and 37 people were killed by asphyxiation, respectively. In 2002, the volcano Nyiragongo erupted, and lava flowed into Lake Kivu. It was feared that the hot lava could sink to the deeper strata of Lake Kivu, and heat up water to an extent that it could rise to a level where the hydrostatic pressure is no longer sufficient to keep the gases trapped. However, Lorke et al. (2004) showed that the influence of the lava on the lake was minimal, and they concluded that subaerial volcano eruptions are not likely to trigger a limnic eruption of Lake Kivu. In May 2021, Nyiragongo erupted again but without apparently disturbing the lake's stratification, although the observations indicate the formation of a new dyke beneath the lake (Jones, 2021). Even if these recent volcanic activities did not significantly affect the lake, they highlight how geologically dynamic and potentially unstable this region is.
In order to prevent further gas eruptions, Lakes Nyos and Monoun, as well as Kabuno Bay, which is a small, 150 m deep subbasin of Lake Kivu with only low CH 4 content, are degassed artificially today. Degassing is performed using long pipes which transfer gas-rich deep water by selfsiphoning to the lake surface (Kusakabe, 2017;Halbwachs et al., 2020). In the main basin of Lake Kivu, an artificial decrease of gas concentrations is beneficial for two reasons: i) decreasing gas pressure increases the safety margin, i.e. the distance that a water parcel can rise without outgassing during a mixing event, and ii) the extracted CH 4 can provide energy for both Rwanda and the DRC. In December 2015, the U. S. company ContourGlobal started the operation of the first large-scale CH 4 extraction facility "KivuWatt" with an electric capacity of 26 MW. KivuWatt extracts water from a depth of 350 m, removes most of the CH 4 and reinjects the partially degassed water at ~240 m, right above the main chemocline. The projected expansion of gas extraction at Lake Kivu (Expert Working Group on Lake Kivu, 2009) would lead to a major redistribution of warm, salty and nutrient-rich deep water within the lake. If not managed carefully, such redistribution of water has the potential to change stratification, gas distribution and nutrient availability in Lake Kivu during the next decades. For example, if the reinjected water ends up too close to the lake surface it could provoke harmful events like outgassing or algal blooms. Regular monitoring of stratification and gas concentrations combined with numerical modelling for forecasting are key to prevent such negative consequences (see Schmid et al., 2019 for a detailed discussion on monitoring).
Due to the almost perfect horizontal homogeneity in the stably stratified deep waters, i.e. below ~60 m (see Schmid and Wüest, 2012 for a detailed evaluation of horizontal homogeneity), one-dimensional modelling is an adequate tool to assess the changes discussed above. Note, however, that such modelling is not capable to evaluate inherently three-dimensional processes like surface layer dynamics across basins or the dynamics in the vicinity of groundwater sources or water extraction and reinjection by power plants. 1D modelling of Lake Kivu's water and gas dynamics was previously performed by Schmid et al. (2005) using Aquasim, a tool for simulating diverse aquatic systems (Reichert, 1994). This model was later improved, and used to predict changes in the lake's stratification, gas content and nutrient concentrations for different methane extraction scenarios Schmid et al., 2019). However, the model by Schmid et al. (2005) included several Fig. 1. a) Map of Lake Kivu with bathymetry (contour lines every 100 m), surface tributaries and delineation of the watershed. The part of the watershed, which is not drained by surface tributaries is colored in dark-grey and accounts for 14% of the total watershed. The 3 H sampling site and the Lake Kivu weather station are indicated, as well as locations, where temperature and salinity spikes indicate the intrusion of groundwater (according to Ross et al., 2015a). Warm and salty, and cooler and fresher groundwater are marked with red and blue dots, respectively. b) Schematic cross-section of Lake Kivu with temperature, salinity, gas profiles and groundwater sources (reproduced from Schmid et al. (2021)). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) simplifications and assumptions that could lead to erroneous conclusions if applied to long-term transient simulations. Most importantly, the model could not reproduce feedbacks between changes in lake stratification, and the intrusion depths of the subaquatic springs and turbulent diffusive transport was assumed to be constant in the top 120 m. Both of these simplifications are not adequate when doing transient simulations.
Here, we propose a new modelling approach using the physical onedimensional lake model Simstrat (Goudsmit et al., 2002), which we coupled to the biogeochemical library AED2. An uncalibrated version of Simstrat had already been used to model Lake Kivu's surface waters in a lake model intercomparison study and performed reasonably well (Thiery et al., 2014). In short, the new modelling approach has the following advantages: i) diffusive transport between 0 and 120 m is dynamically computed using a k-ε turbulence model, ii) the transport through double-diffusive staircases below 120 m is parameterized based on recent field observations, and iii) the groundwater inflows are free to plunge and rise as a function of their density compared to the density of the ambient lake water. Due to these advantages, our new model can replace the model developed by Schmid et al. (2005) for longer, transient simulations of Lake Kivu, i.e. for example to model the evolution of gas concentrations and stratification due to large-scale gas extraction projects.
Although our model could be readily used for such simulations, we do not focus on the changes induced by gas extraction in this work. Instead, we use the model to address three specific scientific questions. Firstly, Ross et al. (2015a) suggested that the cool and relatively fresh sources at ~190 m and ~250 m might mainly consist of recently infiltrated rainwater at the northern shore of the lake. There, rainwater probably quickly percolates through the porous volcanic rock. We can test this hypothesis using our model and recently acquired tritium ( 3 H) measurements. Secondly, based on observations in the sedimentary record, several authors (Ross et al., 2015b;Votava et al., 2017;Uveges et al., 2020) suggest that Lake Kivu experienced a large-scale mixing event ~1000 years ago with nutrient-and gas-rich water reaching the surface layers and creating algae blooms. To assess whether a complete mixing event could have happened at that time, we use our model to determine whether and within which time period today's stratification and gas concentrations could have developed from a homogeneous (i.e. completely mixed) and degassed lake. Finally, we previously observed large depletions of atmospheric noble gases in the deep water of Lake Kivu (Bärenbold et al., 2020b). Using our model, we test our hypothesis that the inflow of depleted groundwater is the most probable cause of these depletions.

Basic model description
The lake model Simstrat was originally developed at Eawag by Goudsmit et al. (2002) and it is available on github (https://github. com/Eawag-AppliedSystemAnalysis/Simstrat). Simstrat solves the one-dimensional hydrodynamic equations on a staggered grid, using an implicit Euler numerical scheme. Turbulent diffusion is resolved by a k-ε closure, which means that the turbulent kinetic energy k and its dissipation ε are explicitly modelled as state variables. The main sources of turbulent kinetic energy are direct wind forcing and buoyant instability of the water column. However, Simstrat also allows for turbulent energy production by indirect wind forcing, i.e. the dissipation of internal waves (Goudsmit et al., 2002). A certain fraction of the wind energy (represented by the seiche energy parameter α s ) is redistributed to internal seiche energy, which is subsequently dissipated through friction at the boundaries and supplies mixing energy distributed over the entire depth of the lake. As discussed by Goudsmit et al. (2002), this process is a major source of energy for turbulent mixing in the hypolimnia of lakes, and neglecting it leads to systematic underestimation of the vertical turbulent mixing.
The original version of Simstrat has been updated several times in the past. The most important upgrades were the use of a different heat flux parameterization (Schmid and Köster, 2016), the implementation of density-driven plunging of inflows (Gaudard et al., 2017), and the possibility to simulate lake ice cover (Gaudard et al., 2019). In order to explicitly simulate the formation of CH 4 in Lake Kivu, we coupled Simstrat v2.2 with the biogeochemical model library AED2 v1.3.5 (https://github.com/AquaticEcoDynamics/libaed2), which is developed by the University of Western Australia. AED2 consists of a set of biogeochemical modules (carbon, nitrogen, phosphorus, phytoplankton etc.), which each can be switched on and off, or adjusted according to specific needs. In this work, we use the carbon and oxygen modules of AED2 to simulate the carbon cycle, including the formation and consumption of CH 4 .
The calculation of the radiation budget in Simstrat requires some basic physical parameters. Water albedo is determined as a function of latitude and seasonality by using a linear interpolation of the modified values of Grishchenko in Cogley (1979), which are implemented in Simstrat v2.2. We use a standard value of 0.35 for the fraction of shortwave radiation absorbed as heat in the surface boundary layer. The remaining shortwave radiation is absorbed in each water layer with an extinction coefficient of 0.27 m -1 (Thiery et al., 2014). In addition, we need to define the geothermal heat flux, which is not well constrained in the Lake Kivu region. Degens et al. (1973) recorded a few measurements in the sediment which indicate a range of 0.016-0.18 W m − 2 . For our simulations, we chose a value of 0.13 W m − 2 , which is well within this range and leads to the best agreement between observed and simulated deep water temperatures. All Simstrat settings and parameters are summarized in the supporting information (Tables S1 and S2).
In the AED2 carbon module, the CH 4 sediment flux P CH4 in mmol m − 2 d − 1 is modelled taking into account oxygen inhibition and temperature dependence: where F sed,CH4 is the maximum sediment CH 4 flux at 20 • C, K sed,CH4 is the half saturation constant for the inhibition of CH 4 production by oxygen, [O 2 ] is the oxygen concentration and θ sed,CH4 is the Arrhenius temperature multiplier for CH 4 production. Similarly, CH 4 oxidation by methanotrophs Ox CH4 in mmol m − 3 d − 1 is modelled as where R ox,CH4 is the maximum oxidation rate at 20 • C, K ox,CH4 is the half saturation constant for CH 4 oxidation and vT ox,CH4 is the Arrhenius temperature multiplier for CH 4 oxidation. Due to the narrow range of temperature in Lake Kivu (23-26 • C), we neglect the temperature dependence in both equations, and equation (1) then reduces to a constant production rate F sed,CH4 in the anoxic waters below ~60 m (see Table S3 for numerical values). CH 4 oxidation by methanotrophs (Equation (2)) on the other hand happens right at the oxycline and therefore is a function of oxycline depth. Note that for computational efficiency, we did not include the simulation of nutrients and phytoplankton, so oxygen consumption is underestimated and the oxycline depth is probably overestimated, which slightly affects the depth at which CH 4 disappears.

Modifications of the model for Lake Kivu
Due to exceptionally high CH 4 and CO 2 concentrations and the presence of double-diffusive staircases, the simulation of Lake Kivu requires some changes to the original Simstrat-AED2. This adjusted version is called Kivu-Simstrat (see Software Availability section), and it is used throughout this work. The respective changes are described in this section and depicted in the modelling flow chart in Fig. 2.

Parameterization of double-diffusive transport
One of the unique features of Lake Kivu are the double-diffusive staircases, which exist below a depth of 120 m and all the way down to the lake bottom at 485 m (Newman, 1976;Sommer et al., 2013). These staircases consist of ~300 thin mixed layers with an average thickness of ~70 cm and sharp gradients in between (Sommer et al., 2013). Double diffusion is caused by the coexistence of two substances with reciprocal effect on density and very different molecular diffusion coefficients. In Lake Kivu, the temperature increase with depth destabilizes the water column, whereas the increase with depth of salinity and CO 2 (somewhat weakened by the CH 4 increase) stabilizes it. Consequently, the strong stratification suppresses turbulent diffusion below 120 m, but molecular transport is enhanced by the step-wise profile compared to a smooth profile. By measuring the temperature and salinity gradients between double-diffusive layers, Sommer et al. (2019) concluded that the apparent temperature diffusivity (calculated as the heat flux divided by the smoothed temperature gradient) is usually between 10 − 6 and 10 − 5 m 2 s − 1 . Therefore, double diffusion in Lake Kivu enhances the heat flux by a factor of ~10-100 compared to a smooth profile with only molecular heat diffusion (~1.4×10 − 7 m 2 s − 1 ). Conversely, the salinity flux, which has a molecular diffusion coefficient of 1.2×10 − 9 m 2 s − 1 , is increased by a factor of 30-300 due to the steeper interface gradients (Sommer et al., 2013).
Double-diffusive convection is an inherently three-dimensional process and therefore not readily replicable in one-dimensional models. We initially chose a very high depth resolution (grid cells <1 cm) in order to correctly represent the thin interfaces. However, the simulated interfaces were too steep and the mixed layers too large and thus, the resulting double-diffusive heat flux was on the order of 0.3 W m − 2 instead of 0.1-0.15 W m − 2 (determined by Sommer et al., 2013). Because of this disagreement with data and because such a high spatial resolution is computationally very expensive, we decided to resort to a coarser resolution, which produces a smoothed profile between 120 and 485 m without double-diffusive staircases. We found that a resolution of 2 m is sufficient to correctly reproduce the smoothed profile and that it allows simulations over thousands of years within a few days time. The double-diffusive transport of heat and salinity through this smoothed profile is expressed as apparent diffusion, and we parameterize it as a function of stability N 2 [s − 2 ], which is given by where g is the gravitational acceleration, ρ is water density, ρ 0 is the density of freshwater (1 kg L − 1 ) and z is the vertical axis. N 2 is thus proportional to the density gradient and positive values indicate a stably stratified water column, while negative values represent unstable conditions and thus convective instability. Fig. 3a and b shows that the relationships between the apparent diffusivities of temperature K T,app and salinity K S,app (data from Sommer et al., 2019) and N 2 correlate well in the log-log domain.
The diffusivities can therefore be expressed as power laws according to where F = 0.74 is a dimensionless correction factor. This correction factor accounts for the fact that the parameterization tends to lead to steeper temperature and salinity gradients in the model compared to observations of Sommer et al. (2019) and thus slightly overestimates vertical transport. At the logarithmic scale, this correction factor consists of an offset of the y-axis, but preserves the slope of the regression line. At high values of N 2 , the lower limit for apparent diffusion is the respective molecular diffusion (Fig. 3). Conversely, when N 2 is very close to zero (i.e. <10 − 6 s − 2 ), our parameterization simulates convective instability and can produce very high diffusivities. To avoid unphysically high mixing rates, we use a threshold of 10 − 10 s − 2 for N 2 , which yields upper bounds of ~8×10 − 4 and ~2×10 − 4 m 2 s − 1 for the diffusivities of temperature and salinity, respectively. The molecular diffusion coefficients of CO 2 , CH 4 (Sommer et al., 2013), noble gases (Wise and Houghton, 1968) and water isotopes 3 H and 2 H (Mills, 1973) are of similar order of magnitude compared to salinity, and therefore the double-diffusive transport of each of these variables is approximated as where D i,mol and D S,mol are the molecular diffusion coefficients of variable i and salinity, respectively.

Influence of dissolved gases on density
Due to their high concentrations, salinity and dissolved gases need to be taken into account in water density calculations for Lake Kivu. In our model, we compute water density ρ as a function of temperature T according to Chen and Millero (1986) and we add the effects of salinity S, CO 2 and CH 4 similarly to Schmid et al. (2004): where β S , β CO2 and β CH4 are the contraction coefficients of salinity, CO 2 and CH 4 respectively (see Schmid et al., 2004 for numerical values). In equation (7) and throughout this work, S is computed from conductivity measurements by using ionic composition in Lake Kivu according to Wüest et al. (1996).

Total alkalinity
In the AED2 carbon module, total alkalinity (TA) is derived as a regression of dissolved inorganic carbon (DIC) and S with different builtin regression functions to choose from. For Lake Kivu, Wüest et al., 2009 showed that TA in mmol m − 3 is a linear function of S in g kg − 1 and suggested the following relation: We implemented this option in AED2 and used it throughout our simulations.

Lakeatmosphere interactions of water isotopes ( 3 H and 2 H) and noble gases
In order to model the water isotopes 3 H and 2 H, we take into account their interaction with the atmosphere using the Craig-Gordon model (see for example Aeschbach-Hertig, 1994). The net flux F C of the isotope concentration C out of the water column is given by the combined effect of evaporation and precipitation: where C P is the concentration in precipitation, C L the surface concentration in the lake, E the evaporation rate, R the precipitation rate, h w the relative humidity with respect to surface water temperature, α the equilibrium fractionation coefficient and Δε the kinetic enrichment. For 3 H, we assume an equilibrium fractionation factor α of 0.89 and negligible kinetic fractionation (Aeschbach-Hertig, 1994). For 2 H, α is computed according to Majoube (1971), and the kinetic enrichment factor Δε is derived from relative humidity (Gonfiantini, 1986). Finally, we need to determine the precipitation isotope signals. For 3 H, we use the database of the Global Network of Isotopes in Precipitation (GNIP). Because data is scarce in Africa and there is no measurement station close to Lake Kivu, we use the average of surrounding stations to derive a time series with yearly data from 1953 to 2018 (IAEA/WMO, 2020). Note that there is no data available after 2000 and thus, we calculate the values from 2000 to 2018 using exponential extrapolation of the trend from 1990 to 2000 (see Figs. S1 and S2 in the supporting information). We express 2 H in δ notation, i.e. as excess with regard to its ratio with the Vienna Standard Mean Ocean Water (VSMOW), and we use a constant precipitation value of − 4.5‰ because: i) the value is in the range of observed values (IAEA/WMO, 2020), ii) data is too scarce for deriving seasonal variation, and iii) the value leads to a good agreement between modelled δ 2 H signals and observations at the lake surface. Note that we also use these 3 H and δ 2 H precipitation values for the surface inflows (tributaries) into Lake Kivu, thus neglecting evaporation during transport to the lake. The assumption of neglecting evaporation effects on the isotopic composition of the tributaries is motivated by the fact that the watershed is small and thus rainwater ends up in the lake quickly. In reality, there is substantial evaporation despite the small watershed and this in turn would lead to higher δ 2 H values in tributaries compared to precipitation. We neglected this evaporation effect in our study (i.e. assumed the same δ 2 H in both tributaries and precipitation) because we do not simulate 2 H dynamics in the surface waters.
The exchange of noble gases across the lake-atmosphere interface was approximated by prescribing air saturated water at 25 • C at the lake surface due to virtually constant concentrations of noble gases in the Fig. 3. Regression between N 2 and apparent temperature and salinity diffusivities (data from Sommer et al., 2019). The corrected regression line includes a correction factor, but preserves the slope. Molecular diffusion is used as the lower boundary of apparent diffusivity. a) Apparent temperature diffusivity K T,app , b) apparent salinity diffusivity K S,app .
atmosphere and because the surface water temperature of Lake Kivu is usually close to 25 • C.

Meteorological forcing
Simstrat requires air temperature (T a ), horizontal wind speed at 10 m above the water surface (UV 10 ), wind direction (φ), atmospheric vapor pressure (e a ), incoming short-wave radiation (SW in ), and incoming long wave radiation (LW in ). These meteorological variables are monitored continuously every 30 min by an automatic weather station on a platform on Lake Kivu (see Rooney et al., 2018 for details). W. Thiery from Vrije Universiteit Brussel kindly provided continuous data from October 2012 to September 2019 from this weather station. To allow for the simulation of longer timescales, we use the output of a regional climate model, which is used for scenario simulations within the ISIMIP project (Inter-Sectoral Impact Model Intercomparison Project, see Frieler et al. (2017) for scenario descriptions). The climate model output consists of daily modelled values of all important meteorological variables from the year 1660-2300 at Lake Kivu (grid cell 1.75 • S/29.2 • E). We chose MIROC5 from the four available regional models but the choice of climate model is not important, as long as interannual variability is well represented. The air pressure and air temperature of the climate model are corrected by adjusting their average to the measured average by the weather station (Fig. 2), while wind speed and incoming radiation are adjusted using Simstrat's built-in correction factors (Table 1). Note that the difference in air pressure and temperature (Table 1) is caused by the higher elevation of the climate model grid cell compared to Lake Kivu. We also noticed that wind directions differ between observations and the regional climate model, probably due to local effects below the grid scale of the climate model. However, simulation results are only weakly sensitive to wind direction in the lake model because wind direction mostly influences the timing and less the magnitude of mixing events. In addition, the magnitude of mixing events is scaled by the wind and seiche energy fit parameters of the model to match the observed temperature and salinity profiles. Finally, the difference in wind speed could be shown to be caused by the different time resolutions of the weather station (half-hourly) and the climate model (daily averages).
Using the corrected climate modelling output of MIROC5, we created two forcing time series: i) the control dataset "piControl" from 1660 to 2300 which represents climate forcing with constant, pre-industrial greenhouse gas concentrations, and ii) a combination of a historical dataset from 1860 to 2005 and a projection from 2006 to 2100, with a net radiative forcing of 6.0 W m − 2 in 2100, which will be called "RCP6" (representative concentration pathway 6.0). The goal of the "piControl" time series is to provide a realistic, but stochastic distribution of seasonal mixing depth for the steady-state and long-term simulations (see Results section), while RCP6 is only used for the transient simulation of 3 H from 1953 to 2018.

Water balance
The water balance of Lake Kivu can be written as with precipitation P, surface inflow I, groundwater inflow G, evaporation E and outflow O. Precipitation and evaporation are almost equal in Lake Kivu (Schmid and Wüest, 2012;Muvundja et al., 2014), meaning that the discharge of inflowing tributaries and groundwater is mostly balanced by the outflowing Ruzizi River. In our model, we use the values by Schmid and Wüest (2012), who suggest average values of 63 m 3 s − 1 for surface inflow, 41 m 3 s − 1 for groundwater inflow and 95 m 3 s − 1 for the outflow. We add the net difference between P and E to the outflow, which then amounts to 108 m 3 s − 1 .

Tritium ( 3 H) measurements
3 H (or T) is a radioactive hydrogen isotope, which is transported in the environment as tritiated water (HTO). It is extremely rare in nature due to low natural production and a relatively fast decay with a half-life of ~12.3 years. However, the atomic bomb tests in the 1950s and 60s released large amounts of 3 H into the atmosphere (up to several hundred tritium units TU, where 1 TU is 1 3 H atom per 10 18 water atoms). Since then, 3 H has been widely used in lake and groundwater environments as a physical tracer (Kipfer et al., 2002).
We took water samples using the pressure-proof sampling equipment described in Bärenbold et al. (2020b), and analyzed them for 3 H content by measuring its decay product 3 He according to Beyerle et al. (2000). In this procedure, samples are degassed under vacuum in order to evacuate the naturally present 3 He, then left on the shelf to allow for sufficient decay of 3 H to 3 He (typically 6-12 months), and subsequently analyzed for newly produced 3 He. Knowing the amount of 3 He produced in the sample and the elapsed time, one can back-calculate the original amount of 3 H present in the sample. Due to the high helium content in Lake Kivu (Bärenbold et al., 2020b), some of our samples had to be degassed twice to remove all naturally present 3 He. The estimated standard deviation of our measurements due to calibration uncertainty, and possible contamination with atmospheric air is in the range of 0.03-0.10 TU. In addition to this data, we added two unpublished 3 H measurements taken by our own group during an expedition in 2004.

Seasonal mixing dynamics (0-60 m)
Seasonal mixing in Lake Kivu usually occurs during the dry season and affects the top ~30-60 m. It is strongest during July and August due to high wind speeds and lower air temperatures, which lead to convective instability and enhanced mixing (Schmid and Wüest, 2012). Interannual variability of seasonal mixing depth is therefore mainly driven by meteorological conditions. The extent of the surface layer (where conductivity is nearly constant) gives an impression of varying mixing depth in different years (Fig. 4a). The typical range of surface temperatures in Lake Kivu is ~23-26 • C (Thiery et al., 2014;Morana et al., 2015) with warmer temperatures in the stratified rainy season and cooler temperatures in the windier dry season.
In the model, average mixing depth and surface temperature are mainly determined by the fit parameters of wind speed, short wave radiation and long wave radiation. By manually tuning these three parameters (Table S2 for numerical values), we can adequately reproduce mixing depth, which varies between ~20 and 60 m (Fig. 4b). Modelled water temperatures exhibit a range of ~23-27 • C (Fig. 4c), with peak water temperatures during stratified periods occasionally 1-2 • C higher than typical observations (Thiery et al., 2014;Morana et al., 2015). However, as long as the temporal distribution of the mixed layer extent Table 1 Averages of meteorological data from an automatic weather station on Lake Kivu (Rooney et al., 2018) and the output of the regional climate model MIROC5. All variables of the weather station are measured half-hourly and the daily averaged wind speed is shown for comparison (in parentheses). Wind speed is actually measured at a height of 7.2 m, and thus was corrected to 10 m assuming a logarithmic wind profile. is realistic, these higher surface temperatures will not influence our model results.

Turbulent mixing at intermediate depths (60-120 m)
Diffusive transport between ~60 and 120 m is governed neither by seasonal mixing nor by double-diffusive convection. This depth range is characterized by a strong increase in salinity with depth and therefore a strong density gradient. It was therefore put forward by Schmid et al. (2005) that turbulent mixing in this depth range is weak. In contrast, the stable isotope measurements by Ross et al. (2015a) indicate that mixing might be stronger than previously thought. In fact, their measurements  seem to show that the evaporative isotope signal at the lake surface is mixed down against the groundwater-induced upwelling.
Our simulation results suggest that direct wind-induced turbulent mixing is indeed very weak in Lake Kivu below 60 m. The wind-induced currents are not able to penetrate below the surface layer, and without considering the energy transfer via internal seiches, the modelled salinity gradient below the surface mixed layer becomes much steeper than observed (Fig. 5a). If the seiche energy parameter α s is tuned to 0.0027 (i.e. 0.27% of wind energy is transformed to seiche energy), our model can reproduce the observed salinity gradient very well (Fig. 5b), except for a small jump at ~120 m. This jump is an artefact caused by the transition from the k-ε mixing scheme to the parameterization of double diffusion. We were able to reduce the size of the jump by tuning the calibration parameter q NN but not to eliminate it. The parameter q NN determines how the mixing energy from internal seiching is distributed vertically as a function of stability N 2 . We conclude that internal seiching can explain why there is still a certain amount of turbulent mixing (10 − 7 -10 − 6 m 2 s − 1 ) in a depth region where the density gradient is very strong.

Steady-state simulation
While diffusive transport is computed using the k-ε closure (0-120 m) and a parameterization of double-diffusive convection (120-485 m), vertical advection is a function of the discharge, inflow depth and properties (T, S, CO 2 , CH 4 ) of the different groundwater springs present in Lake Kivu. Unfortunately, the properties of the groundwater sources are inherently difficult to measure in the field. However, despite the lack of data, we can use indirect evidence to shed light on the influence of the groundwater sources on Lake Kivu.
Using the water balance of Lake Kivu, Schmid and Wüest (2012) suggested that the total groundwater input into Lake Kivu is between 32 and 48 m 3 s − 1 . In addition, Ross et al. (2015a) were able to identify several warm and saline signals below 250 m, as well as a strong cool and less saline signal right at the chemocline (~250-255 m) and a more diffusive cool signal between ~90 and 195 m. Based on this information and the current vertical structure of the lake, we assume that: i) four hydrothermal point sources enter the lake below 250 m and plunge due to their high density, ii) two cooler point sources enter the lake at 253 and 190 m, and iii) one diffusive groundwater source, which has the same properties as the point source at 190 m, enters the lake between 90 and 195 m depth (Fig. 2 in Ross et al., 2015a).
We summarized this information in Table 2. Throughout this work, we refer to the ground-water sources at 190 and 253 m as cool sources 1 and 2, respectively; and the sources below 250 m as hydrothermal sources.
The water introduced as point sources is free to plunge or rise according to its density in Simstrat, thereby entraining ambient lake water. The final stratification depth of these groundwater inflows is not only influenced by their temperature, salinity and gas content, but also by the current lake density profile, which leads to a highly complex and interdependent calibration problem with many unknown variables. In order to calibrate inflow depth and discharge, but also temperature, salinity and CO 2 content of the groundwater sources, we chose to tune our model with the goal of producing a steady-state with temperature, salinity and gas profiles close to today's observations (Fig. 2). This approach is motivated by the findings of Bärenbold et al. (2020a) who conclude that CH 4 and CO 2 concentrations did not show a clear temporal trend during the last ~45 years. Conversely, it has been shown (Katsev et al., 2014;Sommer et al., 2019) that Lake Kivu experienced a warming trend on the order of 0.01 • C/yr between 2004 and 2015 throughout the stratified deep water below ~ 60 m. We chose to neglect this warming rate in our simulations for two reasons: i) the observed warming trend seems variable in time (i.e. no warming observed below 250 m from 1974 to 2004) and is comparably slow, and ii) we focus on the present state and the past evolution of Lake Kivu in this work, not on its future evolution. The causes for the recent lake warming have not been clearly identified. The warming near the surface is likely due to climate change and thus probably a recent phenomenon (Katsev et al., 2014). In contrast to this recent trend, the near-surface lake temperatures could have varied to some extent in both directions during the simulated period of the past 2000 years depending on past climate conditions. The causes for the recent warming of the deep water are unknown, and it needs to be further evaluated in the future whether they represent a persistent change in the thermal dynamics of the lake that requires adaptation of the model for future projections. Note, that a warming trend could be introduced into the model by slightly modifying the temperature of the hydrothermal groundwater inflows or the geothermal heat flux.
While most of the CO 2 stems from groundwater inflow this is not true for CH 4 , which was shown to contain a mixture of modern and old, 14 Cdead carbon (Schoell et al., 1988;Pasche et al., 2011). Using 14 C analysis, Pasche et al. (2011) showed that the degradation of organic matter can only provide a maximum of ~50% of CH 4 below 250 m. The remaining 50% stem either from reduction of old, geogenic CO 2 or from the direct inflow of geogenic CH 4 . The absolute CH 4 fluxes derived in Pasche et al. (2011) lead to an increase in CH 4 concentrations of around 7-8% within 30 years. However, such an increase is clearly not validated by the most recent measurements (Bärenbold et al., 2020a), which indicate approximately constant concentrations since the measurements of Tietze in 1974 (Tietze, 1978). Therefore, we assume that the actual CH 4 production below 250 m is close to a steady-state production of ~28 gC m − 2 y − 1 . According to Pasche et al. (2011), ~50% of this production is supplied by degradation of organic matter, and therefore we set F sed, CH4 in AED2 to 3.2 mmol m − 2 d − 1 CH 4 (14.0 gC m − 2 y − 1 ) for the whole lake depth. For the remaining 50% below 250 m, we will explore two scenarios: i) the reduction of geogenic CO 2 to CH 4 within the lake (implemented analogously to CH 4 production from organic matter), and ii) the direct groundwater-mediated inflow of geogenic CH 4 . In both cases, an additional ~3.2 mmol m − 2 d − 1 of CH 4 is added below 250 m, thus yielding a total CH 4 production of ~6.4 mmol m − 2 d − 1 (28 gC m − 2 y − 1 ) below 250 m. Fig. 6 shows simulation results for temperature, salinity, CO 2 and CH 4 for a simulation duration of 500 years (exact groundwater properties used in the simulation are provided in Table S4). The initial conditions used in this simulation are the temperature and salinity background profile of Ross et al. (2015a) and the CH 4 and CO 2 concentrations measured in 2018 (Bärenbold et al., 2020a). Our model Table 2 Summary of groundwater inflows into Lake Kivu used in this work. The inflow depths and discharges are calibrated in order to produce a steady-state close to today's profiles. assumptions and tuned inflows lead to a steady-state after 400-500 years for temperature, salinity and CO 2 concentrations, with CO 2 concentrations agreeing with observations well within measurement uncertainty. The agreement between data and simulation is also excellent for salinity with a root mean square error (RMSE) of 0.014‰, but not for temperature. While the general structure of the temperature profile is generally well reproduced (RMSE of 0.024 • C), modelled temperature is consistently ~0.3 • C too low between 250 and 350 m and ~0.1 • C too high below 400 m. Despite a large calibration effort on the groundwater properties, these differences could not be eliminated. The main difference between temperature and the other main agents in Lake Kivu is that the transport of salinity and gases is dominated by vertical advection (i. e. the upwelling caused by groundwater inflow), while for temperature both advection and diffusive transport are important (Sommer et al., 2019). As a consequence, it seems that the properties of the groundwater sources are well calibrated, while the parameterization of heat transport through the double-diffusive staircases only works satisfactorily. Indeed, Fig. 3a shows that the proposed linear regression reflects well the general relationship between stability and apparent temperature diffusivity. However, the variability is rather large, thus indicating that the prediction capacity of stability to deduce double-diffusive transport in Lake Kivu is limited.
For CH 4 , we explored two scenarios with different origins for the 14 Cdead part of CH 4 . Scenario 1 implies that the remaining 50% of CH 4 are produced by CO 2 reduction within the lake (Fig. 6d), while in scenario 2 the remaining 50% stem from groundwater (Fig. 6e). Scenario 1 indicates that CH 4 concentrations should reach a clear maximum between 300 and 400 m and decrease towards the lake bottom due to the inflow of CH 4 -free groundwater there. This decrease is supported by the CH 4 data of Tietze in 1974(Tietze, 1978 and the Eawag dataset in Bärenbold et al. (2020a), but not by the data of Tochon in 2003 (Schmid et al., 2005) and the UFZ dataset in Bärenbold et al. (2020a). Scenario 2 also shows slightly decreasing concentrations at the lake bottom, but this decrease is small, and could probably be eliminated by further tuning the CH 4 concentrations in the groundwater sources (see Table S4 for values used in scenario 2). Scenario 2 agrees better with observed CH 4 concentrations and thus, our results suggest that at least part of the remaining 50% of CH 4 are introduced by the hydrothermal groundwater sources. In both scenarios, CH 4 production from the decomposition of organic matter is assumed constant with depth. In reality, sediment focusing or the reduction of degradable organic matter at greater depths due to longer settling times could potentially lead to varying CH 4 production with depth. Fig. 6. Depth profiles at different simulation times of a 500 year steady-state simulation compared to observations. The temperature and salinity simulations are compared to the background data by Ross et al. (2015a) and the CO 2 and CH 4 simulations to Bärenbold et al. (2020a) (average of Eawag and UFZ datasets therein, and the uncertainty is the maximum uncertainty of both datasets). The simulation is initialized with the observations, and forced by constant groundwater inflows (see Table S4) and meteorological input (piControl). a) Temperature profile; the horizontal bars represent the depths where groundwater inflows stratify, b) salinity profile, c) CO 2 profile, d) CH 4 profile with 50% of CH 4 produced by each degradation of organic matter and CO 2 reduction, e) CH 4 profile with 50% of CH 4 produced by each degradation of organic matter and groundwater inflow.

Water tracer simulations
In this section, we use the calibrated steady-state model to reproduce the 3 H and δ 2 H measurements from this work and by Ross et al. (2015a) respectively.

3 H simulations: water age and origin
There are no direct 3 H measurements available from the groundwater sources. In general, 3 H in groundwater can have two origins, namely a low, rather constant background signal due to interaction with rocks (Mamyrin and Tolstikhin, 1984), and a variable signal due to gas exchange with the atmosphere. Therefore, we make the following assumptions for our model runs: i) the old, hydrothermal groundwater in the Lake Kivu basin has a constant background 3 H signal of 0.5 TU due to interaction with rocks, ii) young groundwater stems from the recent infiltration of precipitation and has an atmospheric 3 H signal and, iii) all groundwater sources in Lake Kivu can be a mixture of old and young groundwater. Note that we assume that the time from infiltration of precipitation to groundwater recharge is short, i.e. there is no lag that would decrease the 3 H concentration due to radioactive decay. This assumption is justified by the fact that the northern shore, where the groundwater sources were found (Ross et al., 2015a) consists of permeable volcanic rock, which allows fast infiltration of rainwater. Fig. 7a shows simulation results under the assumption that every groundwater source is entirely composed of old hydrothermal groundwater compared to field observations from 2004 to 2018 (see Table S5 for numerical values). The observations show a decrease in time in the seasonal mixed layer from 2.2 TU in 2004 to 1.6 TU in 2018 and a decrease with depth from 1.6 TU to close to 0 TU at 200 m. The decrease in time can be explained by decreasing 3 H concentrations in the atmosphere due to radioactive decay of the remaining bomb 3 H, while the decrease with depth shows that the water below 60 m is predominantly old (i.e. no recent interaction with the atmosphere), and that vertical mixing is weak. The assumption that all groundwater sources are entirely composed of old groundwater agrees very well with the observations in the mixed layer, as well as with the observation below 300 m. However, between 100 and 200 m, the measured 3 H concentrations are higher than predicted by the simulation, thus indicating that the cool groundwater sources do also contain atmospheric 3 H. Indeed, if we adjust the fraction of young groundwater in the cool sources 1 and 2 to 45% and 20%, respectively, we get a very good agreement with the observations in the corresponding depth range (Fig. 7b).
The finding that the cool groundwater sources are partially young is supported by the absence of surface tributaries in the north of Lake Kivu where Ross et al. (2015a) identified the cool sources. This northern part of the catchment has a surface area of ~686 km 2 , which accounts for ~14% of the whole Lake Kivu catchment area of 4940 km 2 . We can estimate a lower bound for the amount of water which infiltrates and feeds the groundwater sources by using the surface inflow of 63 m 3 s − 1 (Schmid and Wüest, 2012). As 86% of the catchment produce ~63 m 3 s − 1 of run-off, the remaining 14% would yield a potential groundwater recharge of ~10 m 3 s − 1 . However, precipitation is expected to infiltrate quickly through the porous volcanic rock in the north and thus, there is less water lost to evaporation than in the remaining watershed. An upper bound of the infiltrated water can be derived by using the average precipitation of 1470 mm yr − 1 (Muvundja et al., 2014) which would give ~32 m 3 s − 1 over volcanic part of the catchment. The amount of young water required by our calibrated model to fit the observed 3 H profile is ~15.5 m 3 s − 1 with a range of ~11.5-18 m 3 s − 1 (range estimated from uncertainty of 3 H) and is therefore well within the range of 10-32 m 3 s − 1 given by run-off and precipitation data.
We conclude that the cooler groundwater sources are probably a  Table S6 for numerical values), c) simulation results with δ 2 H concentrations in cool sources 1 and 2 as measured by Ross et al. (2015a), d) simulation results with δ 2 H concentrations in cool sources 1 and 2 adjusted to fit the background observations of Ross et al. (2015a), (see Table S6 for numerical values).
mixture of old and more recent groundwater, with ~45 and ~20% of recently infiltrated water in cool sources 1 and 2, respectively. Conversely, our only observation below 300 m suggests that the deep groundwater sources predominantly consist of old, hydrothermal water.

δ 2 H simulations contradict field observations
In contrast to 3 H, there are field observations available for δ 2 H for at least some of the groundwater sources (Ross et al., 2015a). Below 200 m, Ross et al. (2015a) found very variable δ 2 H groundwater signals ranging from − 2.8 to 17.4‰. In our model, we cannot use these values because our implemented hydrothermal sources do not directly correspond to the hydrothermal sources found by Ross et al. (2015a). Therefore, the δ 2 H concentrations of the hydrothermal sources are adjusted to reach the best agreement between simulation and observed background observations by Ross et al. (2015a). In contrast to the warm sources, our cool sources 1 and 2 were clearly identified by Ross et al. (2015a) and thus, Fig. 7c shows a simulation where the δ 2 H concentrations of cool groundwater sources 1 and 2 correspond to the values found by Ross et al. (2015a), while in 7d we adjust all sources so that the model fits the observed background profile of Ross et al. (2015a). In both Figures, we also show the δ 2 H observations of Tietze (1978) and Katsev et al. (2014). Comparing the three datasets illustrates the seasonal variability of the δ 2 H values in the surface layer, as well as the generally good agreement below 100 m. Due to better data consistency (i.e. less scattering of measurements at similar depths), we use the data of Ross et al. (2015a) for the calibration of δ 2 H values in the inflows (Fig. 7d). To reproduce the observed δ 2 H values in the surface layer, the δ 2 H concentration in precipitation was fixed to − 4.5‰, meaning that the model cannot reproduce seasonal or interannual variations at the surface.
When using the δ 2 H observations of Ross et al. (2015a) as groundwater input, the observed δ 2 H values above 250 m are considerably higher compared to simulated δ 2 H (Fig. 7c). The δ 2 H values in groundwater necessary to provide a good agreement of observations and simulation are around 15 and 30% higher for cool sources 1 and 2, respectively (Fig. 7d, see Table S6 for numerical values). These differences are well beyond measurement uncertainty and in addition, both δ 2 H groundwater observations consist of several measurements with rather low variability (pers. comm. K. A. Ross). The results of the δ 2 H simulations are further corroborated by δ 18 O simulations (not shown in this work), which exhibit the same systematic deviation from observations.
We conclude that our model is able to reproduce the observed δ 2 H profiles within the uncertainty of almost all data points of Ross et al. (2015a) if the inflow concentrations are freely adjustable and allowed to have an evaporative δ 2 H signal. In this case, the observed stable isotope profiles are not in contradiction with low turbulent mixing in Lake Kivu. However, the observed δ 2 H concentrations close to cool sources 1 and 2 are not high enough to explain the observed lake profiles. Potentially, a certain discrepancy could be because δ 2 H in the groundwater inflows (or the groundwater inflows themselves) are variable in time. However, it is unlikely that this could explain a difference of 15 and 30% in inflow concentrations, and we are unable to fully understand the reason for this disagreement between model and observations.

Origin and maintenance of Lake Kivu's density stratification
It is generally assumed that Lake Kivu has experienced several (partial) degassing events in the past (Haberyan and Hecky, 1987;Zhang et al., 2014;Ross et al., 2015b;Votava et al., 2017;Uveges et al., 2020). This may have happened by spontaneous ebullition due to gas oversaturation in the deep waters. Another possibility is the creation of instability within the water column due to either subaquatic volcanism (Ross et al., 2015b) or hyperpycnal flows due to exceptional rain events (Zhang et al., 2014). Ross et al. (2015b) and Uveges et al. (2020) both estimate that the last large mixing event may have occurred ~800-1000 years ago, with several smaller events between 800 years BP and the present. In this context, we can use our model to evaluate how long it takes Lake Kivu to reach its current state. Because we don't know when and to what extent Lake Kivu has mixed in the past, we assume that i) the lake is completely mixed and gas-free at the start of the simulation, ii) the groundwater discharge and properties (i.e. temperature, salinity and gas concentrations) are constant in time and the same as today, and iii) the average climatic conditions are constant and the same as today.
Our simulated salinity and gas concentrations agree well with today's state after roughly 1500 years (Fig. 8), thus confirming that the density-driven inflow of several groundwater sources can explain today's vertical lake structure. This result is in contrast to arguments by Hirslund (2020) who suggested that the inflow of groundwater sources cannot be responsible for Lake Kivu's density stratification. For temperature, it takes more time to reach a steady-state. The reason is that double-diffusive transport is much faster for temperature compared to salinity and gases due to the much higher molecular diffusion coefficient. While the vertical transport of salinity and gas concentrations is dominated by advection, double diffusion contributes around half of the vertical heat flux. The double-diffusive heat flux is parameterized as a function of N 2 and therefore, a stable temperature profile will only develop once the density stratification (mainly determined by salinity and CO 2 ) is fairly stable after ~1500 years. At this point, the overshooting temperature below 350 m decreases again and converges to ~26 • C.
Our model results suggest that we should see slightly decreasing deep water temperatures if the lake is close to a steady-state. In contrast, recent temperature observations indicate an increase of ~0.01 • C yr − 1 (Sommer et al., 2019). There are several possible answers to this problem: i) the lake is not at a steady-state and thus, temperature, salinity and gas concentrations are still slightly increasing, ii) the lake was close to a steady-state, but there was an increased discharge of hydrothermal water during the last decades which counteracts the temperature decrease, and iii) our parameterization of double-diffusive heat flux is not suitable for the simulation of a near-homogeneous lake at the start of the simulation and thus the simulated overshooting is a model artefact. Note that answer ii) would also imply a slight steepening of the main chemocline at 250 m due to increased upwelling as it was observed recently (Hirslund, 2012). However, this steepening stopped again after 2012 (pers. comm. M. Schmid).
Assuming constant groundwater inflows throughout the simulation duration, our simulations show that it is improbable that today's lake structure has developed within only 1000 years from an initially homogeneous and degassed lake. The time which is necessary to reach a steady-state close to today's observed temperature, salinity and gas concentrations rather seems to be around 2000 years. Therefore, a potential mixing event ~1000 years ago, as suggested by Ross et al. (2015b) and Uveges et al. (2020), could not have mixed the entire lake. However, a partial mixing event ~1000 years ago or even smaller mixing events during the last few centuries are well compatible with our simulations.

Atmospheric noble gas depletions in the deep water
In an earlier publication, we described that the atmospheric noble gases (i.e. the noble gases whose main source is the atmosphere) Ne, 36 Ar and Kr are depleted by 55-70% in the deep water of Lake Kivu (Bärenbold et al., 2020b). Therein, we discussed three different mechanisms to explain the observed noble gas depletions, namely i) continuous outgassing, ii) the inflow of noble gas depleted groundwater, and iii) a relic of a past large outgassing event. Mechanism i) was subsequently excluded because the isotope ratios 20 Ne/ 22 Ne and 40 Ar/ 36 Ar did not show the expected depletion patterns (Bärenbold et al., 2020b). We further argued that mechanism ii), the inflow of depleted groundwater is the most probable reason for the observed profiles.
In this section, we use the model setup of the previous section to test and evaluate hypotheses ii) and iii). For hypothesis ii), we adjusted the Fig. 8. Depth profiles at different simulation times for a long-term transient simulation compared to observations. The temperature and salinity simulations are compared to the background data of Ross et al. (2015a) and the CO 2 and CH 4 simulations to Bärenbold et al. (2020a) (average of Eawag and UFZ datasets therein, and the uncertainty is the maximum uncertainty of both datasets). The simulation is initialized with a completely homogeneous and gas-free lake, and is forced by constant groundwater inflow (see Table S4 for groundwater properties) and meteorological input (piControl). a) Temperature profile; the horizontal bars represent the depths where groundwater inflows stratify at 1000 and 2000 years respectively, b) salinity profile, c) CO 2 profile, d) CH 4 profile with 50% of CH 4 produced by each degradation of organic matter and CO 2 reduction, e) CH 4 profile with 50% of CH 4 produced by each degradation of organic matter and groundwater inflow. Fig. 9. Simulated Neon depth profiles over 2000 years, starting with a homogeneous, degassed lake, compared to observations from Bärenbold et al. (2020b): a) depletion of hydrothermal inflows is 40-50% and depletion of cool sources 1 and 2 is 5 and 10% compared to air saturated water (ASW), respectively (see Table S6 for numerical values), b) all inflows contain air saturated water at 25 • C.
concentrations in the groundwater sources to fit the observed noble gas profiles in Bärenbold et al. (2020a). Conversely, in hypothesis iii) the noble gas concentrations in the groundwater are fixed to those of air saturated water (ASW) at 25 • C. The resulting Ne profiles indicate that we can exclude that the currently observed profiles are a result of a past outgassing event only (Fig. 9b). In contrast, the observed concentrations can be perfectly reproduced by the inflow of noble gas depleted groundwater (Fig. 9a). The best agreement with observations is reached if the hydrothermal sources are depleted by 40-50% (Table S6) and the cool sources 1 and 2 by 5% and 10%, respectively. The depth profiles of the other atmospheric noble gases measured in Bärenbold et al. (2020b) can be successfully reproduced as well, but are not shown here.

Model limitations
Compared to past one-dimensional modelling attempts on Lake Kivu (Schmid et al., 2005), one of the main advantages of our model is the fact that turbulent transport is dynamically calculated throughout the whole water column. From 0 to 120 m, we use Simstrat's built-in k-ε turbulence closure, and below, we apply a newly derived parameterization, which is based on measurements by Sommer et al. (2013). This situation creates an artificial transition point at 120 m. However, while both the transition at 120 m and the parameterization are valid for today's lake structure back to at least 1976 (Newman, 1976), we have no information about the presence and extent of double-diffusive staircases during the past centuries. Particularly, during periods with a weaker density stratification, we do not know how well our parameterization performs. Nevertheless, it can be assumed that an initially weak density stratification quickly becomes stronger due to the constant deep inflow of salts and CO 2 . As a consequence, although turbulent transport might initially be underestimated by our parameterization, this probably does not significantly alter the simulations results.
Another important simplifying assumption of our modelling exercise is the fact that climatic conditions, as well as the discharges and properties of all groundwater sources, are assumed to be constant in time. Part of the cooler groundwater sources seems to have an atmospheric 3 H signature, and thus we expect the discharge and properties of these sources to vary with climate conditions (i.e. higher discharge in a wet climate, lower discharge during droughts). Furthermore, volcanic activity can greatly influence hydrological parameters like evaporation and rock permeability and thus lead to variable infiltration and groundwater properties over time. For example, Votava et al. (2017) suggest that the intermittent deposition of CaCO 3 found in lake sediments indicate a corresponding intermittency of deep groundwater inflow into Lake Kivu. However, for our long-term simulations it is not crucial to know whether climate or inflow do vary over shorter timescales due to the long water residence times in the stably stratified deep waters.
Finally, it is clear that our model is overparameterized, which is unavoidable given the lack of direct observations of the properties of the groundwater sources. Besides the calibration parameters of Simstrat-AED2 (including meteorological and mixing fit parameters), the properties of every groundwater source (temperature, salinity, gas content, stable isotopes etc.) as well as their inflow depths and discharges can be freely adjusted. Thus, even if the simulations are constrained by observed vertical profiles and the water balance, likely several parameter sets exist which yield comparable results. For example, the discharge of the hydrothermal groundwater sources could be lower, but its temperature, salinity and gas content could be higher to keep equal the heat, salt and gas input by the sources. However, such changes to groundwater properties also influence density and thus the depth where the source stratifies, which makes the search for alternative solutions difficult. Our suggested set of parameters and inflow properties might not be unique, but it strongly suggests that to form and maintain the current state, several warm and dense groundwater sources below 260 m and two cooler sources at ~190 and 250 m depth are required.

Conclusions
We developed and calibrated a one-dimensional model for Lake Kivu, exploiting the fact that the lake is well mixed horizontally below ~60 m. In order to simulate the effect of dissolved gases CO 2 and CH 4 on lake physics, we first coupled the physical lake model Simstrat with the biogeochemical model library AED2. We then developed a parameterization of vertical transport of heat, salts and gases across the doublediffusive staircases of Lake Kivu (below 120 m) as a function of lake stability N 2 . After manually calibrating temperature, salinity and gas content of the inflowing groundwater sources, we were able to reproduce today's lake profile by a steady-state simulation. Using this steadystate model and our 3 H dataset allowed to show that the cooler groundwater sources, which enter the lake between 100 and 250 m, are probably a mixture of modern and old, hydrothermal water.
Most importantly, we showed that the evolution of Lake Kivu from a completely homogeneous and gas-free state to today's density stratification within 2000 years is successfully reproduced using a onedimensional lake model with constant groundwater inflows. Our results demonstrate that the stratification depth of groundwater inflows markedly shifts as a function of the evolving density structure of the lake. As the duration of 2000 years is an upper bound, which is valid for an initially completely mixed lake, our simulations indicate that a partial mixing event ~1000 years BP (Zhang et al., 2014;Ross et al., 2015b;Uveges et al., 2020) is in good agreement with the currently observed state of the lake. Furthermore, our simulation results suggest that the depletion of atmospheric noble gases in the deep water of Lake Kivu (Bärenbold et al., 2020b) are not a relic from a past outgassing event, but have to be sustained by an active process (i.e. inflow of noble gas depleted groundwater).
In this work, we used our model to focus on finding a steady-state solution for today's lake structure, and on reproducing the past evolution of Lake Kivu. In the future, our model can be a valuable tool to simulate the effect of a changing climate on lake temperature, stratification and nutrient availability. In addition, it can facilitate the decisionmaking concerning the construction and operation of ongoing and future methane harvesting projects in Lake Kivu. In order to better constrain the discharge and properties of the groundwater inflows in future model applications, more research on the hydrogeological setting in the Lake Kivu region would be beneficial. Moreover, the sampling of inflowing groundwater sources using a remotely operated underwater vehicle (ROV) could allow the determination of temperature, salinity and gas concentrations of groundwater sources without contamination with ambient lake water.
The Kivu-Simstrat model is derived from the coupling of the physical lake model Simstrat and the biogeochemical library AED2. This coupled model (without the specific changes for Lake Kivu) represents a computationally efficient and physically sound tool to simulate geochemical tracers and biochemical processes in lakes and reservoirs worldwide in the future.

Funding
This work was supported by the Swiss National Science Foundation [grant number 200021_160114].

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. organization of the field work and Reto Britt, Michael Plüss, Maximilian Schmidt and the whole team from LKMP for help with tritium sampling in Lake Kivu. We also want to thank Alexandra Lightfoot and Matthias Brennwald from the noble gas laboratory at ETH Zurich for help with tritium sample analysis. Adrien Gaudard and Davide Vanzo provided helpful insights into modelling and the implementation of new Simstrat features. Particularly, Adrien's manual was of great help. We also would like to acknowledge Wim Thiery from Vrije Universiteit Brussel who made the data of the Lake Kivu weather station available to us.