On the risk of a dissolved gas-triggered limnic eruption in Lake Kivu

Abstract

Lake Kivu is distinguished by several unique characteristics that set it apart from other lakes around the world. One of the notable features is a temperature increase with depth, accompanied by unusual staircase-like patterns in the thermodynamic and environmental parameters. The lake also experiences suppressed vertical mixing due to stable density stratification, with its deep water separated from the surface water by chemoclines. Additionally, Lake Kivu contains high concentrations of dissolved methane (CH₄) and carbon dioxide (CO₂), and there is no standard method for measuring their concentrations. The lake is also recognized as a renewable energy source due to its continuous supply of CH₄, and it demonstrates a quadruple-diffusive convection transport mechanism. These factors contribute to the lake's distinctiveness. The occurrence of catastrophic limnic eruptions at Lakes Nyos and Monoun, along with the structural similarities between these lakes and Lake Kivu, raises serious concerns about the likelihood of a similar disaster in Lake Kivu in the future. The scale of threats posed in Lake Kivu can be orders of magnitude greater than the other two lakes, given its 3000 times larger size, two to four orders of magnitude higher content of dissolved CO₂, containing substantial quantities of CH₄ in addition to CO₂ in solution, and holding a far denser population living in its much wider catchment area. The present study aims to assess the probability of a future gas outburst in this giant lake by numerical modeling of its hydrodynamics over the next half a millennium. The turbulent transport is calculated using the extended k–ε model. An implicit Euler method is applied to solve the governing partial differential equations on a vertically staggered grid system, discretized using a finite-volume approach. Since the previously calibrated model successfully reproduces the measured lake profiles, the same tuned parameter values are used in this study, assuming a stable steady-state condition in the future. The results of our simulations effectively address common concerns regarding the risk of a gas burst in the lake due to buoyancy instability-triggered overturn and/or supersaturation of the water column.

---

Environmental significance

The occurrence of catastrophic limnic eruptions in Lakes Nyos and Monoun, along with the structural similarities between these lakes and Lake Kivu, raises serious concerns about the probability of a similar disaster in this lake in the future. We assess the probability of a future gas outburst in this giant lake by numerical modeling of its hydrodynamics over the next half a millennium. Given that the Lake Kivu catchment area is densely populated with around 5.7 million residents, the findings of this study hold significant implications for the region's environmental sustainability and societal well-being.

1 Introduction

Lake Kivu is situated on the western side of the East African Rift Valley, at the border between the Democratic Republic of Congo (DRC) and the Republic of Rwanda, near the equator, between 1°34′–2°30′S and 28°50′–29°23′E. This lake and several other lakes filled large local depressions formed by subsidence in the Rift Valley because of continental plate divergence. Being of volcanic origin and located at an elevation of 1463 m above sea level, the lake has a maximum depth of about 485 m, an average depth of 240 m, a surface area of around 2370 km², and a volume of approximately 550 km³.^1–4 This lake is one of the deepest freshwater reservoirs, being ranked the world's thirteenth and twentieth deepest by mean and maximum depths, respectively.

Lake Kivu lies in a highly active tectonic and volcanic region, which has shaped the current structures in the lake and its surrounding areas. Nyamuragira and Nyiragongo are the most active volcanic mountains bordering the north of the lake. The lake's depth was only about 85 m in the past due to northward drainage of lake water into Lake Edward. The formation of the Nyamuragira and Nyiragongo mountains about 20 000 years ago produced lava flows that created a natural dam, blocking the outflow of the former shallow lake toward Lake Edward. This process led to the formation of Lake Kivu, as the water level gradually rose over approximately 4000 years to fill the valley. After rising to almost 122 m above the current lake level, the water body found a new drainage path towards the south about 10 000 years ago, directing the outflow of Lake Kivu towards the north of Lake Tanganyika in Central Africa through the 117 km-long Ruzizi River.^5,6 Therefore, Lake Kivu is a relatively young and consequently unstable aquatic system from a geological point of view, at least compared to other nearby lakes. For example, the latest eruption of Nyiragongo volcano and the subsequent surface fracturing and associated seismic events occurred on May 22, 2021.⁷ Therefore, the presence of many geologically recent features in the active volcano-tectonic areas surrounding the lake implies that similar volcanic and seismic events are likely to occur in this dynamic system in the future.

The complex shape of Lake Kivu was subdivided into five basins with different physicochemical features, distinct characteristics, and fluid inputs, as shown in Fig. 1. These basins include a large one, which is called Main Basin, and four smaller ones, which are called Kabuno Bay, Kalehe, Ishungu, and Bukavu. Among these basins, the Main Basin—the largest and deepest part of the lake—is the most important. Kabuno Bay, with a maximum depth of about 105 m, exhibits the most distinct limnological characteristics among the five basins, as it is almost completely isolated from the rest of the lake by a few islands.^8,9 The presence of basins with distinct features constituting the lake can be related to the combined effects of different processes such as biological cycles, rock–water interactions, and deep subaquatic inflows.^3,10–15 It is worth noting that, in addition to these natural factors, the water chemistry in this lake has also been influenced by human-related activities such as untreated domestic sewage discharge from nearby communities, outflows of wastewater treatment plants, runoff of chemical fertilizers consumed by farmers, and effluents from various industries and factories in the region.^16,17 Considering the fact that the catchment area of Lake Kivu, with a population of 5.7 million, is densely inhabited, the crucial impact of human emissions on the lake's composition is further underscored.

Fig. 1 Map of Lake Kivu and its five basins. The white lines in this map indicate the boundaries of the neighboring countries of the Democratic Republic of Congo, Rwanda, Uganda, Burundi, and Tanzania. Reprinted (adapted) from ref. 18. Copyright 2023 Elsevier Ltd.¹⁸

Lake Kivu possesses several unique characteristics that distinguish it from other lakes worldwide. First, the temperature of the lake's water column is controlled by both the heat exchange of the lake surface with the atmosphere and that of the water body with subsurface groundwater inflows. These two heat exchange sources increase temperature with depth, a rare phenomenon not commonly observed in other lakes. Second, because of the inflow of warmer, more saline, and gas-rich subaquatic groundwater streams, thermodynamic and environmental parameters such as salinity, temperature, concentrations of nutrients (i.e., ammonium, phosphate, and silica), and dissolved gases (i.e., CH₄ and CO₂) are strongly and stably stratified, increasing with depth. This means that the profiles of the mentioned parameters exhibit peculiar staircase structures, composed of a series of very short depth intervals with drastic gradients sandwiched between fairly homogeneous, low-gradient thick layers. Third, the lake exhibits an unusual vertical density profile, driven by the joint effects of stratified temperature, salinity, dissolved gas gradients, and groundwater inflows. Indeed, the so-called staircase stratification consists of several alternating well-mixed regions with relatively uniform density over a long distance separated by thin interfaces with very steep density gradients. Unlike other lakes, vertical mixing in Lake Kivu has been suppressed by stable density stratification, which inhibits turbulent transport. Although density stratification is common in natural water bodies, it is typically disrupted by strong wind events or seasonal cooling that drive vertical mixing. In contrast, Lake Kivu is exceptional because its deep waters remain persistently isolated by strong stratification, which effectively constrains vertical exchange. Fourth, another large-scale distinctive feature of Lake Kivu is the effective decoupling of the deep water body from the open surface layer by the main density gradients, commonly referred to as chemoclines. Occurring at depths of approximately 85, 191, 257, 313, and 390 m, the chemoclines play a crucial role in regulating vertical flow within the stratified ecosystem and shaping the present configuration of the lake. Indeed, the presence of chemoclines has significantly diminished the vertical mass fluxes of dissolved gases by blocking emission pathways into the atmosphere, accumulating enormous volumes of CH₄ and CO₂ over time. Fifth, the high CH₄ and CO₂ concentrations relative to other lakes and oceans make this lake particularly special among aquatic systems. The first estimation of the total dissolved gas volumes in the lake, which is still used in the evaluations, was performed by Tietze in 1978.² According to Tietze's calculations, the amounts of CH₄ and CO₂ stored in the lake over hundreds of years are approximately 60 km³ and 300 km³, respectively, at a standard temperature of 0 °C and a standard pressure of 1 atm. Sixth, due to its high dissolved gas concentrations compared to other lakes, there is currently no standard procedure for accurately measuring them in Lake Kivu. Contrary to other CH₄ reservoirs around the globe, Lake Kivu can be recognized as a renewable source because CH₄ is continuously supplied to the lake by the conversion of dead organic matter to CH₄ and CO₂ by methanogens, the reduction of CO₂ by H₂, and subaquatic springs carrying CH₄. These features place Lake Kivu among the most extraordinary lakes on Earth and make it an exceptional aquatic ecosystem for studies by limnologists, ecologists, climatologists, volcanologists, geologists, geochemists, and microbiologists.

In addition to the unique limnological, geochemical, physical, and architectural characteristics mentioned above, Lake Kivu is also a special system from a fluid-dynamics perspective. In this lake, four factors compete to determine the water density profile and thus affect density stratification: temperature, salinity, and CH₄ and CO₂ concentrations. Although all these parameters increase with depth, their contributions to the evolution of the staircase structure differ (Fig. 2). While increases in temperature and CH₄ concentration with depth decrease water density and destabilize the water column, salinity and CO₂ concentration enhance stratification stability by increasing density with depth. Since the stabilizing effects are approximately six times larger than the destabilizing ones, the system maintains its overall gravitational stability. The nonlinear dynamics arising from the competition between the stabilizing and destabilizing factors induce double-diffusive (thermohaline) convection processes, consisting of coexisting transport mechanisms of double diffusion through the sharp interfaces and convection phenomenon in the mixed layers. Newman¹⁹ was the first researcher to describe the role of this mixing process in Lake Kivu and reported the presence of 150 vertically-stacked mixed layers. Double-diffusive convection occurs in systems with gradients of two properties, such as temperature and salinity, where molecular diffusivities differ markedly by orders of magnitude (i.e., two or more) and have opposite effects on the vertical density profile.^20–23 While this crucial transport mechanism is commonly observed in natural settings such as oceans and saline lakes,²⁰ Lake Kivu stands out as uniquely distinct for fluid dynamic studies due to the contributions of four, rather than the usual two, density-relevant components in shaping its staircase structure. Therefore, the commonly used term “double diffusion” or “thermohaline”, which refers to the effects of two diffusion components only, needs to be revised to “quadruple diffusion” in the context of Lake Kivu; this change properly reflects the interplay of CH₄ and CO₂ concentrations, along with dissolved salts and heat, in shaping the lake's stratification dynamics and transport processes.

Fig. 2 Vertical profiles of temperature, salinity, and dissolved CH₄ and CO₂ concentrations in Lake Kivu. Reprinted (adapted) from ref. 24. Copyright 2020 The Authors.²⁴

From the standpoint of CH₄ concentration and vertical density stratification, Lake Kivu can be subdivided into four relatively homogeneous zones, roughly interfaced by the chemoclines. These major zones, from the surface to the bottom of the lake, are called the Biozone, the Intermediate Zone, the Potential Resource Zone, and the Resource Zone. The Biozone, also known as the mixolimnion or oxygenated layer, hosts all aerobic organisms, including fish species and zooplankton. It is supplied directly with oxygen from the atmosphere and from algae for various biological activities. This zone, extending from the surface to a depth of about 60–65 m, contains negligible dissolved gas concentrations, as almost all CH₄ that reaches this zone is consumed by bacteria under oxic conditions. The Intermediate Zone, located within the approximate depth range of 60–200 m below the lake surface, experiences continuous dilution from groundwater streams and upwelling water. This flushing process results in a low concentration of CH₄ in this layer. With an average CH₄ concentration of about 5.7 mol m⁻³, the Potential Resource Zone serves as a transition from the overlying diluted layers to the underlying gas-rich deep water. As the name implies, this zone has the potential to become a source of CH₄ for local communities if its dissolved CH₄ content gradually increases to a level that is commercially harvestable, or if advancements in harvesting technologies significantly enhance the economic viability of CH₄ extraction. Starting from a depth of approximately 260 m and ending at the bottom of the lake, the Resource Zone, with an average CH₄ concentration of about 16.7 mol m⁻³, contains most of the total and almost all the exploitable solution CH₄-in-place. It is worth noting that the dissolved CH₄ concentration in this zone is, on average, approximately three times as large as that in the overlying Potential Resource Zone.

The zonation pattern of CH₄ is also observed for CO₂, with the sharpest concentration increase occurring at a depth of approximately 257 m, coinciding with the onset of the Resource Zone. Similar to CH₄, the CO₂ concentration in the Biozone is almost zero, since there is no chemocline in this zone to impede gas escape to the atmosphere. The average CO₂ concentrations in the Intermediate Zone, Potential Resource Zone, and Resource Zone are around 10.6, 26, and 81.9 mol m⁻³, respectively. As shown in Fig. 2 and 3, the substantially higher CO₂ concentrations in the Resource Zone compared to shallower layers have resulted from a combined influence of subaquatic groundwater inflows and the lake's meromictic stratification. Indeed, the warm, saline, and CO₂-rich groundwater streams enter the lake at depth, reinforcing the chemoclines that suppress vertical mixing and inhibit the release of dissolved gases to the atmosphere. These conditions have allowed CO₂ to accumulate over long timescales, particularly in the deep water body. As discussed in subsequent sections, these processes have led to a pronounced enrichment of CO₂ in the Resource Zone. The similar concentration profiles of CH₄ and CO₂ depicted in Fig. 2 and 3 reflect the distinct yet coupled dynamics of the two gases.

Fig. 3 Comparison between dissolved CH₄ and CO₂ concentrations measured by the 2018 campaign and those measured in 1974, 2003, and 2004. Reprinted (adapted) from ref. 24. Copyright 2020 The Authors.²⁴

Despite its unique characteristics, Lake Kivu shares some features with other lakes, two of the most infamous of which are Lakes Nyos and Monoun in Cameroon. Although these two lakes are much shallower and smaller than Lake Kivu, they are similar in having density stratification and containing large amounts of dissolved gases. For this reason, lessons from their history can be applied to the conservation, strategic planning, and sustainable management of the Lake Kivu ecosystem. In separate events during the mid-1980s, Lakes Nyos and Monoun experienced gas bursts that suddenly released massive amounts of dissolved CO₂ into the atmosphere, killing 1746 and 37 people, respectively, along with numerous wild animals and cattle.²⁵ The lethal effects of these natural catastrophes were mainly due to inhalation, leading to asphyxia. At high concentrations, CO₂ is rather toxic, causing unconsciousness and death at levels exceeding 100 000 ppm. These rare, tragic eruption incidents, which had not been reported in the literature previously, were attributed to the collapse of the density stratification structures in the water columns. To date, two scientific hypotheses have been proposed to explain the causes of the disastrous overturns in these lakes.^25–27 The first one, which attributes the disruption of the stratified structures to explosive processes such as volcanic eruptions at the bottom of the lakes, is inconsistent with the absence of geochemical signatures, such as elevated levels of sulfur and chlorine, typically expected from volcanic emissions. On the other hand, a more likely and well-supported hypothesis is the overturn of water columns, exsolution of CO₂ due to its displacement from the deep waters to shallower depths with lower hydrostatic pressure, release and coalescence of gas bubbles, and, consequently, gas bursts at the surface. Although rockfalls and a landslide were speculated to have triggered the overturns in Lakes Nyos and Monoun, respectively, other potential triggering mechanisms include seismic activity, the inflow of cold-water streams, and the generation of internal waves in the lakes.

Although gas-charged lakes are rare, they are not exclusive to Central Africa. Globally, only a few lakes are known to accumulate hazardous concentrations of dissolved gases under stable stratification, with Lakes Kivu, Nyos, and Monoun representing the most prominent cases. Other volcanic crater lakes, such as Lake Quilotoa in Ecuador and Lake Albano in Italy, also exhibit elevated CO₂ levels, but their scale and hazard potential are far less significant than those of the three African lakes. Thus, the occurrence of gas-rich, eruption-prone lakes is geographically concentrated in the East African Rift and Cameroon Volcanic Line, making these systems both exceptional and of considerable scientific interest.

The exceptionally high concentrations of dissolved gases in Lakes Kivu, Nyos, and Monoun result from a combination of geological and biological processes. Since these lakes are located in volcanically active regions, magmatic gases, particularly CO₂, seep upward through hydrothermal groundwater inflows. In Lake Kivu, CH₄ has been produced both geogenically and biogenically, the latter driven by microbial degradation of organic matter under anoxic conditions in the deep water. These inputs and in situ gas generation, compounded with strong and persistent density stratification maintained by saline groundwater, have inhibited seasonal mixing and allowed progressive accumulation of dissolved gases over centuries.

The occurrence of the catastrophic limnic eruptions at Lakes Nyos and Monoun, and the structural similarities between these lakes and Lake Kivu, raise serious concerns about the likelihood of a similar future disaster in Lake Kivu. It is worthwhile to note that the scale of threats posed by a gas burst of Lake Kivu can be orders of magnitude vaster than the Cameroon lakes in view of its 3000 times larger size,²⁸ two to four orders of magnitude higher content of dissolved CO₂,²⁵ containing substantial quantities of CH₄ in addition to CO₂ in solution, and holding a far denser population living in its much wider catchment area. Therefore, assessing the potential risk of a limnic eruption in this giant lake and implementing appropriate preventative measures is of the utmost importance. The present study aims to analyze changes in dissolved gas concentrations, the stability of the lake's water column, and the probability of a future gas outburst through numerical modeling of its hydrodynamics over the next half-millennium.

2 Debates on the potential risk of a limnic eruption

Almost all previous assessments of the danger of a gas-driven eruption in Lake Kivu were based on the measurements of gas concentrations and their changes over time. The success of researchers who attempted to monitor gas concentrations in the lake has depended on the design and accuracy of the sampling equipment and the measurement techniques used. The first measurements of gas concentrations in this lake date back to 1935, when Damas²⁹ documented CO₂ and H₂S levels. Despite his remarkable work in gathering comprehensive data from the lake, the recorded concentrations are not sufficiently reliable since Damas degassed more than half of the dissolved CO₂ into the atmosphere during the measurements. Schmitz and Kufferath³⁰ were the first researchers to measure CH₄ concentrations in addition to CO₂ by releasing dissolved gases under atmospheric conditions. Although this measurement technique yielded reasonably accurate results for CH₄, it underestimated CO₂ concentrations because the solubility of CO₂ in the water samples was higher than that of CH₄, retaining appreciable amounts of CO₂ in solution at low pressures. Degens et al.^1,31 also applied a similar method for the concentrations of dissolved CH₄, CO₂, H₂S, and N₂. However, their measurements face the challenges that were discussed above. The first comprehensive and reliable measurements of CH₄ and CO₂ concentrations in the lake were carried out by Tietze,² who measured both the exsolved gas and the portion that remained in solution in water samples using a special-purpose apparatus.

The most controversial datasets that raised serious concerns about the potential threat of a gas burst in Lake Kivu were those measured by Halbwachs and Tochon, and published by Schmid et al.³² Based on their analysis of the datasets, CH₄ and CO₂ concentrations in the deep water, respectively, increased by about 15–20% and 10% over 30 years from the 1974 measurements by Tietze;² conjecturally, the increases in gas concentrations led them to speculate that the concentrations would reach saturation levels in the present century if the observed trends continue in the coming years, triggering a severe gas eruption. Despite that a later detailed analysis of the carbon cycle by Pasche et al.³³ in 2011 revealed that CH₄ concentrations are increasing at a slower rate than that hypothesized before, the notionally-critical safety warning of Schmid et al.³² in 2005 was taken seriously by the governments of the countries sharing the lake such that an intercalibration campaign was initiated by the Lake Kivu Monitoring Program (LKMP) through Energy Development Corporation Limited (EDCL) as a subsidiary of the government-owned Rwanda Energy Group (REG) with the aims of measuring CH₄ and CO₂ concentrations in the lake using different accurate methods, estimating an uncertainty range for measured concentrations and the total solution gas-in-place, and predicting and monitoring the temporal evolution of gas concentrations in the future. The research teams that participated in the campaign include the French National Centre for Scientific Research (CNRS), the Swiss Federal Institute of Aquatic Science and Technology (Eawag), KivuWatt Ltd., and the Helmholtz Centre for Environmental Research (UFZ).

The 2018 campaign shed new light on the evolution of gas concentrations in Lake Kivu. The participants obtained vertical concentration profiles using different measurement methods at approximately the same time to compare the results of various methods for determining the same quantity. The results of all measured datasets for CH₄ and CO₂ concentrations, along with the previous ones, are presented in Fig. 3. By comparing the curves in this figure, it can be seen that the concentration data obtained over a long period of time using various methods agree fairly well, falling within the measurement uncertainties of 10% and 5% for CH₄ and CO₂ concentrations, respectively.²⁴ In other words, since the differences between the measured data of gas concentrations over the past half-century are within the expected inaccuracy limits of the applied methods, no firm conclusion can be drawn about the temporal changes of these parameters in the lake. Therefore, the potential danger of a spontaneous gas outburst triggered by supersaturation of CH₄ and CO₂ in the lake cannot be unambiguously evaluated solely based on the available field data.

The use of numerical modeling to predict the temporal evolution of density stratification, stability, and dissolved gas concentrations in Lake Kivu is limited to the study by Schmid et al.,³² which used the AquaSim software,³⁴ a tool for data analysis and simulation of aquatic systems. Although their model incorporates most of the key transport processes governing the lake dynamics, it suffers from several inherent assumptions and simplifications. These limitations can lead to misinterpretations, particularly in long-term simulations. However, they also offer opportunities for refinement to improve simulation accuracy and ensure robust interpretations. The most significant of these assumptions include parameterizing transport through the staircases using apparent thermal diffusivity values, which underestimate the vertical heat flux by factors of 2–20,³⁵ performing a static calculation of turbulent transport for a fixed density stratification, and not stratifying the subaquatic springs over depth intervals as a function of water density. Furthermore, the unknown parameters of their one-dimensional model were determined by fitting the profiles observed by Tietze² and Schmid et al.,³² assuming that CH₄ and CO₂ concentrations steadily increased over 30 years. This assumption did not fully account for the above-mentioned uncertainties associated with different measurement methods. Consequently, the simulation results from their model, which simplifies certain physical conditions, indicated that the observed increases in concentrations will continue in the coming decades, potentially leading to a gas eruption triggered by supersaturation. Thus, both the available measured data and the AquaSim model have critical limitations in predicting the future evolution of the lake's dissolved gas concentrations and stability.

3 Problem statement and research question

As explained above, the majority of previous studies on the risk of a gas outburst in Lake Kivu have focused on comparative analyses of field-measured data to identify temporal trends. A comparison of the data reported by different research teams participating in the intercalibration campaign with the previous measurements demonstrates no clear evidence of substantial changes in the dissolved gas concentrations in the lake over the past half-century. In fact, all measured data for CH₄ and CO₂ concentrations lie well within the uncertainty limits associated with differences in the design and operation of the various measurement devices and the underlying physics. Therefore, the evolution of gas concentrations and the potential threat of a supersaturation-triggered gas outburst cannot be reliably understood from available field data. On the other hand, the simplified model of Schmid et al.³² may be less suitable for long-term simulations spanning several human generations and may not fully capture the lake's real-world setting, as it was calibrated by matching simulation results to the seemingly increasing observed data. These limitations motivated us to revisit the temporal evolution of the density stratification, stability, dissolved gas concentrations, and the risk of a limnic eruption in Lake Kivu using a representative model that incorporates modifications to address the aforementioned shortcomings and accounts for the uncertainties associated with the measured data.

4 The lake model

To address the identified challenges and ensure a more realistic simulation of the lake's hydrodynamics, the modeling approach used by Bärenbold^36,37 was adopted in this study. This model, which was initially developed by Goudsmit et al.³⁸ at Eawag and later modified by Bärenbold, is called Simstrat. Over the past two decades, Simstrat has been updated several times, incorporating various key features into the original model.^39–41 Recently, Bärenbold^36,37 coupled the Simstrat model with the open-source biogeochemical library Aquatic EcoDynamics (AED2) to enable explicit simulation of the carbon cycle, and called the resulting model “Simstrat-AED2”. Then, the complexity of double-diffusive convection processes and the extraordinarily high concentrations of CH₄ and CO₂ in Lake Kivu prompted them to implement a range of modifications to Simstrat-AED2 and to name the resulting model “Kivu-Simstrat”. This model offers several advantages over the AquaSim model,³² including the application of the k–ε turbulence model for dynamic computation of transport at depths shallower than 120 m, parameterization of double diffusive transport through the staircases within the depth range of 120–485 m using real field measurements of salinity and temperature profiles, and allowing subaquatic inflows to plunge or rise depending on their density relative to that of the surrounding water body. All these improvements convinced us to choose Kivu-Simstrat as an adequate and representative numerical model for simulating the hydrodynamic processes in Lake Kivu. Given the size and complexity of this model, it cannot be presented in full detail here. Therefore, it is only briefly discussed in this section, and interested readers are referred to the relevant literature^36,37 for further details.

4.1 Dimensions of the model

Except for Kabuno Bay, vertical profiles recorded at different locations in other bays in Lake Kivu confirm that the distributions of various parameters, such as gas concentrations, are nearly homogeneous in the horizontal plane throughout the lake,⁴² indicating much faster horizontal mixing compared to vertical transport processes. Given the negligible horizontal gradients, a one-dimensional model is deemed sufficient to capture transport mechanisms in the system under consideration. Accordingly, this study emphasizes vertical processes in the lake, while noting that higher-dimensional approaches (i.e., 2D or 3D) may be more appropriate for sub-basins with distinct morphologies, hydrodynamic characteristics, biogeochemical features, and subaquatic inflows. However, developing such comprehensive models would require extensive datasets from all lake basins, which are currently limited.

4.2 Governing differential equations

As the extended k–ε model, which is the basis for the calculations of mixing and vertical density profile in Kivu-Simstrat, has been described elsewhere,⁴³ it is thus not presented here. With the vertical z-axis positive upwards, the one-dimensional partial differential equations governing the distributions of temperature, horizontal velocity components, turbulent kinetic energy, and turbulent dissipation rate are, respectively, as follows:³⁸where T is the temperature in °C, t is the time in s, z is the depth in m, A is the areal extent of the lake at depth z, ν′ and are, respectively, the molecular and turbulent thermal diffusivities in m² s⁻¹, c_P and ρ₀, respectively, denote the reference specific heat and density of the lake water in J kg⁻¹ K⁻¹ and kg m⁻³, H_sol stands for the shortwave solar radiation received by the lake in W m⁻², H_geo represents the heat flux from geothermal sources in W m⁻², u and v are the mean velocity components in the horizontal directions in m s⁻¹, ν and ν_t are, respectively, the molecular and turbulent viscosities in m² s⁻¹, f is the Coriolis parameter in s⁻¹, k is the turbulent kinetic energy in J kg⁻¹, ν_k and ν_ε, respectively, represent turbulent diffusivities of kinetic energy and dissipation rate in m² s⁻¹, P is the shear stress production in W kg⁻¹, P_Seiche denotes the turbulent kinetic energy produced by internal seiching in W kg⁻¹, B is the buoyancy production in W kg⁻¹, ε is the turbulent dissipation rate in W kg⁻¹, and c_ε1, c_ε2, and c_ε3 are the dimensionless constants of the model.

To fully describe the lake's dynamics, the above equations must be solved simultaneously with three fundamental conservation equations: the conservation of mass (continuity) equation, the conservation of momentum (Navier–Stokes) equation, and the solute and heat advection–diffusion equations. These three kinds of differential equations are, respectively, written as follows:

where ρ is the water density in kg m⁻³, w is the velocity component in the vertical direction in m s⁻¹, p is the pressure in Pa, g is the acceleration due to gravity in m s⁻², C_i is the concentration of the solute i in mol m⁻³, D_i is the mass diffusivity of the solute i in m² s⁻¹, λ is the thermal conductivity in W m⁻¹ K⁻¹, and other parameters are as defined previously.

4.3 Boundary conditions

As described previously, there are several major hydrological inflows and outflows to/from the lake that must be imposed as boundary conditions. To efficiently solve the system of linear equations resulting from the discretization using a tridiagonal matrix solver, a no-flux condition was applied to the boundaries, and fluxes across the boundaries were implicitly treated as sources and sinks, rather than being explicitly incorporated in the modeling.

4.4 Grid size and geometry

The lake's bathymetry, defined as the areal extent as a function of underwater depth, determines the size and geometry of the model's grid cells. As presented in Fig. 4, the plot of cross-sectional area variation with depth in the Kivu-Simstrat model was created from previously published bathymetric data.^4,5 High-resolution bathymetric survey and seismic data reveal that the lake is characterized by a steep-walled depression with a deep central basin.⁵ For hydrodynamic modeling, this morphology can be idealized as a bowl-shaped structure, capturing the essential features of steep margins enclosing a deep basin. As schematically illustrated in Fig. 5, this simplification provides a tractable geometric representation of the otherwise complex lake morphology. Considering this, the optimal approach for its numerical simulation is to discretize it vertically. To achieve this, the entire lake's structure was modeled as a series of staggered trapezoids, with the areas of the top and bottom bases derived from Fig. 4, and a height of 2 m. The volume of the segment at a depth of z in this discretization scheme is calculated using the trapezoidal rule:where V_z is the volume of the segment, h_z is its height, and A_z and A_z+1 are the areas of the bases.

Fig. 4 The plot of cross-sectional area variation with depth in the Kivu-Simstrat model.

Fig. 5 Schematic of the staggered grid discretization of the lake medium used in the present study.

4.5 Numerical scheme

An implicit Euler method was used in the Kivu-Simstrat model to solve the above-mentioned partial differential equations on a vertically staggered grid system discretized using a finite-volume approach. Fig. 5 presents a schematic of the staggered-grid discretization of the lake medium used in the model. In this type of discretization scheme, the water body was vertically subdivided into n_v non-equidistant volumes from the bottom of the lake to the surface, and they were indexed from 1 to n_v. The indices of the bottom and top faces of the i^th volume are i and i + 1, respectively. While quantities such as u, v, and T, which represent the mean values within the vertical intervals, are located at the centers of the volume elements, the turbulence model parameters k and ε are calculated at the interfaces between the volumes. The choice of time discretization in this model is not essential due to the fully implicit treatment of the problem; therefore, a constant time step (i.e., 300 s) was used to guarantee high numerical stability and low computational cost.

4.6 Numerical resolution requirements for long-term simulations

To ensure that the governing equations remain valid over the large spatial scale of Lake Kivu and the multi-century duration of our simulations, particular attention was given to the choice of numerical resolution. Following Bärenbold,^36,37 we adopted a vertical resolution of 2 m, which adequately reproduces the smoothed background profiles of temperature, salinity, and dissolved gas concentrations while remaining computationally efficient for long-term simulations extending over centuries. Although finer depth resolutions with grid cell sizes less than 1 cm were tested in an attempt to better capture the sharp gradients between thin mixed layers, these simulations were both computationally prohibitive and inconsistent with the heat flux data, resulting in unrealistically steep gradients at the interfaces between mixed layers.^36,37 Consequently, we used a vertical resolution of 2 m for all long-term simulations in this study, as it strikes an effective balance between physical accuracy and computational feasibility. Equally important, the chosen resolution is consistent with the meter-scale depth intervals of the observed data used for model calibration and validation, as temperature, salinity, and dissolved gas concentrations in Lake Kivu have routinely been measured at vertical spacings on the order of meters. Thus, the numerical discretization matches the effective resolution of the available field data, ensuring that the model neither under-resolves observed gradients nor introduces artificial detail unsupported by measurements. This alignment between model resolution and sampling interval further supports the adequacy of the 2 m grid size for reproducing the large-scale stratification and long-term dynamics considered in this study.

4.7 Water density calculation

The effects of salinity and dissolved gas concentrations on the water density of volume elements were accounted for by applying the following linear corrections:⁴⁴

ρ(T,S,C_CO2,C_CH4) = ρ(T)(1 + β_SS + β_CO2C_CO2 + β_CH4C_CH4) (11)

where ρ is the density in g cm⁻³, T is the temperature in °C, S is the water salinity in g kg⁻¹, and C_CO2 and C_CH4 are, respectively, the concentrations of CO₂ and CH₄ in g kg⁻¹. In this equation, β_S, β_CO2, and β_CH4 are the correction coefficients due to the contributions of dissolved salts, CO₂, and CH₄, respectively, with values of 0.75 × 10⁻³, 0.284 × 10⁻³, and −1.25 × 10⁻³ kg g⁻¹. The temperature-dependent water density is approximated by the following correlation:⁴⁵

ρ(T) = 0.9998395 + 6.7914 × 10⁻⁵T − 9.0894 × 10⁻⁶T² + 1.0171 × 10⁻⁷T³ − 1.2846 × 10⁻⁹T⁴ + 1.1592 × 10⁻¹¹T⁵ − 5.0125 × 10⁻¹⁴T⁶ (12)

where ρ(T) is the lake water density at the temperature T in g cm⁻³. Although other dissolved gases are present in the lake in addition to CH₄ and CO₂, they were excluded from the water density calculation owing to their negligible concentrations, as discussed in the following sections.

4.8 Calibrated parameters

The unknown and uncertain parameters, especially the characteristics of subaquatic groundwater sources, biogeochemical properties in the AED2 model, and constants of the governing differential equations, were manually tuned by Bärenbold^36,37 salinity, temperature, and gas concentrations in the lake. These parameters include the depth, rate, salinity, temperature, and dissolved gas concentrations of inflowing groundwater streams, the influx of CH₄ from the lake-sediment interface, short- and long-wave radiations, wind speed above the lake surface, and seiche energy. Since this calibrated model successfully reproduces the measured profiles, the same calibrated parameters were employed in this study, assuming a stable steady-state condition in the future. Interested readers are referred to the aforementioned references for a detailed description of the history-matching procedure and the resulting values of the tuning parameters.

5 Simulation results and discussion

The calibrated model was run for 500 years, from 2024 to 2524, because this time scale is long enough to allow constituents to diffuse throughout the system and capture the full range of behavior of any slow processes (e.g., diffusion), given the lake's large size in all directions.

5.1 Long-term salinity evolution and lake stability

The lake salinity in the calibrated model was calculated from electrical conductivity and ionic composition data. In the prediction phase, the salinity profile at any simulation time results from a balance between the inputs and outputs through groundwater seepage and lake surface tributaries, respectively, natural vertical upwelling, and precipitation onto the lake surface. The simulated vertical salinity profiles after 100, 200, 300, 400, and 500 years are presented in Fig. 6. In this and all subsequent figures, the initial condition, when displayed, refers to the year 2024. This figure indicates that the lake salinity increases only slightly over half a millennium of simulation. This suggests that the lake dynamics will be very close to a steady-state condition in the future. In other words, the salinity profiles indicate that the system under study will not reach an exact steady state, but a quasi-steady state assumption can be considered reasonable. The negligible change in the vertical salinity profiles will favor maintaining lake stratification as eqn (11) mathematically describes the strong dependence of density on salinity. Another interesting point that can be drawn from Fig. 6 is that the total mass flow rate of salt entering the lake via the groundwater discharges nearly balances the rate of salt outflow by the Ruzizi River (i.e., approximately 3.88 × 10⁶ t yr⁻¹); otherwise, remarkable salinity increase trends in the lake would be expected as a result of salt accumulation over such a long period of time. A theoretical justification for this balance is the dominance of natural upwelling induced by the groundwater streams over diffusion, such that salts delivered to the lake by the subaquatic seepage are efficiently transported to the surface by upward advection. It is worth noting that the simulation results indicate negligible salt sedimentation in the lake. The almost-invariant distribution of lake salinity is also particularly important from the standpoint of ecological safety, as it will continue to allow local people who rely on salinity-sensitive marine-based industries such as fishing and agriculture to sustain their long-term livelihoods and economic well-being.

Fig. 6 Predicted future vertical salinity profiles in the lake.

5.2 Long-term temperature evolution and density structure in the lake

Besides salinity, temperature is one of the main factors affecting the staircase structure of the lake density. The future evolution of this parameter was also predicted using the numerical model and shown in Fig. 7. The calculated temperature profile at any given time was determined by an energy balance between the regional climate, the geothermal heat flux, the double-diffusive heat flux, and the groundwater inflows. The competition between these different types of heat sources and transfer mechanisms determines spatio-temporal temperature patterns in the lake. Therefore, the temperature at a given depth may increase, decrease, or even remain constant over time, depending on the relative fluxes of the mentioned heat sources and the values of thermal transport properties throughout the entire lake. This line of reasoning clearly explains why the future temperature profiles depicted in Fig. 7 remain almost unchanged in the Intermediate and Potential Resource Zones, whereas the simulated distributions in the Resource Zone exhibit a shifting trend toward lower values. On the other hand, the simulated profiles in the Biozone, which are mainly controlled by the interaction with the changing regional atmosphere and turbulent mixing, do not follow any regular pattern. It is also interesting to note that the decline in the evolution of temperature distributions in the Resource Zone is more pronounced over the next 300 years, after which it becomes almost time-invariant. The primary reason for the decrease in temperature with time at the bottom of the water column is the cooling effect of the subaquatic discharges that have somewhat lower temperatures than the ambient water in the Resource Zone; this is while the above-mentioned sources of heat have already reached a relatively stable balance in the overlying zones, such that a quasi-steady state was attained. From the safety margin viewpoint, the overall decrease in lake temperature over time will be favorable for producing denser deep water and, thus, greater buoyant stability, as water density is inversely related to temperature, according to eqn (12).

Fig. 7 Predicted future vertical temperature profiles in the lake.

5.3 Future dissolved CO₂ concentration profiles and lake stability

The simulation results for future CO₂ concentration profiles are presented in Fig. 8. Similar to the lake salinity profile, the current vertical distribution of CO₂ concentration is predicted not to evolve over the next half a millennium. This evidence further supports the hypothesis that the lake will have a quasi-steady state. The unchanged pattern in Fig. 8 is attributed to a balance between the inflow of CO₂-rich water through the subaquatic springs, its consumption via geogenic reduction to CH₄, its degassing into the atmosphere, and its outflow via the Ruzizi River. The microbial reduction of CO₂ to CH₄ prevents the accumulation of CO₂ in the lake; while this situation is favorable for keeping the system away from a CO₂ supersaturation state, it is unfavorable for the density stratification, which could potentially become stronger in the absence of CO₂ reduction. From a comparative perspective, the simulated future distributions of CO₂ concentration using Kivu-Simstrat are completely different from the significantly increasing ones previously predicted by the AquaSim model. As discussed in the previous sections, this discrepancy stems mainly from the calibration of the AquaSim model parameters to the apparently increasing measured data, which were actually almost constant over 30 years, given the uncertainties involved in different measurement methods for gas concentrations. For this reason, along with the aforementioned model improvements, the nearly invariant distribution of CO₂ concentration predicted in this study is a more reliable projection of future lake conditions. Therefore, based on the simulation results, the potential risk of a gas outburst triggered by CO₂ supersaturation in the coming centuries is very low.

Fig. 8 Predicted future vertical profiles of CO₂ concentration in the lake.

5.4 Future dissolved CH₄ concentration profiles and lake stability

Contrary to temporally constant distributions of CO₂ concentration, steadily increasing trends are predicted for future CH₄ profiles in the lake. The gradual buildup of CH₄ concentration, which is evident in Fig. 9, predominantly occurs below the main chemocline at a depth of about 257 m. In other words, the Resource Zone is the only region in the lake expected to experience a remarkable temporal increase in CH₄ concentration over the following centuries. This observation strongly underscores the critical role of chemoclines, especially the main one, in vertical transport through the stratified water column and in shaping the future configuration of the lake. Indeed, the main chemocline is so strong that it significantly reduces the vertical mass fluxes of dissolved gases by blocking emission pathways into the overlying zones. The suppression of the upward transport, combined with the supply of CH₄ by subaquatic sources, microbial methanogenesis from organic matter, and the reduction of CO₂, gives rise to the accumulation of a vast amount of dissolved CH₄ in the Resource Zone over time. The simulation results shown in Fig. 9 indicate that the gradual buildup of CH₄ increases its average concentration in the Resource Zone from the initial value of 16.7 mol m⁻³ in 2024 to about 19.5, 21.4, and 22.2 mol m⁻³ after 200, 400, and 500 years, respectively. Considering safety margins, the rising levels of CH₄ in the lake negatively affect buoyant stability. This is due to the inverse relationship between water density and CH₄ concentration, as described in eqn (11).

Fig. 9 Predicted future vertical profiles of CH₄ concentration in the lake.

5.5 Future trends of other dissolved gas concentrations in the lake

In addition to CH₄ and CO₂ as the main dissolved gases, the temporal variations of other dissolved gases such as N₂, O₂, He, Ne, Ar, and Kr were also examined to more accurately estimate their contributions to the saturation level of the water column. Since N₂ concentration has not yet been measured in the lake, its future distributions cannot be predicted using the model. For other dissolved gases, their simulated future profiles are shown in Fig. 10. According to the plots in this figure, while the deeper water is richer in He, the concentrations of other gases (i.e., O₂, Ne, Ar, and Kr) decrease with depth; these opposite trends can be attributed to the volcanic and atmospheric origins of He and other gases, respectively. One of the key features inferred from Fig. 10 is that the future profiles of He concentration exhibit shifting trends towards lower values. In contrast, the distributions of O₂, Ne, Ar, and Kr concentrations remain almost constant over time. Such different temporal behaviors stem from a much higher diffusion coefficient of He in water compared to other gases, resulting in an upward transport rate that exceeds its magmatic inflow flux, leading to gradual decreases in the future He profiles in the lake.

Fig. 10 Predicted future vertical profiles of O₂, He, Ne, Ar, and Kr concentrations in the lake.

5.6 Transport mechanisms

From a mechanistic perspective, heat and mass transport in Lake Kivu play a central role in shaping the temporal evolution of temperature and solute distributions. To better understand the relative contributions of various transport mechanisms in the lake, advective and diffusive flux distributions of heat and salt are presented in Fig. 11 and 12, respectively. Several important facts can be drawn from the analysis of this figure. First, advective heat and salt fluxes below a depth of 400 m are three orders of magnitude larger than those at shallower depths. This is mainly because of the bowl-shaped structure of the lake basin, which leads to much larger values of both heat and mass fluxes at the bottom of the lake than at the top, despite lower rates of deeper subaquatic inflows. Second, there is a depth interval of about 50 m from 257 to 307 m, where advective fluxes are negative. As depicted in Fig. 13, this depth range corresponds precisely to the interval with negative upwelling flow rates. This downward flow in the mentioned range can be attributed to the presence of the major chemocline at a depth of 257 m, which acts as a physical barrier, blocking the advective transport of heat and solutes from the underlying to the overlying water column. However, the chemocline does not eliminate vertical transport. Instead, buoyancy forces within this layer generate downward velocities (Fig. 13), resulting in negative rather than zero advective fluxes (Fig. 11 and 12). In other words, while the chemocline suppresses the upward exchange of heat and solutes, the velocity field continues to transport properties downward. This behavior can be interpreted as a restorative sinking effect, whereby displaced parcels are returned downward by density stratification, producing net downward advection. Therefore, the negative advective flux below the major chemocline reflects downward transport rather than the absence of transport. Although the main chemocline significantly reduces the advective pathway, it is interesting to note that it has little effect on diffusive fluxes. The principal chemocline's key role in interrupting the vertical continuity of the lake is evident in the distributions of temperature, salinity, and gas concentrations, as depicted in Fig. 6–9. Third, a comparison of the relative contributions of advection and diffusion mechanisms to heat and mass transfer in the lake indicates that advection predominantly governs the vertical transport of dissolved salts and gases. However, both advection and diffusion mechanisms are responsible for heat transfer at comparable levels. This is because the thermal diffusivity in water is, on average, two orders of magnitude greater than the diffusion coefficients of dissolved species. The buoyancy-induced heat and mass transfer occurs when the lake's surface loses buoyancy and becomes denser than the deep water. A brief examination of Fig. 14 suffices to conclude that buoyancy does not influence heat and mass transport in Lake Kivu, as there will be no loss of buoyancy over time.

Fig. 11 Predicted future vertical distributions of (a) advective and (b) diffusive heat fluxes.

Fig. 12 Predicted future vertical distributions of (a) advective and (b) diffusive salt fluxes.

Fig. 13 Predicted future vertical distributions of upwelling flow rate.

Fig. 14 Predicted future vertical profiles of water density in the lake.

6 Environmental implications for the lake stability and safety

The simulation results can also be utilized to analyze public safety in the vicinity of Lake Kivu. The safety margin of the lake can be assessed from the perspectives of buoyancy-driven instability and water column supersaturation. The stability of stratification, determined by the balance between stabilizing and destabilizing factors on the vertical density profile, is commonly quantified by the square of the buoyancy (Brunt‐Väisälä) frequency as follows:where N is the Brunt‐Väisälä frequency in s⁻¹, g is the acceleration due to gravity in m s⁻², ρ₀ is the density of fresh water, ρ(z) is the water density at any particular depth of the lake, and z is the depth increasing downwards from the lake surface in m. The positive, zero, and negative values of N² correspond to the water body's stable, neutral equilibrium, and unstable states, respectively. The stability N² profiles of the stratification at different times were calculated using the future density distributions and presented in Fig. 15. The positive N² profiles in this figure suggest that the water column in Lake Kivu is strongly density-stratified. The stability analysis also indicates that the strength of the stratification will remain approximately the same after 500 years of simulation. Therefore, concerns about the potential overturn of Lake Kivu's water column, arising from density-driven buoyancy instability, sudden release of dissolved gases, and, consequently, subsequent surface eruptions, can be reasonably mitigated.

Fig. 15 Predicted future vertical profiles of stability N² in the lake.

The buoyancy-driven lake overturn discussed above can be interpreted as a form of Rayleigh–Taylor (RT) instability, which occurs when a denser fluid overlies a lighter one in a gravitational field. In this regard, RT instability has long been recognized as a fundamental driver of turbulence and mixing. Therefore, it links environmental buoyancy instabilities to a broader class of hydrodynamic processes that also include Richtmyer–Meshkov and Kelvin–Helmholtz modes of instability. Moreover, RT instability has been extensively studied at various scales and different contexts,^46–49 providing a framework for interpreting buoyancy-triggered overturn. A detailed discussion of different instability forms and their coupled modes is beyond the scope of this study, and interested readers are referred to the relevant literature for further information.⁵⁰ Nevertheless, identifying the potential buoyancy-triggered overturn phenomenon in Lake Kivu as a form of RT instability underscores its broader scientific significance.

Although N² is a practical indicator for quantifying the stability of the lake's stratified structure, it does not provide any insight into the proximity of the water column to the saturation level. For this reason, the total dissolved gas pressure (TDGP) was calculated at each lake depth and compared with the corresponding hydrostatic pressure to determine the difference between TDGP and hydrostatic pressure, indicating the “safety margin”. The larger the difference, the wider the safety margin, and the lower the risk of an eruption. Indeed, dissolved gas bubbles begin to form and coalesce when the TDGP at a given depth exceeds the hydrostatic pressure. In the present work, the exerted TDGP was calculated for dissolved gases CH₄, CO₂, N₂, O₂, He, Ne, Ar, and Kr by applying the concepts of partial pressure, Henry's law, and solubility, as described in the SI. The simulated profiles of TDGP over the next 500 years are shown in Fig. 16, along with the hydrostatic pressure of the water column. As illustrated in this figure, the depths near the lake surface are the only regions in the column where the TDGP approaches the hydrostatic pressure. To evaluate and understand the contributions of different dissolved gases to the TDGP, the estimated partial pressures of these gases after 500 years of simulation were plotted together in Fig. 17. This figure indicates that the main reason for such a proximity in the shallow depths of the lake is the high concentrations of dissolved N₂ in water. Since N₂, rather than CH₄ and CO₂, is mainly responsible for the TDGP in the top part of the water body, it is of no practical importance. Fig. 16 also illustrates that the safety margin widens as we move from the surface to the bottom of the lake. In other words, the hydrostatic pressure in the Potential Resource and Resource Zones, which serve as the main storage spaces for CH₄ and CO₂, is too high to be exceeded by the TDGP. Therefore, the risk of sudden degassing due to the supersaturation of the water body in Lake Kivu is predicted to remain very low over the next half-millennium.

Fig. 16 Predicted profiles of TDGP and the hydrostatic pressure of the water column over time.

Fig. 17 Simulated contributions of different dissolved gases to the TDGP after 500 years of simulation.

7 Limitations and potential time-varying external factors

In this work, we assumed that external drivers, including atmospheric temperature, precipitation, inflow/outflow rates, tectonic activity, volcanic inputs, and anthropogenic interventions, remain constant over the half-millennium time span of the simulations. This assumption enabled us to isolate the intrinsic dynamics of the lake under present-day baseline conditions. Nevertheless, it is important to point out that long-term changes in these drivers could significantly alter lake stability. In particular, factors such as climate warming, shifts in precipitation and runoff, rising atmospheric CO₂, changes in wind regimes, tectonic or volcanic disturbances, water level fluctuations, and human-related activities (e.g., gas extraction and wastewater discharge) have the potential to affect the future evolution of dissolved gas concentrations and stratification stability in Lake Kivu. Although a comprehensive evaluation of this extensive parameter space with inherently uncertain future trajectories is beyond the scope of the present study, our results and conclusions apply under the assumption of stationary conditions or minor changes in external factors.

8 Conclusions

The following conclusions can be drawn from the simulation results and analyses performed in the present study:

• The lake salinity is predicted to increase only slightly over half a millennium of simulation. The invariant distribution of salinity means that the total mass flow rate of salt entering the lake via the groundwater discharges nearly balances the rate of salt outflow by the Ruzizi River because of the dominance of natural upwelling over the diffusion mechanisms. This negligible change in the vertical salinity profiles will favor maintaining stratification and, consequently, the ecological safety of the lake.

• The future temperature profiles remain almost unchanged in the Intermediate and Potential Resource Zones, whereas the simulated distributions in the Resource Zone exhibit shifting trends toward lower values. On the other hand, the simulated profiles in the Biozone, which are mainly controlled by the interaction with the changing regional atmosphere and turbulent mixing, do not follow any regular pattern. From the safety margin viewpoint, the overall decrease in lake temperature over time will be favorable for generating denser deep water, thereby enhancing buoyant stability.

• Similar to the lake salinity profile, the current vertical distribution of CO₂ concentration is predicted not to evolve over the next half a millennium. Therefore, based on the simulation results, the potential risk of a gas outburst triggered by CO₂ supersaturation in the coming centuries is very low.

• Contrary to temporally constant distributions of CO₂ concentration, steadily increasing trends are predicted for future CH₄ profiles in the lake. The gradual buildup of CH₄ concentration predominantly occurs below the main chemocline in the Resource Zone. The accumulation of a vast amount of dissolved CH₄ in the Resource Zone over time can be attributed to a combination of factors, including the suppression of the upward transport by the principal chemocline, the continuous supply of CH₄ by subaquatic sources, microbial methanogenesis from organic matter, and the reduction of CO₂ to CH₄.

• While the deeper water is richer in He, the concentrations of other gases (i.e., O₂, Ne, Ar, and Kr) decrease with depth; these opposite trends can be attributed to the volcanic and atmospheric origins of He and other gases, respectively. Also, the simulation results indicate that future He concentration profiles exhibit shifting trends towards lower values. In contrast, the distributions of O₂, Ne, Ar, and Kr concentrations remain almost constant over time. Such different temporal behaviors are primarily due to the much higher diffusion coefficient of He in water compared with that of other gases.

• Advective heat and salt fluxes below a depth of 400 m are three orders of magnitude larger than those at shallower depths, mainly because of the bowl-shaped structure of the lake basin. Also, the simulation results demonstrate that there is a depth interval of about 50 m exactly below the major chemocline, where upwelling flow rates and, consequently, advective fluxes are negative. Although the main chemocline significantly reduces the advective pathway, it has little effect on diffusive fluxes. Another important conclusion from this study is that large values of thermal diffusivity lead to comparable contributions from both advection and diffusion mechanisms to heat transfer, whereas the vertical transport of dissolved substances is mainly dominated by advection. Additionally, our simulations revealed no role of buoyancy in the transport of heat and mass in Lake Kivu over the next 500 years.

• Finally, in response to the main question of this research study, the risk of sudden degassing due to the buoyancy instability-driven overturn and/or supersaturation of the water column in Lake Kivu was predicted to remain very low over the next half a millennium.

Conflicts of interest

The authors declare that there are no conflicts of interest.

Abbreviations

AED2	Aquatic EcoDynamics
CH4	Methane
CNRS	French National Centre for Scientific Research
CO2	Carbon dioxide
DRC	Democratic Republic of Congo
Eawag	Swiss Federal Institute of Aquatic Science and Technology
EDCL	Energy Development Corporation Limited
LKMP	Lake Kivu Monitoring Program
REG	Rwanda Energy Group
RT	Rayleigh–Taylor
TDGP	Total dissolved gas pressure
UFZ	Helmholtz Centre for Environmental Research

Nomenclature

Variables

Az	Areal extent of the lake at depth z, [L^2]
A1, A2, and A3	Correlation constants, [−]
a	A constant that accounts for the effect of polar contribution on the second virial coefficient, [–]
B	Buoyancy production, [L^2 T^−3]; the second virial coefficient, [L^3 N^−1]
B1, B2, and B3	Correlation constants, [–]
b	A constant that accounts for the effect of dimerization on the second virial coefficient, [–]
CCH4	Concentration of CH4, [M M^−1]
CCO2	Concentration of CO2, [M M^−1]
Ci	Concentration of the solute i, [N L^−3]
cε1, cε2, and cε3	Constants of the model, [–]
cP	Reference specific heat of the lake water, [L^2 T^−2 Θ^−1]
Di	Mass diffusivity of the solute i, [L^2 T^−1]
f	Coriolis parameter, [T^−1]
fi	Fugacity of the gas species i, [M L^−1 T^−2]
f^(0), f^(1), and f^(2)	Functions of the reduced temperature, [–]
g	Acceleration due to gravity, [L T^−2]
Ki	Bunsen solubility coefficient of the gas species i, [N M^−1 L^−2 T^2]
k	Turbulent kinetic energy, [L^2 T^−2]
hz	Height of the trapezoid at depth z, [L]
Hgeo	Heat flux from geothermal sources, [M T^−3]
Hsol	Shortwave solar radiation received by the lake, [M T^−3]
N	Number of species; The Brunt-Väisälä frequency, [T^−1]
P	Total pressure of the system, [M L^−1 T^−2]
Pair	Atmospheric pressure, [M L^−1 T^−2]
Pc	Critical pressure, [M L^−1 T^−2]
Pi	Partial pressure of the dissolved species i, [M L^−1 T^−2]
PN2	Partial pressure of N2, [M L^−1 T^−2]
PSeiche	Turbulent kinetic energy produced by internal seiching, [L^2 T^−3]
p	Pressure, [M L^−1 T^−2]
R	Universal gas constant, [M L^2 T^−2 Θ^−1 N^−1]
S	Water salinity, [M M^−1] or [–]
T	Temperature, [Θ]
Tc	Critical temperature, [Θ]
Tr	Reduced temperature, [–]
t	Time, [T]
u	Mean velocity component in the x-direction, [L T^−1]
Vz	Volume of the segment at depth z, [L^3]
v	Mean velocity component in the y-direction, [L T^−1]
vi	Partial molar volume of the dissolved gas i, [L^3 N^−1]
w	Velocity component in the vertical direction, [L T^−1]
z	Depth, [L]

Greek

βCH4	Correction coefficient due to the contribution of CH4, [M M^−1]
βCO2	Correction coefficient due to the contribution of CO2, [M M^−1]
βS	Correction coefficient due to the contribution of dissolved salts, [M M^−1]
ε	Turbulent dissipation rate, [L^2 T^−3]
λ	Thermal conductivity, [M L T^−3 Θ^−1]
ν	Molecular viscosity, [L^2 T^−1]
ν′	Molecular thermal diffusivity, [L^2 T^−1]
νε	Turbulent diffusivity of dissipation rate, [L^2 T^−1]
νk	Turbulent diffusivity of kinetic energy, [L^2 T^−1]
νt	Turbulent viscosity, [L^2 T^−1]
	Turbulent thermal diffusivity, [L^2 T^−1]
ρ	Water density, [M L^−3]
ρ0	Reference density of water, [M L^−3]
ω	Acentric factor, [–]

Data availability

The data supporting the findings of this study are available within the article or its supplementary information (SI). Supplementary information: the equations used to calculate TDGP for dissolved gases in Lake Kivu, along with tabulated constants. See DOI: https://doi.org/10.1039/d5em00513b.

Acknowledgements

The authors acknowledge financial support from Alfajiri Energy Corp. and Mitacs through the Mitacs Accelerate Postdoctoral Fellowship Program. The support of the Department of Chemical and Petroleum Engineering at the Schulich School of Engineering and Angelo Nwigwe, the Mitacs Advisor at the University of Calgary, is also greatly appreciated. Finally, the fruitful discussions and meetings with researchers from Eawag, particularly Dr Martin Schmid and Dr Fabian Bärenbold, are gratefully acknowledged.

References

1. E. T. Degens, R. P. von Herzen, H.-K. Wong, W. G. Deuser and H. W. Jannasch, Lake Kivu: structure, chemistry and biology of an East African rift lake, Geol. Rundsch., 1973, 62(1), 245–277.
2. K. Tietze, Geophysikalische Untersuchung des Kivusees und seiner ungewöhnlichen Methangaslagerstätte - Schichtung, Dynamik und Gasgehalt des Seewassers, Doctoral thesis, Christan-Albrechts-Universität zu Kiel, 1978.
3. R. H. Spigel and G. w. Coulter, Comparison of Hydrology and Physical Limnology of the East African Great Lakes: Tanganyika, Malawi, Victoria, Kivu and Turkana (with reference to some North American Great Lakes), in Limnology, Climatology and Paleoclimatology of the East African Lakes, ed. T. C. Johnson and E. O. Odada, Routledge, 1st edn., 1996.
4. O. S. A. E. Lahmeyer, Bathymetric Survey of Lake Kivu. Final report, Kigali, 1998.
5. K. A. Ross, B. Smets, M. De Batist, M. Hilbe, M. Schmid and F. S. Anselmetti, Lake-level rise in the late Pleistocene and active subaquatic volcanism since the Holocene in Lake Kivu, East African Rift, Geomorphology, 2014, 221, 274–285.
6. K. A. Ross, M. Schmid, S. Ogorka, F. A. Muvundja and F. S. Anselmetti, The history of subaquatic volcanism recorded in the sediments of Lake Kivu; East Africa, J. Paleolimnol., 2015, 54(1), 137–152.
7. G. Boudoire, S. Calabrese, A. Colacicco, P. Sordini, P. Habakaramo Macumu, V. Rafflin, S. Valade, T. Mweze, J.-C. Kazadi Mwepu, F. Safari Habari, T. Amani Kahamire, Y. Mumbere Mutima, J.-C. Ngaruye, A. Tuyishime, A. Tumaini Sadiki, G. Mavonga Tuluka, M. Mapendano Yalire, E.-D. Kets, F. Grassa, W. D’Alessandro, S. Caliro, F. Rufino and D. Tedesco, Scientific response to the 2021 eruption of Nyiragongo based on the implementation of a participatory monitoring system, Sci. Rep., 2022, 12, 7488.
8. M. Isumbisho, H. Sarmento, B. Kaningini, J.-C. Micha and J.-P. Descy, Zooplankton of Lake Kivu, East Africa, half a century after the Tanganyika sardine introduction, J. Plankton Res., 2006, 28(11), 971–989.
9. H. Sarmento, M. Isumbisho and J.-P. Descy, Phytoplankton ecology of Lake Kivu (eastern Africa), J. Plankton Res., 2006, 28(9), 815–829.
10. E. T. Degens and G. Kulbicki, Hydrothermal origin of metals in some East African Rift Lakes, Miner. Depos., 1973, 8(4), 388–404.
11. K. Tietze, M. Geyh, H. Müller, L. Schröder, W. Stahl and H. Wehner, The genesis of the methane in Lake Kivu (Central Africa), Geol. Rundsch., 1980, 69(2), 452–472.
12. K. A. Haberyan and R. E. Hecky, The late Pleistocene and Holocene stratigraphy and paleolimnology of Lakes Kivu and Tanganyika, Palaeogeogr. Palaeoclimatol. Palaeoecol., 1987, 61, 169–197.
13. M. Schoell, K. Tietze and S. M. Schoberth, Origin of methane in Lake Kivu (East-Central Africa), Chem. Geol., 1988, 71(1–3), 257–265.
14. K. Tietze, Lake Kivu Gas Development and Promotion-Related Issues: Safe and Environmentally Sound Exploitation, Final Report, 2000.
15. F. Tassi, O. Vaselli, D. Tedesco, G. Montegrossi, T. Darrah, E. Cuoco, M. Y. Mapendano, R. Poreda and A. D. Huertas, Water and gas chemistry at Lake Kivu (DRC): Geochemical evidence of vertical and horizontal heterogeneities in a multibasin structure, G-Cubed, 2009, 10(2), 1–22.
16. F. Hategekimana, J. D. Ndikuryayo, E. Habimana, T. Mugerwa, K. Christian and R. D. Ed. Rwabuhungu, Lake Kivu Water Chemistry Variation with Depth Over Time, Northwestern Rwanda, Rwanda J. Eng. Sci. Technol. Environ., 2020, 3(1) DOI:10.4314/rjeste.v3i1.5.
17. F. Hategekimana, T. Mugerwa, C. Nsengiyumva, F. V. Byiringiro and D. E. R. Rwatangabo, Geochemical Characterization of Nyamyumba Hot Springs, Northwest Rwanda, AppliedChem, 2022, 2(4), 247–258.
18. M. F. Amisi, M. P. Mulungula, K. T. Kisse, B. C. Muhigirwa, P. Natacha, H. B. Lwikitcha, M. R. Eric, A. B. Désiré, N. Déo, A. Z. Migeni, S. Smith, A. Wüest and T. Lawrence, Current status and strategic way forward for long-term management of Lake Kivu (East Africa), J. Great Lakes Res., 2023, 49(6), 102024.
19. F. C. Newman, Temperature Steps in Lake Kivu: A Bottom Heated Saline Lake, J. Phys. Oceanogr., 1976, 6(2), 157–163.
20. R. W. Schmitt, Double diffusion in oceanography, Annu. Rev. Fluid Mech., 1994, 26, 255–285.
21. D. E. Kelley, H. J. S. Fernando, A. E. Gargett, J. Tanny and E. Özsoy, The diffusive regime of double-diffusive convection, Prog. Oceanogr., 2003, 56(3–4), 461–481.
22. B. Ruddick and A. E. Gargett, Oceanic double-infusion: introduction, Prog. Oceanogr., 2003, 56(3–4), 381–393.
23. M. Schmid, M. Busbridge and A. Wüest, Double-diffusive convection in Lake Kivu, Limnol. Oceanogr., 2010, 55(1), 225–238.
24. F. Bärenbold, B. Boehrer, R. Grilli, A. Mugisha, W. von Tümpling, A. Umutoni and M. Schmid, No increasing risk of a limnic eruption at Lake Kivu: Intercomparison study reveals gas concentrations close to steady state, PLoS One, 2020, 15(8), e0237836.
25. H. Sigurdsson, Gas Bursts from Cameroon Crater Lakes: A New Natural Hazard, Disasters, 1988, 12(2), 131–146.
26. H. Sigurdsson, Lethal gas bursts from Cameroon Crater Lakes, Eos, Trans. Am. Geophys. Union, 1987, 68(23), 570–573.
27. H. Sigurdsson, J. D. Devine, F. M. Tchoua, T. S. Presser, M. K. W. Pringle and W. C. Evans, Origin of the lethal gas burst from Lake Monoun, Cameroun, J. Volcanol. Geotherm. Res., 1987, 31(1–2), 1–16.
28. A. Nayar, Earth science: A lakeful of trouble, Nature, 2009, 460, 321–323.
29. D. S. H. Damas, La stratification thermique et chimique des lacs Kivu Édouard et Ndalaga (Congo Belge), SIL Proc., 1938, 8(3), 51–68.
30. D. M. Schmitz and J. Kufferath, Problèmes posés par la présence de gaz dissous dans les eaux profondes du lac Kivu, Bull. Séances Acad. R. Sci. Colon., 1955, 1, 326–356.
31. E. T. Degens, W. G. Deuser, R. P. von Herzen, H.-K. Wong, F. B. Wooding, H. W. Jannasch and J. W. Kanwisher, Lake Kivu expedition : geophysics, hydrography, sedimentology (preliminary report), 1971.
32. M. Schmid, M. Halbwachs, B. Wehrli and A. Wüest, Weak mixing in Lake Kivu: New insights indicate increasing risk of uncontrolled gas eruption, G-Cubed, 2005, 6(7), 1–11.
33. N. Pasche, M. Schmid, F. Vazquez, C. J. Schubert, A. Wüest, J. D. Kessler, M. A. Pack, W. S. Reeburgh and H. Bürgmann, Methane sources and sinks in Lake Kivu, J. Geophys. Res. Biogeosci., 2011, 116(G3), 1–16.
34. P. Reichert, AQUASIM-A tool for simulation and data analysis of aquatic systems, Water Sci. Technol., 1994, 30(2), 21–30.
35. T. Sommer, M. Schmid and A. Wüest, The role of double diffusion for the heat and salt balance in Lake Kivu, Limnol. Oceanogr., 2019, 64(2), 650–660.
36. F. Bärenbold, Managing Lake Kivu: moving from a steady-state to a dynamic modelling approach, Doctoral thesis, ETH Zurich, 2021.
37. F. Bärenbold, R. Kipfer and M. Schmid, Dynamic modelling provides new insights into development and maintenance of Lake Kivu’s density stratification, Environ. Model. Software, 2022, 147, 105251.
38. G.-H. Goudsmit, H. Burchard, F. Peeters and A. Wüest, Application of k-ϵ turbulence models to enclosed basins: The role of internal seiches, J. Geophys. Res. Oceans, 2002, 107(C12), 1–13.
39. M. Schmid and O. Köster, Excess warming of a Central European lake driven by solar brightening, Water Resour. Res., 2016, 52(10), 8103–8116.
40. A. Gaudard, R. Schwefel, L. R. Vinnå, M. Schmid, A. Wüest and D. Bouffard, Optimizing the parameterization of deep mixing and internal seiches in one-dimensional hydrodynamic models: a case study with Simstrat v1.3, Geosci. Model Dev., 2017, 10(9), 3411–3423.
41. A. Gaudard, L. R. Vinnå, F. Bärenbold, M. Schmid and D. Bouffard, Toward an open access to high-frequency lake modeling and statistics data for scientists and practitioners – the case of Swiss lakes using Simstrat v2.1, Geosci. Model Dev., 2019, 12(9), 3955–3974.
42. M. Schmid and A. Wüest, Stratification, Mixing and Transport Processes in Lake Kivu, in Lake Kivu: Limnology and Biogeochemistry of a Tropical Great Lake, ed. J.-P. Descy, F. Darchambeau and M. Schmid, Springer Dordrecht, 1st edn, 2012, vol. 5, pp. 13–29.
43. H. Burchard, O. Petersen and T. P. Rippeth, Comparing the performance of the Mellor-Yamada and the κ-ε two-equation turbulence models, J. Geophys. Res. Oceans, 1998, 103(C5), 10543–10554.
44. M. Schmid, K. Tietze, M. Halbwachs, A. Lorke, D. McGinnis and A. Wüest, How hazardous is the gas accumulation in Lake Kivu? Arguments for a risk assessment in light of the Nyiragongo Volcano eruption of 2002, Acta Vulcanol., 2002, 14(1–2), 115–122.
45. C.-T. Chen and F. J. Millero, The use and misuse of pure water PVT properties for lake waters, Nature, 1977, 266, 707–708.
46. Y. Zhou, Rayleigh–Taylor and Richtmyer–Meshkov instability induced flow, turbulence, and mixing. I, Phys. Rep., 2017, 720–722, 1–136.
47. Y. Zhou, Rayleigh–Taylor and Richtmyer–Meshkov instability induced flow, turbulence, and mixing. II, Phys. Rep., 2017, 723–725, 1–160.
48. Y. Zhou, R. J. R. Williams, P. Ramaprabhu, M. Groom, B. Thornber, A. Hillier, W. Mostert, B. Rollin, S. Balachandar, P. D. Powell, A. Mahalov and N. Attal, Rayleigh–Taylor and Richtmyer–Meshkov instabilities: A journey through scales, Physica D, 2021, 423, 132838.
49. Y. Zhou, J. D. Sadler and O. A. Hurricane, Instabilities and mixing in inertial confinement fusion, Annu. Rev. Fluid Mech., 2025, 57, 197–225.
50. Y. Zhou, Hydrodynamic Instabilities and Turbulence: Rayleigh–Taylor, Richtmyer–Meshkov, and Kelvin–Helmholtz Mixing, Cambridge University Press, 2024.

---