<figure><img src="image_1.png" /></figure>

# Primary and net ecosystem production in a large lake diagnosed from high-resolution oxygen measurements

Bieito Fernández Castro<fnref n="1" />,<fnref n="2" />, Hannah Elisa Chmiel<fnref n="1" />, Camille Minaudo<fnref n="1" />, Shubham Krishna<fnref n="1" />, Pascal Perolo<fnref n="3" />, Serena Rasconi<fnref n="4" />, Alfred Wüest<fnref n="1" />,<fnref n="5" />

## Key Points:

* Metabolic rates diagnosed with high-resolution free-water oxygen measurements in Lake Geneva

* Two new methods proposed to filter the imprint of vertical dislocations that interferes with the diel oxygen signal

* The summer metabolism of Lake Geneva is net autotrophic and shows both rich vertical structures and temporal dynamics

<footnote n="1">¹Physics of Aquatic Systems Laboratory, Margaretha Kamprad Chair, Institute of Environmental Engineering, École Polytechnique Fédérale de Lausanne, Lausanne, Switzerland.</footnote>

<footnote n="2">²Ocean and Earth Science, National Oceanography Centre, University of Southampton, Southampton, UK</footnote>

<footnote n="3">³Institute of Earth Surface Dynamics, University of Lausanne, Lausanne, Switzerland</footnote>

<footnote n="4">⁴Université de Savoie Mont-Blanc, INRAE, CARRTEL, Thonon les Bains, France</footnote>

<footnote n="5">⁵Eawag, Swiss Federal Institute of Aquatic Science and Technology, Surface Waters Research and Management, Kastanienbaum, Switzerland.</footnote>
**Abstract**

The rates of gross primary production (GPP), ecosystem respiration (R) and net ecosystem production (NEP) provide quantitative information about the cycling of carbon and energy in aquatic ecosystems. In lakes, metabolic rates are often diagnosed from diel oxygen fluctuations recorded with high-resolution sondes. This requires that the imprint of ecosystem metabolism can be separated from that of physical processes. Here, we quantified the vertical and temporal variability of the metabolic rates of a deep, large, mesotrophic lake (Lake Geneva, Switzerland-France) by using a 6-months record (April – October 2019) of high-frequency, depth-resolved (0–30 m) dissolved oxygen measurements. Two new alternative methods (in the time and frequency domain) were used to filter low-frequency basin-scale internal motions from the oxygen signal. Both methods proved successful and yielded consistent metabolic estimates showing net autotrophy (NEP = GPP - R = 55 mmol m⁻² d⁻¹) over the sampling period and depth interval, with GPP (235 mmol m⁻² d⁻¹) exceeding R (180 mmol m⁻² d⁻¹). They also revealed significant temporal variability, with at least two short-lived blooms occurring during calm periods, and a vertical partitioning of metabolism, with stronger diel cycles and positive NEP in the upper ~10 m and negative NEP below, where the diel oxygen signal was dominated by internal motions. The proposed methods expand the range of applicability of the diel oxygen technique to large lakes hosting energetic, low-frequency internal motions, offering new possibilities for unveiling the rich spatio-temporal metabolism dynamics in these systems.

**Plain Language Summary**

Quantifying the rates of production and degradation of organic matter is of fundamental importance for understanding the cycling of carbon in lakes, and for assessing and managing the health of these ecosystems. Owing to the increased availability of high-resolution sondes, this quantification is now frequently achieved by recording changes of oxygen concentration directly in free-water, allowing for an unprecedented temporal coverage. In small, shallow, productive lakes, the application of this method is relatively straightforward because oxygen changes driven by biological processes are larger than vertical transport. Instead, large, deep lakes host large-amplitude vertical motions, which can significantly impact the oxygen record, hindering the calculation of metabolic rates. In this study, we quantified the metabolism of Lake Geneva during the productive period of 2019 with continuous oxygen measurements at different depths. To obtain reliable metabolic rates, we developed new analytical techniques allowing for the isolation of the biological signal. This study unveils the rich spatio-temporal metabolic dynamics of Lake Geneva, which would remain unnoticed with traditional sampling techniques, and provides new insights and mathematical tools for improving the application of the free-water oxygen technique in large and energetic systems.

**1 Introduction**

The term *metabolism* describes the carbon and energy transfer in, out and within a living system, which may refer to the level of individual organisms as well as to the scale of entire ecosystems. In the aquatic environment, the functioning of ecosystem metabolism is often described by the rates of gross primary production (GPP) and ecosystem respiration (R), which together build the metabolic balance of a system, termed as net ecosystem production (NEP = GPP - R). Quantifying the magnitude and variability of these rates has been on the focus of aquatic research for nearly a century, in order to understand how carbon is transferred, processed and allocated within the aquatic environment (Staehr, Christensen, et al., 2012). It is also instrumental for assessing the health of aquatic ecosystems suffering from anthropogenic pressures, such as eutrophication, and for monitoring their recovery when restoration measurements are applied.

Accepted Article
For both monitoring and research purposes, metabolic rates have traditionally been quantified as the rate of change of dissolved oxygen (DO) in the aquatic solution or as the assimilation rate of isotopically labelled carbon (¹⁴C or ¹³C) into biomass in bottle incubations of natural water samples (Gaarder & Gran, 1927; Steeman Nielson, 1951; Hama et al., 1983). Alternatively, metabolic rates are obtained from *in situ* free-water measurements of oxygen concentration using the so-called *diel DO technique* (Odum, 1956). This method relies on the light dependency of photosynthetic primary production to break down temporal changes of DO during day- and night-time into GPP, R, and NEP rates. In recent years, the use of the free-water technique has gained popularity owing to the increased availability of affordable and reliable sensors that allow for the acquisition of *in situ* DO measurements at high-temporal resolution from the moored instruments (Marcé et al., 2016) and autonomous underwater vehicles (Briggs et al., 2018). These developments enable the quantification of metabolism with an unprecedented temporal coverage and resolution, fostering the development of an entire field of research around the diel oxygen technique (Cole et al., 2000; Gelda & Effler, 2002; Staehr & Sand-Jensen, 2007; Staehr et al., 2010).

In addition to the improved temporal coverage, a major advantage of free-water DO measurements over bottle incubation is that the undesired effects of enclosing the plankton unnaturally in a small flow-limiting bottle (Quay et al., 2010) and associated error propagation can be avoided (Staehr, Christensen, et al., 2012). However, free-water measurements come with their own methodological challenges, as physical processes are blended with the biological imprint in the measured oxygen records (Staehr et al., 2010; Antenucci et al., 2013; Dugan et al., 2016). In stratified water bodies, the main physical terms contributing to the DO signal are the gas exchange of oxygen through the air-water interface, the vertical transfer of DO by turbulent diffusion, and the vertical dislocation of water masses by internal baroclinic motions (Coloso et al., 2008; Valerio et al., 2019; Peeters et al., 2019). In order to separate the physical and biological components in the DO signal, the most common approach is to model the physical processes and subtract their contribution from the observations (Staehr et al., 2010; Staehr, Christensen, et al., 2012; Peeters et al., 2016). This approach, which is suited for physical processes that are tractable with relatively simple mathematical formulations, like air-water gas exchange, becomes of limited applicability for more unpredictable processes like heterogeneous turbulent mixing (Coloso et al., 2011), and internal motions (Coloso et al., 2008).

Because most the previous studies on lake metabolism have relied on DO readings from sensors located within the mixed layer, the parameterization of the gas transfer velocity has been a major concern regarding the uncertainty of metabolic rate estimates (Hanson et al., 2008; Cremona et al., 2014; Dugan et al., 2016). Recent studies have extended the mixed layer approach to the entire water column, and showed that meta- and hypolimnetic processes are important for the whole-lake metabolism, particularly in clear-water lakes (Coloso et al., 2008; Staehr, Christensen, et al., 2012; Obrador et al., 2014; Giling et al., 2017). Quantifying metabolic rates below the mixed layer involves dealing with the more intricate signal of internal motions and turbulent diffusion. Although often uncertain, turbulent diffusion can be modelled by applying some reasonable estimate of the vertical diffusion coefficient derived from existing parameterizations (Hondzo & Stefan, 1993), or from the heat-budget method (Powell & Jassby, 1974). Because turbulent diffusion is a continuous process acting in metalimnetic layers, it has been found that the metabolic calculations at the diel scale are often not very sensitive to this parameter (Sadro et al., 2011; Obrador et al., 2014; Giling et al., 2017).

The vertical dislocations caused by internal motions are challenging to model. In small lakes, they appear as high-frequency noise (a few minutes to a few hours), and can be filtered numerically (Coloso et al., 2008; Batt & Carpenter, 2012; Staehr, Testa, et al., 2012; Obrador et al., 2014). The combination of numerical filtering and inverse mod-

Accepted Article
elling approaches to calculate metabolic rates are thus successful in small, weakly en-
ergized lakes, i.e. exposed to light winds and featuring small vertical motions (Obrador
et al., 2014; Staehr, Testa, et al., 2012; Phillips, 2020). However, as the lake size increases,
the time scales of internal motions cover a broader range, including basin-scale seiches
and waves, with time scales ranging from the buoyancy period (a few minutes) to sev-
eral days (Csanady, 1975). These motions potentially overlap with the diel frequency
of GPP and cannot be treated as high-frequency noise, further challenging the calcula-
tion of metabolic rates. In the absence of a well-established method for dealing with this
problem, the application of the DO diel method below the epilimnion of large lakes re-
mains so far very limited (Giling et al., 2017).

In this study, we quantified lake metabolism from a 6-months, depth-resolved, oxy-
gen dataset collected on the large, deep, clear-water Lake Geneva (Switzerland-France).
In order to obtain a meaningful estimate of metabolism in such a large lake, depth-resolved
DO measurements were recorded in the upper 30 m with moored sensors during spring
and summer 2019. Because of its large size, the time scales of internal motions overlap
with the diel cycle, hindering the direct estimation of biological rates, particularly in the
metalimnion. The main goals of this study are: (1) to develop and test new methodolo-
gies that aim to disentangle the biological signal from DO records in the presence of low-
frequency, large internal motions, and (2) to quantify the magnitude and the vertical and
temporal variability of ecosystem metabolism in Lake Geneva during the summer pro-
ductive season. We propose and evaluate the performance of two different approaches,
applied in the frequency and time domain, respectively, to separate the biological from
the physical signals. Both methods rely on the assumption that oxygen fluctuations due
to lake internal motions are chiefly governed by vertical dislocations and, therefore, they
must be correlated with temperature fluctuations.

## 2 Methods

### 2.1 Study site

Lake Geneva is a 582 km² large meso- to eutrophic lake with mean and maximum
depths of 154 and 309 m, respectively. It is located at the border of Switzerland and France
(46° 27 N, 6° 32 E) and consists of a large and deep main basin in its eastern part, and
a small and shallower basin in the western part. The lake is oligo-monomictic and the
thermal stratification persists from spring to late autumn and deepest seasonal mixing
is reached in late winter (February). Complete turnover, however, only occurs about once
every decade (Schwefel et al., 2016). The major tributary to the lake is the Rhne River,
which is responsible for 70–75% of the total inflow. It enters the lake at its eastern part,
and provides an annual mean discharge of ~182 m³ s⁻¹. About 250 m³ s⁻¹ leave the lake
at its outflow in the west, and the theoretical water residence time is estimated at ~11
years. Due to its large size, the currents are affected by Coriolis force and the most en-
ergetic modes of internal motions are Poincaré waves in open waters (Lemmin et al., 2005;
Lemmin, 2020), and Kelvin waves in near-shore regions (Bouffard & Lemmin, 2013), with
periods of ~16 hours and ~4–8 days, respectively. The nearshore areas are also affected
by coastal upwellings, which can persist for several days (Reiss et al., 2020).

The shorelines of Lake Geneva are densely populated and to improve the lake wa-
ter quality, phosphorous supply to the lake has been strictly regulated over the past 40
years. Following severe eutrophication in the 1980's, the lake has ever since recovered
and decreased its annual mean values of total phosphorous concentration to current lev-
els of 12 and 46 µg L⁻¹ in the upper 0–30 m and in the lower 250-309 m, respectively (CIPEL,
2019). Photosynthetic production rates and phytoplankton biomass, however, did not
follow their expected trends, and the temporal dynamics of algae blooms in Lake Geneva
remain a highly debated topic (Tadonléké et al., 2009; Anneville et al., 2019). For sci-
entific purpose, a floating platform (LXPLORE, https://lexplore.info) was installed in

Accepted Article
February 2019. This 100 m² structure is located about 570 m off the northern lake shore, near the town of Pully, where the maximum water depth reaches 110 m (46° 30 N, 6° 39 E).

## 2.2 Moored instruments

In order to monitor biogeochemical and physical dynamics in the productive surface layer (top 30 m), a sensor chain was installed next to the platform. Temperature, DO, and photosynthetically active radiation (PAR) were measured continuously from the sensor chain. Temperature was measured every minute by Minilog-II-T sensors (VEMCO, resolution 0.01°C) at every 2.5 m depth. PAR was measured every minute by six LI-COR sensors (LI-193 SA) located at 0.5, 2.5, 5, 10, 20, and 30 m depth. DO concentration and saturation were measured every 5 minutes by seven optical sensors from different providers: two miniDOT sensors (PME) were placed at 0.5 m and 30 m depth, one multiparameter sonde (EXO2, YSI) containing a DO sensor was attached at 2.5 m depth, two RBR coda-ODO sensors were placed at 5 and 10 m, and two Aanderaa Oxygen optodes at 15 and 20 m depth. The YSI sonde and the RBR-coda sensors were equipped with automated wiper systems, that cleaned the sensor surface prior to each measurement, preventing biofouling at the sensor surface. Furthermore, the YSI sonde and the miniDot at 0.5 m contained a copper shield for protection against biofouling. The calibration of the sensors was checked regularly during maintenance work (every 4-6 weeks) by immersing the sensors together in 100% saturated water under controlled temperature condition in the laboratory. Data were post-corrected for drift, when necessary, assuming linear drift between maintenance dates. Here we focus on a data subset obtained during 19 Apr – 7 Oct 2019. Prior to the metabolic calculations, the mooring time-series were averaged hourly and the DO data was interpolated in the same vertical grid as temperature.

## 2.3 Microstructure and CTD profiles

Additional water properties down to 80–100 m depth were recorded approximately once per week using microstructure profilers. Three different instruments by Rockland Scientific Instruments (Canada) were used: two small self-contained microCTD profilers (SN310, SN158) and a VMP-500 (SN028). In addition to high-resolution microstructure sensors, all three profilers were equipped with accurate conductivity-temperature-depth (CTD) and optical sensors for chlorophyll-a fluorescence and turbidity. Between March 2019 and March 2020, a total of 401 profiles were collected from the LXPLORE platform. Full-depth (down to ~105 m) CTD profiles were carried out with a Sea&Sun CTD75M device, which includes the following sensors: TDI Pt100 (Temperature), 7-pole electrode conductivity cell, OxyGuard dissolved oxygen sensor, pH, Seapoint turbidity meter, Seapoint chlorophyll fluorometer. In total, 99 profiles were collected on a regular basis (first monthly until June 2019 and then weekly or bi-weekly) between February 2019 and April 2020.

## 2.4 Meteo- and hydrological data

A Campbell Scientific Automatic Weather Station, installed on the LXPLORE platform, recorded meteorological parameters in 10 minutes intervals since 8 April 2019. Variables recorded include: air temperature, relative humidity, shortwave radiation, wind speed, wind direction and atmospheric pressure. Radiative shortwave and longwave heat flux calculations were performed using the method proposed in Fink et al. (2014), while turbulent fluxes for momentum and heat were computed following the similarity theory (Zeng et al., 1998; Woolway et al., 2015). Hourly values of temperature and flow of the Rhne River at the inflow were obtained from the station of Porte-du-Scex and provided by the Swiss Federal Office for the Environment (https://www.bafu.admin.ch/).

Accepted Article
## 3 Determination of metabolic rates

### 3.1 Rationale for the correction of baroclinic vertical dislocations

The main challenge we faced for the calculation of depth-resolved metabolic rates from high-frequency (DO) signals in this large lake was to filter out the physical contribution to DO variations due to internal motions. The main assumption underlying the two methods proposed below is that DO fluctuations due to internal motions are dominated by vertical (baroclinic) displacements, and consequently, they are correlated with temperature fluctuations as shown in this section.

223

224

225

226

227

228

229

230

231

232

Using a first-order Taylor expansion, the temperature fluctuations ($T'(z, t) = T - \langle T \rangle$, where $\langle T \rangle$ is the low-passed temperature time-series and represents a mean equilibrium value), at depth $z$ and time $t$, caused by a vertical dislocation $\delta$ due to internal motions (IM) can be approximated by:

$$T'_{\text{IM}}(z, t) \approx -\delta \frac{\partial \langle T \rangle}{\partial z}. \qquad (1)$$

233

234

235

The first order approximation is adequate for temperature because it is a monotonic function of depth for almost the entire lake volume. The vertical dislocation can be estimated from a temperature record as

$$\delta(z,t) \approx -T'_{\text{IM}} \left( \frac{\partial \langle T \rangle}{\partial z} \right)^{-1}. \qquad (2)$$

236

237

238

Similarly, also the DO fluctuations can be approximated by Taylor expansion. Because the DO profile is sometimes not monotonic (it usually exhibits maxima or minima in lakes), a second order approximation may be more adequate in some cases:

$$\mathrm{DO}'_{\mathrm{IM}}(z,t) \approx -\delta \frac{\partial \langle \mathrm{DO} \rangle}{\partial z} + (-\delta)^2 \frac{\partial^2 \langle \mathrm{DO} \rangle}{\partial z^2} \qquad (3)$$

239

240

241

By introducing Eq. 2 into 3 and re-arranging, the DO fluctuations can be expressed as a linear (first-order approximation) or quadratic (second order) function of temperature fluctuations:

$$\mathrm{DO}'_{\mathrm{IM}}(z,t) \approx a_1 T' + a_2 (T')^2 \qquad (4)$$

242

243

244

where $a_1$ and $a_2$ are constants related to the ratio of mean temperature and oxygen vertical gradients. The two methods proposed here use this approximation to filter the DO signal in the frequency (scale separation method) and time domains, respectively.

245

### 3.2 Scale separation method

246

#### 3.2.1 NEP from fortnight-scale DO budget

In this first approach, which we call *scale separation method*, we combined a fortnight-scale oxygen budget to derive NEP, with a spectral approach to remove the physical signal from the diel DO fluctuations and derive GPP and R. For the first step, we recovered the low-frequency seasonal DO signal and eliminated the oscillatory contribution of internal lake motions by low-pass filtering the temperature and DO records with a Gaussian filter with a time scale of 15 days. We consider that DO variations over these time scales are governed by the interplay between NEP, vertical turbulent diffusion, gas exchange with the atmosphere, and the intrusion of the well-oxygenated Rhône River into the lake thermocline during the stratified period (Halder et al., 2013). Due to this intrusion, and to the fact that surface waters leave the lake at the outflow, mass conservation imposes a continuous basin-wide hydraulic uplift.

247

248

249

250

251

252

253

254

255

256

257

Taking this into account, the seasonal NEP for a layer $i = 1, \dots, N$ was computed as the residual between the rate of change of DO concentration on fortnight scales ($\partial \mathrm{DO}/\partial t$)

258

259

Accepted Article
As a boundary condition at the air-lake interface (N+0.5), the diffusive flux was equated to:

$$F_{\text{diff}}^{N+0.5} = A_0 F^{\text{air-lake}} = A_0 k_w (\text{DO}^N - \text{DO}_{\text{eq}}^N) \quad (8)$$

where $A_0$ is the lake surface area and $k_w$ is a parameterized piston velocity calculated using the surface-renewal model by Soloviev et al. (2007), which includes the contribution of wind stress, thermal convection and bubbles to gas exchange. In order to assess the uncertainty associated with the spread of $k_w$ among the available parameterizations (Dugan et al., 2016), a set of parameterizations including the most popular ones was tested. The air-lake fluxes were calculated with non-smoothed values of DO concentration and instantaneous $k_w$, and smoothed afterwards. Turbulent fluxes $F_{\text{diff}}$ between the air-lake interface and the base of the mixed layer were interpolated linearly in depth for each time step, assuming instantaneous homogenization within the mixed layer.

A mixed layer deepening (i.e. entrainment) term was also added to the diffusive flux when a deepening of the mixed layer was observed. Considering that a mixed layer deepens from $h_{\text{mix}}^j$ at time step j to $h_{\text{mix}}^{j+1}$ at time step j+1, the maximum entrainment flux is:

$$F_{\text{ent}}^{\text{max}} = \frac{\langle DO \rangle_{\text{ML}}^j - \langle DO \rangle_{\text{ML}}^{\text{ent}}}{\delta t} \cdot V_{\text{ML}}^j \quad (9)$$

where $\langle DO \rangle_{\text{ML}}^j$ is the volume-weighted mean DO concentration in the mixed layer at step j, $V_{\text{ML}}^j$ is the mixed layer volume at the same time step, and $\Delta t$ is the time interval between the two steps. The mean DO concentration in the mixed layer after entrainment was calculated as:

$$\langle DO \rangle_{\text{ML}}^{\text{ent}} = \frac{\langle DO \rangle_{\text{ML}}^{j} V_{\text{ML}}^{j} + \langle DO \rangle_{\text{TH}}^{j} V_{\text{TH}}^{j}}{V_{\text{ML}}^{j+1}} \quad (10)$$

where TH stands for thermocline, which represents the water volume comprised between $h_{\text{mix}}^j$ and $h_{\text{mix}}^{j+1}$, and, hence, its volume was $V_{\text{TH}}^j = V_{\text{ML}}^{j+1} - V_{\text{ML}}^j$. The entrainment flux is maximal at the bin interface located at $h_{\text{mix}}^j$ and set to decrease linearly towards the surface and $h_{\text{mix}}^{j+1}$.

### 3.2.2 Spectral isolation of the diel signal

The budget analysis allows the quantification of a seasonal-scale NEP, i.e. the fortnight balance between GPP and R. However, to disentangle these terms it is necessary to resolve the day-night oxygen cycle, in which the effect of respiration can be isolated since photosynthesis is only active during daytime (Odum, 1956). We propose a spectral approach to isolate the biological DO signal from physical advection. The spectral approach was introduced by Cox et al. (2015), who determined GPP and R from the amplitude of a DO signal at the diel frequency calculated in the frequency domain. However, with this method, the long-term effect of oxygen accumulation due to NEP is implicitly removed, because the differences in the amplitude of the signal during day and night are not recovered but, at best, modelled by a prescribed function (Cox et al., 2015). This problem is solved in our approach by recovering NEP from the fortnight-scale budget. The newly proposed spectral method yields an estimate of the daily R, and GPP is subsequently computed using the budget-derived NEP as $GPP = NEP + R$ (R is defined positive).

Furthermore, the possible contamination of the diel amplitude by physical signals was not considered by Cox et al. (2015) and is addressed here. The removal of the effect of baroclinic motions was performed with basic spectral analysis techniques (Welch, 1967) in four steps: (i) removal of the fraction of the oxygen variability coherent with temperature variability (Levine & Lueck, 1999), (ii) removal of a spectral baseline around the diel frequency, (iii) rejection of data segments where the de-contaminated spectrum is lower than 1.5 times the baseline, and (iv) rejection of the data segments where the

<figure><img src="image_2.png" /></figure>
diel oxygen increase preceded the light increase by more than 2 hours. R rates were finally computed by integrating the de-noised spectra around the diel frequency. These calculations where done for each time step (1 hour) over 256-d datapoint (~11 days) overlapping segments of DO data at different depths. The details of this procedure are outlined in the Supplementary Material (Text S1).

### 3.3 Diel method with time-domain correction

Alternatively, we introduce a simpler method in which DO variations due to vertical dislocations induced by baroclinic motions are filtered directly in the time domain. With this approach, we could apply the classical diel method (Odum, 1956) directly to a DO signal from which the signature of internal motions has been removed. To do so, the DO records at different depths were split into 24 h segments. For each day, we calculated daily mean temperature and DO profiles (⟨T⟩, ⟨DO⟩) and the T' and DO' anomalies with respect to the mean profiles. Then, DO' was fitted against a second order polynomial of T' (Eq. 4) for each layer, in order to model the fraction of the DO variability that can be explained by vertical dislocations (DO'_{IM}). A clean DO signal for each day and depth was then obtained as DO_{clean} = DO' - DO'_{IM} + ⟨DO⟩. Within each day, the hourly change of DO_{clean} presumably driven by biological production/consumption (hourly NEP, which we denote as nep) was calculated by subtracting the gas exchange and vertical diffusion terms according to Eqs. 5, 7 and 8. Because the fortnight-scale budget indicated that the river contribution and mixed layer deepening are minor (see Results) we ignored these terms here, and the hourly nep for a layer i reduced to:

$$ nep^i = \frac{\partial DO_{clean}^i}{\partial t} - \frac{1}{V^i} (F_{diff}^{i-0.5} - F_{diff}^{i+0.5}) \quad (11) $$

An inverse modelling approach, applied to 24 hours segments of nep, was used to determine the daily rates of GPP and R based on a fit to the light and temperature profiles following Hanson et al. (2008) and Obrador et al. (2014). The determination coefficient (R²) of the DO' signal fit to Eq. 4 was used to reject data segments dominated by physical noise (i.e. with R² > 0.75). For these segments, and for days yielding negative values of GPP or R, we set GPP = NEP (R = 0) when NEP > 0 and R = |NEP| (GPP = 0) when NEP < 0, where NEP was calculated by averaging Eq. 11 over 24 h. This was done to ensure that the oxygen mass is conserved and that the relation between the three metabolic rates (NEP = GPP - R) is always consistent. The details are outlined in the Supplementary Material (Text S2).

## 4 Results

### 4.1 Overview on the dataset

The temperature and oxygen records unraveled the underlying dynamics during the measurement period, which corresponded to the productive season of Lake Geneva (April to October 2019) (Fig. 1). Until mid-May, the water column was weakly stratified with surface temperature exceeding by 2 °C the temperature at 30 m (~8 °C), and the mixed layer was highly dynamic indicating variable and sometimes vigorous mixing conditions (Fig. 1A). From mid-May on, net warming of the water column occurred and near surface temperature increased up to a maximum of ~25 °C in mid-August, increasing water-column stability. The mixed layer depth became shallow varying between 0.5 m (depth of the upper sensor) and 7 m. It was particularly shallow during the last two weeks of May, the last week of June and in mid-July, indicating highly stable conditions. Since mid-August, mixed layer temperatures decreased steadily to ~17 °C at the end of the sampling period. After mid-September, the mixed layer depth began to fluctuate significantly and showed a net deepening.
tially solve this problem but introduce too much noise from advective processes, increas-
ing the uncertainty in the derived NEP rates. We concluded that 15-days smoothing was
the best compromise between removal of the undesired physical noise and preservation
of the biological footprint. The turbulent diffusion term showed a broader range of val-
ues (up to ±30 mmol m⁻³ d⁻¹) and was always divergent (negative) in the upper 10 m,
whereas convergence (positive values) prevailed below (Fig 2D). The contribution of en-
trainment to the diffusive term was overall small (Fig. S4), due to the shallow and weakly
variable mixed layer depth during the summer period. The Rhône intrusion caused DO
changes of up to 10 mmol m⁻³ d⁻¹ in the metalimnion (Fig 2E). The uplift term showed
an opposite pattern due to upward oxygen transport (Fig 2F). Because oxygen concen-
tration was not far from equilibrium at the river intrusion depth (Fig. 1.C), the miss-
balance between the two terms was small.

The diagnosed NEP rates indicated a vertical partition of the lake metabolic bal-
ance into two layers. A positive NEP prevailed in the top ~10 m with typically increas-
ing values towards the surface and maxima of ~25 mmol m⁻³ d⁻¹ (Fig. 3A). Below ~10 m,
a negative NEP dominated, despite some localized positive values at the beginning of
the measurement period in mid-April. A compensation depth of z_c = 8.75 m was de-
termined, defined as the depth where the temporal mean NEP had a value of zero. The
compensation depth was shallower than the euphotic depth (defined by 1% light), which
varied between 15 and 20 m, indicating that, although photosynthesis is still possible at
low light levels, heterotrophic respiration was dominant below z_c. Mean summer NEP
in the top 0–8.75 m layer was 69 ± 18 mmol m⁻² d⁻¹ (± standard deviation). In the
deeper layer, mean NEP equaled –13.6±5.8 mmol m⁻² d⁻¹. The depth-integrated NEP
averaged 56 ± 18 mmol m⁻² d⁻¹, showing that autotrophy dominated the metabolic
balance over the top 30 m of Lake Geneva in summer 2019.

The integration of the DO budget terms in the production-dominated ($z < z_c$), and respiration-dominated ($z > z_c$) layers allows to assess the different contributions to the budget and their temporal variability (Fig. 3B-C). In the upper layer, NEP varied weakly throughout the season, with a minimum of 50 mmol m⁻² d⁻¹ in early May, and several slight local maxima of 90–100 mmol m⁻² d⁻¹ in late May, mid-July and late August (Fig. 2A). Air-lake gas exchange, averaging –65±22 mmol m⁻² d⁻¹, was by far the dominant term, especially since the mixed layer became largely supersaturated during the first NEP maximum by late May. The turbulent diffusion term was important until mid-June (~20 mmol m⁻² d⁻¹), when diffusivities were large, and surface exchange was weaker (Fig. 2B). Downward diffusion through the compensation depth averaged –8.0±7.8 mmol m⁻² d⁻¹, and was the second-most relevant term in the budget. The entrainment term (Eq. 9) contributed little to the turbulent flux. The net rate of DO concentration change averaged –3.1±3.7 mmol m⁻² d⁻¹, representing a minor compensation to the effect of the two dominant terms. The net effect of the Rhône intrusion and hydraulic uplift was very small because of the mutual compensation of the convergence and uplift terms.

In the respiration-dominated layer ($z > z_c$), negative NEP was maximal in June at -25 mmol m⁻² d⁻¹ (Fig. 3B). In this layer, the main contribution to the diagnosed NEP was the negative rate of change during this period, which averaged –9.9±3.1 mmol m⁻² d⁻¹. The layer received oxygen from above due to turbulent diffusion (8.0±7.8 mmol m⁻² d⁻¹), but this process also transferred oxygen downward across the lower edge of the layer (–3.1±7.8 mmol m⁻² d⁻¹), with a resulting net supply of 4.9±6.2 mmol m⁻² d⁻¹, representing the second-most important contribution to the diagnosed negative NEP. This contribution was important before mid-June, and smaller thereafter (Fig. 3C). Again, the net contribution by the Rhône intrusion and uplift was small (–1.2±0.8 mmol m⁻² d⁻¹) due to mutual compensation. However, in this case it was not entirely negligible and represented about 10% of the diagnosed respiration rate.

Accepted Article
## 4.3 Frequency-domain diel method

Figure 4 shows GPP and R rates diagnosed with the spectral method. Both metabolic rates typically decreased with depth. Mean GPP rates were maximal (30 mmol m⁻³ d⁻¹) at the location of the two upper-most sensors at 0.5 and 2.5 m depth (Fig. 4B,C), and then declined strongly to ~7 mmol m⁻³ d⁻¹ at 10 m depth. Below, GPP decreased to a lower degree reaching ~1 mmol m⁻³ d⁻¹ at 30 m. Respiration rates were also maximal at the two upper-most sensors (18 mmol m⁻³ d⁻¹) and decreased with depth but not as sharply as GPP rates, to ~8 mmol m⁻³ d⁻¹ at 10 m. However, R was slightly stronger than GPP in the deeper layers, as dictated by the negative NEP computed with the budget approach. According to the prescribed criteria, the diel biological cycle stood above other sources of variability in ~80% of the segments recorded in the upper 5-10 m. Below, the physical signal became more dominant, and at 20-30 m depth, only ~20%, of the segments showed a defined diel peak.

Mean depth-integrated GPP down to 30 m were 236±118 mmol m⁻² d⁻¹, most of which (184 ± 114 mmol m⁻² d⁻¹) corresponded to depths shallower than the compensation depth ($z_c$ = 8.75 m), and a smaller fraction (22%) took place below (53 ± 47 mmol m⁻² d⁻¹). Ecosystem respiration was also dominated by the upper layer (115 ± 101 mmol m⁻² d⁻¹), but it was more evenly distributed vertically, as the lower layer contributed with 66 ± 49 mmol m⁻² d⁻¹, representing 37% of the total respiration in the upper 30 m (181 ± 108 mmol m⁻² d⁻¹).

Regarding the temporal variability, depth-integrated GPP rates were lowest (~100 mmol m⁻² d⁻¹) at the beginning (before mid-May) and at the end of the sampling period (end of September), with volumetric rates ≤10 mmol m⁻³ d⁻¹ throughout the water column (Fig. 4A,B). Maximum GPP rates coincided with the DO accumulation observed at the end of May and mid/end of June. For these periods, depth-integrated rates were 400–600 mmol m⁻² d⁻¹, and volumetric rates in the shallower sensors exceeded 100 mmol m⁻³ d⁻¹, while at 5 m depth GPP decreased already down to ~30–40 mmol m⁻³ d⁻¹.

## 4.4 Time-domain diel method

Prior to calculating metabolic rates using the classical diel method, the oxygen time-series were de-noised in the time domain. The de-noising procedure involves removing the fraction of the DO signal that can be ascribed to vertical dislocations for each day and sensor, which was obtained by fitting the DO fluctuations to a second order polynomial of temperature fluctuations (see Methods). This is illustrated for a couple of T and DO data segments collected on 22 April at 20 m depth (Fig. 5). This figure shows that the DO fluctuations can be strongly linked to physical processes which, in this case, explained more than 86% of the oxygen variance, even if the signal had the appearance of a diel biological cycle. By applying the correction to the DO record, the cycle pattern completely disappears (Fig. 5C), indicating that biological processes played indeed a minor role in this data segment. The variable imprint of internal motions in the DO diel variability throughout the record is illustrated by the determination coefficient of the polynomial fit (Fig. 6). The variance explained by the temperature fluctuations was low ($R^2 < 0.5$, most of the time) in the upper 2-3 sensors (< 10 m), where primary production was found to be larger. For deeper layers this contribution increased significantly, and baroclinic motions very often dominated ($R^2 > 0.8$) the DO variations. This is consistent with the results from the spectral method indicating only small GPP and few occurrences of diel cycles in this depth range.

Metabolic rates derived from the de-noised DO signal also showed a mostly depth-decreasing distribution, but less sharp compared to the spectral method (Figs. 7 and 8). Time-mean GPP (and NEP) decreased with depth from 22 mmol m⁻³ d⁻¹ (13 mmol m⁻³ d⁻¹) in the two upper-most sensors to 2–3 mmol m⁻³ d⁻¹ (0 to -1 mmol m⁻³ d⁻¹) below 20 m (Fig. 8B). The sharp decrease took place between 2.5 and 10 m depth (Fig. 8B).

Accepted Article
In contrast, the R profile was less variable over depth, and mean R values ranged between 3 and 10 mmol m⁻³ d⁻¹ with a weak maximum located at 5 m depth (Fig. 8B). Consistently with the weaker vertical gradients in GPP and R rates compared to the spectral method, the distribution of the fraction of segments, where a diel cycle was detected according to the prescribed criteria (R² < 0.75, GPP, R > 0), was also flatter (Fig. 8B). Below 15 m, the proportion of "good" segments was very similar for both methods (~20%), but they were less frequent in the upper layers (~40% for the time-domain method, compared to 80% with the spectral method, Fig. 4). This different depth distribution was due to the detection of anomalous cycles, where diagnosed GPP was negative (dark gray steps in Fig. 8B), as the fraction of segments showing weak correlation with temperature (R² < 0.75, light gray steps) closely tracked the vertical distribution of good segments of the spectral method (Fig. 4C).

Mean depth-integrated GPP down to 30 m was 234±100 mmol m⁻² d⁻¹ (± standard deviation of the weekly means). Most of this production (154±80 mmol m⁻² d⁻¹) corresponded to depths shallower than the compensation depth ($z_c$ = 8.75 m), and a smaller fraction (34%) corresponded to deeper layers (80±35 mmol m⁻² d⁻¹). The proportion of deep GPP (34%) was slightly larger in comparison with the spectral method (22%), while the overall GPP was similar. Mean NEP in the upper 30 m depth was 55±44 mmol m⁻² d⁻¹, of which 66±40 mmol m⁻² d⁻¹ corresponded to the depth range above the compensation depth, and -11±42 mmol m⁻² d⁻¹ below. These results are consistent with the scale-separation method.

Regarding temporal variability, maximal GPP (and R) rates were found in the upper 10 m and from mid-May to mid-August (Figs. 7A-B and 8A), with a large maximum at the end of May (weekly mean depth-integral of ~600 mmol m⁻² d⁻¹) and a less pronounced peak (~400 mmol m⁻² d⁻¹) at the end of July. On the other hand, NEP calculated with the time-domain diel method (Fig. 7C) resolved more temporal variability compared to the budget method, but it was also noisier, particularly below 10 m, where significant high-frequency variability was observed. Periods of positive NEP extending beyond the photic zone are probably not realistic and result from physically-driven variability at time scales longer than 24 h.

## 4.5 Sensitivity of the metabolic rates estimates

### 4.5.1 Seasonal budget

The seasonal oxygen budget analysis revealed the crucial role of air-lake gas exchange in the upper productive layer, being by far the largest term in the budget and making the diagnosed NEP rates very dependent on its accurate quantification. However, gas exchange is rarely measured directly and Dugan et al. (2016) showed that different parameterizations of the piston velocity ($k_w$) based on wind speed and other meteorological values can produce diverging results. For the previous analysis we have chosen to use the parameterization proposed by Soloviev et al. (2007), which is based on well established physical theories describing the renewal of the near-surface layer (Lorke & Peeters, 2006), and includes the non-interfacial exchange driven by wave-generated bubbles (Woolf, 1997). Here, we assessed the uncertainty in the derived NEP by resolving the oxygen budget using a set of the most popular piston velocity parameterizations in lake studies (Table 1). Piston velocity, oxygen flux and NEP varied between 3.7 and 9.8 cm h⁻¹, 41 and 114 mmol m⁻² d⁻¹, and 46 and 119 mmol m⁻² d⁻¹, with a mean variability of ~40% around the respective mean values (5.5 cm h⁻¹, 61.8 mmol m⁻² d⁻¹, 66.8 mmol m⁻² d⁻¹). The Vachon and Prairie (2013) parameterization stands out as it produced $k_w$ and flux values about 2-fold larger than the other formulations. By excluding this formulation, the relative errors reduced down to the order of 10%. Due to the unrealistically high $k_w$ yielded by this model, it is not considered in the following discussion.

Accepted Article
### 4.5.2 Frequency-domain diel method

Contrary to NEP, GPP was only weakly sensitive to the choice of the piston velocity parameterization as the relative error of the different parameterizations was only 3.7% (Table 1). In turn, GPP estimates were highly sensitive to the different steps performed to remove physical contribution to the DO signal. This filtering was performed in four different steps, as described in the Methods section and in Text S1. In Table 2, we show how time-mean, depth-integrated GPP rates vary as the filtering steps are applied.

Before any filtering, the mean GPP rate integrated down to 30 m was 601 mmol m⁻² d⁻¹, which is more than twice the final value. Step (i) consisted in removing the fraction of oxygen variance correlated with temperature variance. This first step reduced the unfiltered depth-integrated GPP by 29%. After step (i), we found the DO spectra still showing significant variance around the diel frequency, even in the absence of a clear 24 h peak. This imperfect removal of the physical signal can be due to the statistical uncertainty of the method as well as to any DO variation that is not linearly correlated with temperature, i.e. related to non-monotonous vertical gradients and horizontal heterogeneities. To overcome this deficiency, we performed step (ii) in which we removed the spectral baseline around the 24 h period. With this procedure, we found that in 28.3% of the data segments the diel peak was actually not larger than the background, and GPP rates were reduced by 32.2%.

Step (iii), which imposed a minimal signal-to-noise ratio of 1.5 in the 24 h band, increased the number of segments where the diel signal was considered as non-detectable to 48% of the total, and the GPP rates were reduced by an additional 8%. Finally, segments where the DO signal presented an anomalous timing (oxygen increase during the night, see Fig. 9) were removed during step (iv), increasing the number of bad segments to 58% of the total, and decreasing the integrated GPP rates by an additional 10%. The reduction of GPP rates with the different steps was always largest in the deeper layer (8.75–30 m), where biological rates were weaker and the signature of vertical advection more prominent. At the end of the filtering process, GPP rates integrated in this depth range (52 mmol m⁻² d⁻¹) were a factor of 5.6 smaller compared to the initial value (297 mmol m⁻² d⁻¹), while the reduction in the upper layer was much smaller (factor of 1.7).

A potential limitation of the proposed method, particularly step (i), is that it is sensitive to any correlation between the temperature and oxygen even if it is not driven by vertical dislocations, potentially resulting in spurious over-correction. This could be the case of the uppermost layer (0.5 m), where the temperature also displays a diel cycle as a consequence of direct solar heating (Fig. S5). To evaluate the associated error, we performed a sensitivity test where the correlation correction (step i) was performed for all but the two upper-most sensors (0.5 and 2.5 m, respectively). This resulted in an GPP increase by 6-7%, which is a relatively minor error considering the overall impact of the corrections described above. This is coherent with the weak correlation between *T* and DO for these sensors (Fig. 6), indicating that, even though the *T* and DO diel signals could resemble each other, they are driven by different processes, such that they are not as tightly correlated as when both are driven by internal dislocations.

### 4.5.3 Time-domain diel method

In the time-domain method, the average depth-integrated GPP equaled 476 mmol m⁻² d⁻¹ prior to the filtering steps, which is twice as large as the final value reported above (234 mmol m⁻² d⁻¹). As for the spectral method, GPP was more sensitive to the correction below the compensation depth. Without correction, deep GPP would have been three times larger (251 mmol m⁻² d⁻¹). Other than filtering, we further discarded data by considering that the diel method was not applicable when the DO signal was dominated by physical variability, i.e. when the correlation between DO' and T' exceeded a certain threshold (R² = 0.75). By includ-

597

598

599

600

601

602

603

604

605

606

607

608

609

610

611

612

613

614

615

616

617

618

619

620

621

622

623

624

625

626

627

628

629

630

631

632

633

634

635

636

637

638

639

640

641

642

643

644

645

646

647

Accepted Article
ing this correction, we reduced the overall depth-integrated GPP by 28% from 326 to 234 mmol m⁻² d⁻¹. Again, the correction was stronger below the compensation depth with a 45% reduction. The threshold for R² is somewhat arbitrary and the optimal value is difficult to choose objectively. A more restrictive limit of 50% would further reduce the overall GPP by 17% and the deep GPP by 24%. The number of anomalous diel cycles (computed GPP or R < 0) was remarkably high, ~40%, before filtering the physical signal, but the filtering procedure did not reduce those anomalies. By filtering and restricting the diel calculations to R² < 0.75, the number of usable segment decreased further, and 68.0% of them were discarded. Finally, a sensitivity check, during which the correlation correction was only applied below 2.5 m, indicated that the over-correction of the DO signal due to daytime heating could result in a underestimate of GPP rates by 2-3%.

Overall, this sensitivity analysis highlights the need for removing the influence of baroclinic motions in the DO signal to compute metabolic rates. The rates computed with both the spectral and time-domain methods could be significantly overestimated if the signal was not sufficiently corrected, particularly in the metalimnion, where the biological signal is weak and internal motions are large. Furthermore, it shows that the parameterization of the gas-exchange velocity is the key source of uncertainty for the estimation of NEP at seasonal scales.

# 5 Discussion

## 5.1 The imprint of internal baroclinic motions

We analyzed 172 days of depth-resolved high-resolution observations of dissolved oxygen concentration in the large and deep Lake Geneva with the aim of evaluating the metabolic balance of the lake at different depths during the stratification season. Isolating metabolic processes in the DO signal requires that physical processes are filtered out or modelled and subtracted using mass balance approaches (Staehr, Testa, et al., 2012), in order to derive GPP, R and NEP from the diel signal (Staehr et al., 2010). Previous studies investigating vertical heterogeneities in lake metabolism dealt mostly with three physical processes: air-lake exchange through the lake surface into the epilimnion, changes of the mixed layer depth exchanging water between the epi- and metalimnion, and turbulent diffusion across the mean gradient in the meta- and hypolimnion (Coloso et al., 2008; Staehr, Testa, et al., 2012; Obrador et al., 2014; Giling et al., 2017).

The analysis of our mooring records revealed instead that a large fraction of the temperature and oxygen variability was driven by internal motions (Fig. 1), which caused transient vertical displacements as large as several meters (>10 m during energetic periods). During most of the season, the mixed layer was shallow (<5 m), or even nonexistent. Both temperature and oxygen displayed strong vertical gradients throughout the sampling interval (0–30 m), and the vertical dislocations advected these gradients past the sensors, causing large temporal variations. We estimated that internal motions caused a substantial fraction (60%) of the diel DO variability. This fraction increased with depth from 40–50% in the upper 10 m, to 70–80% in the range of 20–30 m depth. A similar pattern of increasing distortion of the biological DO signal by internal motions with depth was observed in previous studies, due to the light dependency of photosynthesis and the larger stratification in the metalimnion compared to the epilimnion (Coloso et al., 2008; Obrador et al., 2014). However, those studies were carried out in shallow, small and more productive lakes, where internal motions have relatively short periods of a few hours (Staehr, Testa, et al., 2012) and can easily be removed (Coloso et al., 2008). In addition, the biological diel signal in metalimnetic waters of these lakes is sufficiently large and represents a larger fraction of the observed variability (Obrador et al., 2014, up to >50%). Due to the large size of Lake Geneva, the temporal scales of internal motions spanned a broad range from several days to a few minutes, and the vertical extension of the dis-

Accepted Article
locations was much larger, impeding the estimation of metabolic rates with the classical approach, particularly at greater depths.

## 5.2 Seasonal net ecosystem production

In order to recover biological rates from such an intricate signal, we adopted the strategy of separating the time scales in which we computed the different metabolic rates. Because daily NEP depends on a fine miss-balance between GPP and R, i.e. it depends on the oxygen accumulation over the course of one day, its calculation is extremely sensitive to the presence of physical processes affecting the oxygen signal at (sub)daily scales. Due to the presence of Kelvin waves (Bouffard & Lemmin, 2013) and upwelling events (Reiss et al., 2020), the internal motions affected the DO signal at time scales of up to one week or more. For this reason, we used a mass-balance approach applied to a 15-days smoothed DO signal to evaluate NEP on these time scales, which are much longer than the usually employed diel method, but allow obtaining a robust measure of the net lake metabolism on the seasonal scale. NEP was found to decrease with depth, shifting from net autotrophy in the upper ~10 m (the compensation depth was 8.75 m), to net heterotrophy below. The photic layer depth varied between ~12 m in mid-summer and >20 m in spring and early-fall, but was generally deeper than the compensation depth, indicating that respiration dominated over gross production the lower photic zone, despite the availability of light. A similar result was found for the oligotrophic Emerald Lake (California; Sadro et al., 2011), while in more productive lakes the shift to net heterotrophy is often found deeper, closer to the base of the photic zone (Obrador et al., 2014).

The two layers differed not only in the net metabolic balance, but also in the main physical terms contributing to the oxygen budget. The mass-balance in the net autotrophic shallower layer was largely governed by gas exchange, which represented 90% of the computed NEP, whereas in the net heterotrophic layer, the largest terms were the net oxygen decrease in time (72% of the computed net respiration) and turbulent diffusion (36% of the computed net respiration). For the productive layer, these results indicate that the accumulation of oxygen by photosynthesis during the summer season in Lake Geneva is balanced by a net transfer to the atmosphere in time scales of a few days, preventing a seasonal accumulation. Several instances of short-lived subsurface oxygen maxima at ~5 m depth were observed during summer and spring after periods of several days of a persistently shallow mixed layer (Fig. 1B). However, these maxima persisted only for a few days, as they were eroded by entrainment during episodes of mixed layer deepening. Hence, the overall temporal trend indicated rather a net decrease in oxygen concentration in the upper layer of -3.1 mmol m⁻² d⁻¹, which represents a small fraction (<5%) of the computed NEP rate.

The net seasonal oxygen concentration decrease could at least partially be driven by the net warming of the upper layer between April and October. The mean net oxygen flux driven by the seasonal warming can be estimated according to Keeling et al. (1993):

$$F_{\text{DO}}^{\text{thermal, air-lake}} = - \left( \frac{F_{\text{heat}}(z=0) - F_{\text{heat}}(z=z_c)}{c_P} \right) \frac{\partial \text{DO}_{\text{sat}}}{\partial T} \approx 2.6 \text{ mmol m}^{-2} \text{ d}^{-1} \quad (12)$$

where $F_{\text{heat}}(z=0) - F_{\text{heat}}(z=z_c)$ is the net heat flux warming the autotrophic layer from the surface to the compensation depth ($z_c$). The term $\frac{\partial \text{DO}_{\text{sat}}}{\partial T}$ represents the derivative of the DO saturation curve as a function of temperature. The good match between the oxygen rate of change and thermally-driven flux suggests that the relatively small oxygen concentration decrease is driven by warming, whereas the bulk of the gas exchange is coupled to net photosynthesis. The implications of this observation are two-fold. On one side, it designates the possibility of estimating net summer ecosystem production with an accuracy of up to 10% by quantifying the surface oxygen flux alone. The downside of it, however, is that the computed NEP is very sensitive to the gas exchange cal-

699

700

701

702

703

704

705

706

707

708

709

710

711

712

713

714

715

716

717

718

719

720

721

722

723

724

725

726

727

728

729

730

731

732

733

734

735

736

737

738

739

740

741

742

743

744

745

746

747

Accepted Article
culation, in particular to the choice of the piston velocity model. Indeed, we found a large
spread among the NEP estimates derived from different parameterizations. NEP in the
upper autotrophic layer ranged between 46 mmol m⁻² d⁻¹ and 69 mmol m⁻² d⁻¹ de-
pending on the gas exchange model used (i.e. ~15% variation).

Net ecosystem respiration between $z_c$ and 30 m depth during the sampling period (13.6 mmol m⁻² d⁻¹) represented only a small fraction (~20%) of the positive NEP in the productive layer. Hence, the net metabolic balance of the upper 30 m was positive, and amounted to 32–56 mmol m⁻² d⁻¹, depending on the gas exchange model used. Considering 172 sampling days, and an oxygen to carbon molar ratio of 138:106 for both photosynthesis and respiration (Redfield, 1958), this is equivalent to a net carbon accumulation (or export) of 50–89 gC m⁻² per productive season between 0 and 30 m depth (or 73–110 gC m⁻² between the surface and $z_c$). However, this depth range does not account for the full heterotrophic lake metabolism, as hypolimnetic oxygen depletion is also observed in the deeper part of Lake Geneva (Schwefel et al., 2018). By applying the mass-balance budgeting approach to deep oxygen profiles collected with the CTD probe, we found that the seasonal ecosystem respiration between 30 and 100 m depth was 23.2 mmol m⁻² d⁻¹, or 37 gC m⁻² (Fig. S6). The total ecosystem respiration (13.6+23.2 = 36.8 mmol m⁻² d⁻¹) is in very good agreement with the mean summer areal hypolimnetic remineralization derived from depletion rates in the historical time-series (Schwefel et al., 2018; Steinsberger et al., 2020, 2021).

Even when accounting for hypolimnetic respiration down to 100 m, an excess of organic matter production of 13–52 gC m⁻² was found over the productive period. This organic matter is possibly accumulated in the water column in the form of suspended or dissolved matter, or in the sediments, and it could in part be remineralized during the winter season. Sediment oxygen demand and the flux of reduced substances out of the sediment, due to oxic and anoxic respiration, respectively, have been found to contribute with ~10 mmol m⁻² d⁻¹ to summer hypolimnetic oxygen depletion (Schwefel et al., 2018). If sediment remineralization rates are comparable in winter and summer, this process could contribute to close the budget by consuming ~15 gC m⁻² of organic carbon during the winter season. Net organic carbon burial in the sediments (11–20 gC m⁻²; Span et al., 1990; Loizeau et al., 2012; Randlett et al., 2015) and winter water column respiration could well account for the reminding miss-balance. Hence, our summer NEP value is consistent with independent estimates based on hypolimnetic oxygen depletion and sedimentation rates.

## 5.3 Gross primary production and respiration

In order to estimate the rates of GPP and R we built upon the classical diel method that relies on the fact that photosynthesis is only active during the day, while ecosystem respiration is active during both day and night (Odum, 1956). The different timing of the two opposing processes creates a diel fluctuation of oxygen concentrations that can be used to infer GPP and R. We first removed the influence of internal motions on the DO signal around the 24 h period, for which we explored two different approaches in the time and frequency domain, respectively. Both approaches were based on the notion that the DO fluctuations must be correlated to temperature fluctuations if both are due to vertical dislocations caused by baroclinic lake motions.

Both methods yielded a consistent estimate of the mean depth-integrated (0–30 m) GPP rate for the productive season of ~235 mmol m⁻² d⁻¹ (equiv. 370 g C m⁻²) and produced comparable vertical profile shapes, in which GPP decreased over depth. Furthermore, both method results display a similar temporal variability at the seasonal scale, with higher integrated rates between June and August and lower rates in spring and towards fall. However, some differences between both methods were also observed. The vertical decrease in GPP was sharper for the spectral method, for which 22% (52 mmol m⁻² d⁻¹)

Accepted Article
of the integrated GPP took place below $z_c$, instead of 34% (80 mmol m⁻² d⁻¹) for the time-domain method. This difference was larger for respiration. Respiration below $z_c$ represented 36% (66 mmol m⁻² d⁻¹) of the total according to the spectral calculation, in comparison with 56% (103 mmol m⁻² d⁻¹) derived in the time-domain. In addition, some details of the temporal variability were captured differently with the two methods, in particular, the near-surface productivity peak at the beginning of July, found with the spectral method, was weaker when computed in the time-domain.

We anticipate several reasons behind these differences. First, the higher GPP/R rates below $z_c$ derived from the time-domain method occur in a depth range where only a small number (~20%) of 24 h segments display a detectable diel cycle. In order to ensure the conservation of oxygen mass in the cases where no diel cycle was detected, we prescribed GPP = NEP ($R = 0$) when NEP > 0, and $R = |NEP|$ (GPP = 0), when NEP < 0. In the spectral method, we used the NEP derived from the fortnight budget, and the effect of internal motions was effectively removed at all time scales, whereas for the time-domain method, we used the mean NEP in the 24 h de-noised segments. The former estimate is susceptible of containing baroclinic noise at lower frequencies; indeed, NEP estimated with this method is much more variable in the net heterotrophic layer compared with the scale-separation spectral method: the coefficient of variation (CV = standard deviation / seasonal mean) of the depth-integrated NEP was of 1320% with the time-domain method, vs 42% with the spectral method. This noise results in a bias of the GPP/R rates towards high values. The scarcity of observable diel cycles below the upper, most productive layers suggests that the diel method might not be the most appropriate approach to estimate net heterotrophic metabolism in the meta- and hypolimnetic deep layers of lakes, especially if internal motions affect the DO signal. The meta- and hypolimnetic oxygen depletion depends on non-locally produced matter, sinking from upper layers, which decouples R from GPP spatially and explains why respiration can occur in the absence of diel oxygen cycles. Therefore, the imprint of ecosystem respiration is rather manifested as a seasonal decrease in oxygen concentration, partially balanced out by turbulent diffusion. In this sense, the seasonal NEP budget provides a more robust estimate of respiration in the deep layers where GPP is very low.

The lower shallow GPP rates computed with the time-domain method, particularly during the early-July bloom, cannot be explained by the removal of physical noise, which represents a small fraction of the biological signal in this depth range. Instead, they can be linked to the timing of the diel DO signal in relation to the light cycle within the upper layers. The spectral analysis revealed a rather small phase difference between both signals in the upper sensors, of typically less than 2 hours (Fig. 9). This implies that oxygen increases as irradiance increases during the morning, it peaks at noon and decreases over the afternoon and evening even if light is still available, but little further decrease is observed during the night (see an example for the July bloom in Fig. S7). The diel approach is based on the assumption that oxygen is produced as long as there is light (see Eq. 7 in Text S2), and equates nighttime to daytime respiration. As a consequence, a phase difference between the DO and light signal of ~6 h is required for the method to provide non-zero values of the metabolic rates. The physiological mechanism causing the unexpected timing for the DO accumulation requires further investigation. It could be related to photoinhibition of the shallow communities exposed to high irradiance levels around noon (Platt et al., 1980), and to the possible enhancement of respiration during daytime (Pace & Prairie, 2005; Tobias et al., 2007; McNair et al., 2013). If photoinhibition is important, the GPP and R rates derived with the diel method are likely an underestimate, even when using a non-linear light dependency in the photosynthesis model, which accounts for light saturation but not photoinhibition (Hanson et al., 2008). The spectral computation has weaker constrains regarding the timing of the DO signal (we only required the phase difference to be positive, or more exactly, > -2 hours) and possibly produces more reliable estimates of metabolic rates for diel cycles with anomalous timing.

Accepted Article
5.4 GPP comparison with previous estimates

Primary productivity in Lake Geneva has been monitored for decades using bot-
tle incubations of water samples collected bi-weekly or monthly at a fixed, central sta-
tion (SHL2), following the ¹³C/¹⁴C methods (Anneville & Pelletier, 2000; Tadonléké et
al., 2009). More recent estimates derived from the combination of remotely-sensed op-
tical water properties and a bio-optical model allow for a better spatio-temporal reso-
lution of productivity patterns, including short-lived and spatially localized blooms (Soomets
et al., 2019). The annual production of organic carbon derived from bottle incubations
averaged 280±76 gC m⁻² y⁻¹ for the period 1971-2005 (Tadonléké et al., 2009). A sim-
ilar value of 302±20 gC m⁻² y⁻¹ was found for the period 2003-2011, both with bot-
tle incubations and remotely sensed data (Soomets et al., 2019). Our estimate of GPP
for the productive season of 2019 is somewhat larger (370 gC m⁻²), even if it does not
consider wintertime GPP. According to Soomets et al. (2019), wintertime productivity
rates are about half of summertime values, so that their productivity values for the strat-
ified season would correspond to about ~200 gC m⁻². Hence, our GPP is close to twice
as large as the previous estimates. However, these observations were carried out ten years
ago and before, and may not be directly comparable to our observations due to inter-
annual or long-term variations. For a closer examination, we directly compared ¹³C PP
measurements acquired during our sampling period with weekly averages of GPP derived
from the spectral and time-domain methods (Fig. S8). This comparison shows that our
*in situ* estimates are generally larger by about a factor of two or more during the most
productive periods, although they tended to converge towards fall, when the productiv-
ity was lower. Overall, mean summer productivity from the bottle incubation profiles
equaled 170 gC m⁻² y⁻¹, which is about half of the *in situ* estimates for the correspond-
ing weeks.

Disagreement between *in situ* and *in vitro* productivity estimates has been found
in several previous studies in lakes (Melack, 1982; Lauster et al., 2006) and ocean (Williams
et al., 2013). These differences could in part be related to the higher temporal resolu-
tion of sonde-derived data, which allow capturing short-lived blooms missed by the bi-
weekly sampling (Soomets et al., 2019). However, during the most productive periods,
the *in vitro* estimates were generally below the weekly standard deviation of the *in situ*
estimates (Fig. S8), suggesting that there might be other reasons for systematic bias be-
yond temporal resolution. Melack (1982) and Lauster et al. (2006) suggested that pelagic
sonde measurements may implicitly include the contribution of littoral productivity. Our
sampling point was at a 570 m horizontal distance from the lake shore, but the lake slopes
are steep in the area, and there are no shallow areas displaying significant benthic pro-
duction nearby. Therefore, the signal of littoral productivity would have to propagate
from distant areas in significantly less than 24 h in order to be recorded as a diel DO fluc-
tuation, which is highly unlikely due to the large lake size. Nevertheless, Soomets et al.
(2019) found a significant heterogeneity in pelagic productivity in Lake Geneva. The SHL2
sampling point, where water samples are collected, is located in the center of the main
lake basin, which features slightly weaker productivity than the lake average (Soomets
et al., 2019), although this is insufficient to account for a two-fold difference (Kiefer et
al., 2015).

As an alternative justification, it must be noted that *in vitro* and *in situ* methods
are inherently different and they capture different components of the planktonic metabolism.
First of all, *in vitro* measurements involve enclosing the planktonic community in a small
volume, which may be problematic to extrapolate to the system level (Schindler, 1998).
In addition, respiratory processes are accounted for differently. While in the *diel* *in situ*
method, respiration is removed from the *diel* NEP value in order to obtain a GPP es-
timate, the incubation method measures the net fixation of carbon into particulate or-
ganic matter over a few hours of incubation (Perga et al., 2015). Although part of the
heterotrophic community composed of larger organisms is excluded due to the small sam-

Accepted Article
ple size, respiration by autotrophs (and possibly small heterotrophs) still occurs within
the bottles (Marra, 2009), and the derived rates should be interpreted as a measure of
net primary production (NPP, i.e. GPP minus R by autotrophs), and hence it is expected
to be smaller than GPP (Bender et al., 1987; Marra, 2009; Pei & Laws, 2013; Briggs et
al., 2018). Furthermore, because only particulate carbon is measured, the exudation of
recently-fixed dissolved organic carbon (DOC) is also not quantified in the ¹³C/¹⁴C method
(Marra, 2009). Excretion of DOC is typically low in freshwater ecosystems (~15%, Baines
and Pace (1991)), however, it may account for as much as 50% of carbon fixation in olig-
otrophic conditions (Fouilland & Mostajir, 2010; Morana et al., 2014), such as those ex-
perienced in Lake Geneva, where inorganic phosphorous concentrations are drawn down
to the detection limit (Tadonléké et al., 2009).

# 6 Conclusions

In conclusion, we found that dissolved oxygen measurements, recorded in a large and deep lake at multiple fixed depths, were highly affected by internal motions. However, our methods proved to be useful in removing the physical signal at different frequencies and to reveal the underlying metabolic dynamics. We isolated a distinct temporal pattern of productivity with two major bloom phases occurring during the middle of summer, which remained rather unnoticed in the traditional bi-weekly monitoring data set. Discrepancies in method results (i.e. from time-domain, spectral-domain, and bottle incubation methods) require further investigation and might be resolved when metabolic rates are linked together with their environmental drivers. The interpretation of field observations in the light of ecological plankton models offers interesting possibilities to disentangle the existing complex interplay between the environmental drivers and lake metabolism (Krishna et al., 2021).

The results presented here can provide practical guidelines on expectations and sensitivities for researchers aiming to quantify vertically-resolved metabolism, particularly in large energetic lakes with large vertical (baroclinic) dislocations of isolines. First, they highlight the importance of oxygen exchange with the atmosphere, not only at daily scales but, more importantly, for the quantification of seasonal-scale net ecosystem production. Observational, direct quantification of air-lake gas fluxes is limited, and the existing parameterizations can yield scattered results. Further research effort is clearly required on this topic. We also showed that the diel-scale approach is likely not optimal for quantifying heterotrophic metabolism in meta- and hypolimnetic layers, where ecosystem respiration can occur in the absence of distinct diel oxygen cycles. Instead, longer integration time scales were shown to produce more robust estimates. Finally, the methods proposed for filtering the signal of internal motions could be widely used in future studies as they rely on fundamental dynamical considerations and they are relatively simple to implement, particularly in the time-domain. Although the more complex spectral approach performed better in our case, both methods are essentially equivalent, and the performance of the time-domain method could be improved by increasing the length of the data segments over which it is applied, achieving a better removal of low-frequency motions and increased statistic significance of the metabolic model fit.

# Acknowledgments

We acknowledge Sb Sebastian Lavanchy for his outstanding technical support; Natacha Pasche for her commitment in the management and operation of the LéXPLORE Platform; Aurélien Balu, Nicolas Escoffier and Lucas Serra for participating in the mooring maintenance and data collection. We would like to thank the entire team from LXPLORE platform, for their administrative and technical support and for LXPLORE core dataset. We also acknowledge LXPLORE five partner institutions: Eawag, EPFL, University of Geneva, University of Lausanne and CARRTEL (INRAE-USMB). We are grateful to Biel Obrador and Emilio Marañón for their scientific advice. We acknowledge the CIPEL and OLA

Accepted Article
for the ¹³C data. This research was funded by the Swiss National Science Foundation grant 200021-179123 (*Primary Production Under Oligotrophicication in Lakes*). B. F. C. was also supported by the European Union's Horizon 2020 research and innovation program under the Marie Skodowska-Curie grant agreement No. 834330 (SO-CUP). The CTD and weather station datasets is available at https://www.datalakes-eawag.ch/. The temperature, PAR, and oxygen mooring data and the computed rates can be downloaded from https://drive.switch.ch/index.php/s/atGN0s2DCW2QDaZ. ADD CIPEL as: For data concerning Lake Geneva-SHL2 we acknowledge: OLA-IS, AnaEE-France, INRAE of Thonon-les-Bains, CIPEL (Rimet et al., 2020).

References

Anneville, O., Chang, C. W., Dur, G., Souissi, S., Rimet, F., & Hsieh, C.-h. (2019). The paradox of re-oligotrophicication: the role of bottomup versus topdown controls on the phytoplankton community. *Oikos*, *128*(11), 1666-1677. doi: 10.1111/oik.06399

Anneville, O., & Pelletier, J. P. (2000). Recovery of Lake Geneva from eutrophication: quantitative response of phytoplankton. *Arch. Hidrobiol.*, **148**(4), 607-624. doi: 10.1127/archiv-hydrobiol/148/2000/607

Antenucci, J. P., Tan, K. M., Eikaas, H. S., & Imberger, J. (2013). The importance of transport processes and spatial gradients on in situ estimates of lake metabolism. *Hydrobiologia*, **700**(1), 9-21. doi: 10.1007/s10750-012-1212-z

Baines, S. B., & Pace, M. L. (1991). The production of dissolved organic matter by phytoplankton and its importance to bacteria: Patterns across marine and freshwater systems. *Limnol. Oceanogr.*, **36**(6), 1078-1090. doi: 10.4319/lo.1991.36.6.1078

Batt, R. D., & Carpenter, S. R. (2012). Free-water lake metabolism: Addressing noisy time series with a Kalman filter. *Limnol. Oceanogr. Methods*, **10**(1), 20-30. doi: 10.4319/lom.2012.10.20

Bender, M., Grande, K., Johnson, K., Marra, J., Williams, P. J. B., Sieburth, J., ... Heinemann, K. (1987). A comparison of four methods for determining planktonic community production. *Limnol. Oceanogr.*, **32**(5), 1085-1098. doi: 10.4319/lo.1987.32.5.1085

Bouffard, D., & Lemmin, U. (2013). Kelvin waves in Lake Geneva. *J. Great Lakes Res.*, **39**(4), 637-645. doi: 10.1016/j.jglr.2013.09.005

Briggs, N., Guemundsson, K., Cetinić, I., D'Asaro, E., Rehm, E., Lee, C., & Perry, M. J. (2018). A multi-method autonomous assessment of primary productivity and export efficiency in the springtime North Atlantic. *Biogeosciences*, **15**(14), 4515-4532. doi: 10.5194/bg-15-4515-2018

CIPEL. (2019). *Rapports sur les études et recherches entreprises dans le bassin Lémanique* (Tech. Rep.). Thonon les Bains: Conseil scientifique de la Commission internationale pour la protection des eaux du Léman.

Cole, J. J., & Caraco, N. F. (1998). Atmospheric exchange of carbon dioxide in a low-wind oligotrophic lake measured by the addition of SF6. *Limnol. Oceanogr.*, **43**(4), 647-656. doi: 10.4319/lo.1998.43.4.0647

Cole, J. J., Pace, M. L., Carpenter, S. R., & Kitchell, J. F. (2000). Persistence of net heterotrophy in lakes during nutrient addition and food web manipulations. *Limnol. Oceanogr.*, **45**(8), 1718-1730. doi: 10.4319/lo.2000.45.8.1718

Coloso, J. J., Cole, J. J., Hanson, P. C., & Pace, M. L. (2008). Depth-integrated, continuous estimates of metabolism in a clear-water lake. *Can. J. Fish. Aquat. Sci.*, **65**(4), 712-722. doi: 10.1139/F08-006

Coloso, J. J., Cole, J. J., & Pace, M. L. (2011). Short-term variation in thermal stratification complicates estimation of lake metabolism. *Aquat. Sci.*, **73**(2), 305-315. doi: 10.1007/s00027-010-0177-0

Cox, T. J., Maris, T., Soetaert, K., Kromkamp, J. C., Meire, P., & Meysman, F.

<figure><img src="image_3.png" /></figure>
(2015). Estimating primary production from oxygen time series: A novel approach in the frequency domain. *Limnol. Oceanogr. Methods*, *13*(10), 529-552.
doi: 10.1002/lom3.10046

Cremona, F., Laas, A., Nõges, P., & Nõges, T. (2014). High-frequency data within a modeling framework: On the benefit of assessing uncertainties of lake metabolism. *Ecol. Modell.*, *294*, 27-35. doi: 10.1016/j.ecolmodel.2014.09.013

Crusius, J., & Wanninkhof, R. (2003). Gas transfer velocities measured at low wind speed over a lake. *Limnol. Oceanogr.*, *48*(3), 1010-1017. doi: 10.4319/lo.2003.48.3.1010

Csanady, G. T. (1975). Hydrodynamics of large lakes. *Annu. Rev. Fluid Mech.*, *7*(1), 357-386. doi: 10.1146/annurev.fl.07.010175.002041

Dugan, H. A., Woolway, R. I., Santoso, A. B., Corman, J. R., Jaimes, A., Nodine, E. R., ... Weathers, K. C. (2016). Consequences of gas flux model choice on the interpretation of metabolic balance across 15 lakes. *Int. Waters*, *6*(4), 581-592. doi: 10.1080/iw-6.4.8.36

Fink, G., Schmid, M., Wahl, B., Wolf, T., & Wüest, A. (2014). Heat flux modifications related to climate-induced warming of large European lakes. *Water Resour. Res.*, *50*(3), 2072-2085. doi: 10.1002/2013WR014448

Foulland, E., & Mostajir, B. (2010). Revisited phytoplanktonic carbon dependency of heterotrophic bacteria in freshwaters, transitional, coastal and oceanic waters. *FEMS Microbiol. Ecol.*, *73*(3), 419-429. doi: 10.1111/j.1574-6941.2010.00896.x

Gaarder, T., & Gran, H. (1927). Investigations of the production of plankton in the Oslo Fjord. *Cons. Perm. Int. Pour L'Exploration La Mer Rapp. Proces-Verbaux Des Reun.*, *42*, 48.

Gelda, R. K., & Effler, S. W. (2002). Metabolic rate estimates for a eutrophic lake from diel dissolved oxygen signals. *Hydrobiologia*, *485*, 51-66. doi: 10.1023/A:1021327610570

Giling, D. P., Staehr, P. A., Grossart, H. P., Andersen, M. R., Boehrer, B., Escot, C., ... Obrador, B. (2017). Delving deeper: Metabolic processes in the metalinmion of stratified lakes. *Limnol. Oceanogr.*, *62*(3), 1288-1306. doi: 10.1002/lno.10504

Halder, J., Decrouy, L., & Vennemann, T. W. (2013). Mixing of Rhône River water in Lake Geneva (Switzerland-France) inferred from stable hydrogen and oxygen isotope profiles. *J. Hydrol.*, *477*, 152-164. doi: 10.1016/j.jhydrol.2012.11.026

Hama, T., Miyazaki, T., Ogawa, Y., Iwakuma, T., Takahashi, M., Otsuki, A., & Ichimura, S. (1983). Measurement of photosynthetic production of a marine phytoplankton population using a stable 13C isotope. *Mar. Biol.*, *73*(1), 31-36. doi: 10.1007/BF00396282

Hanson, P. C., Carpenter, S. R., Kimura, N., Wu, C., Cornelius, S. P., & Kratz, T. K. (2008). Evaluation of metabolism models for free-water dissolved oxygen methods in lakes. *Limnol. Oceanogr. Methods*, *6*(9), 454-465. doi: 10.4319/lom.2008.6.454

Hondzo, M., & Stefan, H. G. (1993). Lake water temperature simulation model. *J. Hydraul. Eng.*, *119*(11), 1251-1273. doi: 10.1061/(ASCE)0733-9429(1993)119:11(1251)

Keeling, R. F., Najjar, R. P., Bender, M. L., & Tans, P. P. (1993). What atmospheric oxygen measurements can tell us about the global carbon cycle. *Global Biogeochem. Cycles*, *7*(1), 37-67. doi: 10.1029/92GB02733

Kiefer, I., Odermatt, D., Anneville, O., Wüest, A., & Bouffard, D. (2015). Application of remote sensing for the optimization of in-situ sampling for monitoring of phytoplankton abundance in a large lake. *Sci. Total Environ.*, 527-528, 493-506. doi: 10.1016/j.scitotenv.2015.05.011

Krishna, S., Ulloa, H. N., Kerimoglu, O., Minaudo, C., Anneville, O., & Wüest, A. (2021). Model-based data analysis of the effect of winter mixing on primary

<figure><img src="image_4.png" /></figure>
production in a lake under reoligotrophication. *Ecol. Modell.*, **440**, 109401. doi: 10.1016/j.ecolmodel.2020.109401

1065

1066

Lauster, G. H., Hanson, P. C., & Kratz, T. K. (2006). Gross primary production and respiration differences among littoral and pelagic habitats in northern Wisconsin lakes. *Can. J. Fish. Aquat. Sci.*, **63**(5), 1130–1141. doi: 10.1139/F06-018

1068

1069

1070

Lemmin, U. (2020). Insights into the dynamics of the deep hypolimnion of Lake Geneva as revealed by long-term temperature, oxygen, and current measurements. *Limnol. Oceanogr.*, **65**(9), 2092–2107. doi: 10.1002/lno.11441

1071

1072

1073

Lemmin, U., Mortimer, C. H., & Bäuerle, E. (2005). Internal seiche dynamics in Lake Geneva. *Limnol. Oceanogr.*, **50**(1), 207–216. doi: 10.4319/lo.2005.50.1.0207

1074

1075

1076

Levine, E. R., & Lueck, R. G. (1999). Turbulence measurement from an autonomous underwater vehicle. *J. Atmos. Ocean. Technol.*, **16**, 1533–1544.

1077

1078

Loizeau, J. L., Girardclos, S., & Dominik, J. (2012). Taux d'accumulation de sédiments récents et bilan de la matière particulaire dans le Léman (Suisse - France). *Arch. des Sci.*, **65**(1-2), 81–92.

1079

1080

1081

Lorke, A., & Peeters, F. (2006). Toward a unified scaling relation for interfacial fluxes. *J. Phys. Oceanogr.*, **36**(5), 955–961. doi: 10.1175/jpo2903.1

1082

1083

Marcé, R., George, G., Buscarinu, P., Deidda, M., Dunalska, J., De Eyto, E., ... Jennings, E. (2016). Automatic high frequency monitoring for improved lake and reservoir management. *Environ. Sci. Technol.*, **50**(20), 10780–10794. doi: 10.1021/acs.est.6b01604

1084

1085

1086

1087

Marra, J. (2009). Net and gross productivity: Weighing in with 14C. *Aquat. Microb. Ecol.*, **56**(2-3), 123–131. doi: 10.3354/ame01306

1088

1089

McNair, J. N., Gereaux, L. C., Weinke, A. D., Sesselmann, M. R., Kendall, S. T., & Biddanda, B. A. (2013). New methods for estimating components of lake metabolism based on free-water dissolved-oxygen dynamics. *Ecol. Modell.*, **263**, 251–263. doi: 10.1016/j.ecolmodel.2013.05.010

1090

1091

1092

1093

Melack, J. M. (1982). Photosynthetic activity and respiration in an equatorial African soda lake. *Freshw. Biol.*, **12**(4), 381–400. doi: 10.1111/j.1365-2427.1982.tb00631.x

1094

1095

1096

Morana, C., Sarmento, H., Descy, J. P., Gasol, J. M., Borges, A. V., Bouillon, S., & Darchambeau, F. (2014). Production of dissolved organic matter by phytoplankton and its uptake by heterotrophic prokaryotes in large tropical lakes. *Limnol. Oceanogr.*, **59**(4), 1364–1375. doi: 10.4319/lo.2014.59.4.1364

1097

1098

1099

1100

Obrador, B., Staehr, P. A., & Christensen, J. P. (2014). Vertical patterns of metabolism in three contrasting stratified lakes. *Limnol. Oceanogr.*, **59**(4), 1228–1240. doi: 10.4319/lo.2014.59.4.1228

1101

1102

1103

Odum, H. T. (1956). Primary production in flowing waters. *Limnol. Oceanogr.*, **1**(2), 102–117. doi: 10.4319/lo.1956.1.2.0102

1104

1105

Pace, M. L., & Prairie, Y. T. (2005). Respiration in lakes. In *Respir. aquat. ecosyst.* (Vol. 15, pp. 103–121). Oxford University Press. doi: 10.1093/acprof:oso/9780198527084.003.0007

1106

1107

1108

Peeters, F., Atamanchuk, D., Tengberg, A., Encinas-Fernández, J., & Hofmann, H. (2016). Lake Metabolism: Comparison of lake metabolic rates estimated from a diel CO2-and the common diel O2- Technique. *PLoS One*, **11**(12), 1–24. doi: 10.1371/journal.pone.0168393

1109

1110

1111

1112

Peeters, F., Hofmann, H., & Fernández, J. E. (2019). On the calculation of lake metabolic rates: Diel O2 and 18/16O technique. *Water Res.*, **165**, 114990. doi: 10.1016/j.watres.2019.114990

1113

1114

1115

Pei, S., & Laws, E. A. (2013). Does the 14C method estimate net photosynthesis? Implications from batch and continuous culture studies of marine phytoplankton. *Deep. Res. Part I Oceanogr. Res. Pap.*, **82**, 1–9. doi: 10.1016/j.dsr.2013.07.011

1116

1117

1118

1119
1120

1121

1122

1123

1124

1125

1126

1127

1128

1129

1130

1131

1132

1133

1134

1135

1136

1137

1138

1139

1140

1141

1142

1143

1144

1145

1146

1147

1148

1149

1150

1151

1152

1153

1154

1155

1156

1157

1158

1159

1160

1161

1162

1163

1164

1165

1166

1167

1168

1169

1170

1171

1172

1173

1174

Perga, M.-E., Tadonléké, R., & Perney, P. (2015). Measures de la production primaire. Transition des protocoles 14C - 13C. In *Rapp. comm. int. prot. eaux léman contre pollut., camp. 2014* (pp. 163-171).

Phillips, J. S. (2020). Time-varying responses of lake metabolism to light and temperature. *Limnol. Oceanogr.*, 65(3), 652-666. doi: 10.1002/lno.11333

Platt, T., Gallegos, C. L., & Harrison, W. G. (1980). Photoinhibition of Photosynthesis in Natural Assemblages of Marine Phytoplankton. *J. Mar. Res.*, 38(4), 687-701.

Powell, T., & Jassby, A. (1974). The estimation of vertical eddy diffusivities below the thermocline in lakes. *Water Resour. Res.*, 10(2), 191-198. doi: 10.1029/WR010i002p00191

Quay, P. D., Peacock, C., Bjrkman, K., & Karl, D. M. (2010). Measuring primary production rates in the ocean: Enigmatic results between incubation and non-incubation methods at Station ALOHA. *Global Biogeochem. Cycles*, 24(3), GB3014. doi: 10.1029/2009GB003665

Randlett, M. E., Sollberger, S., Del Sontro, T., Müller, B., Corella, J. P., Wehrli, B., & Schubert, C. J. (2015). Mineralization pathways of organic matter deposited in a river-lake transition of the Rhone River Delta, Lake Geneva. *Environ. Sci. Process. Impacts*, 17(2), 370-380. doi: 10.1039/c4em00470a

Read, J. S., Hamilton, D. P., Desai, A. R., Rose, K. C., MacIntyre, S., Lenters, J. D., ... Wu, C. H. (2012). Lake-size dependency of wind shear and convection as controls on gas exchange. *Geophys. Res. Lett.*, 39(9), L09405. doi: 10.1029/2012GL051886

Redfield, A. C. (1958). The biological control of chemical factors in the environment. *Am. Sci.*, 46(3), 205-221. doi: 10.5194/bg-11-1599-2014

Reiss, R. S., Lemmin, U., Cimatoribus, A. A., & Barry, D. A. (2020). Wintertime coastal upwelling in Lake Geneva: An efficient transport process for deepwater renewal in a large, deep lake. *J. Geophys. Res. Ocean.*, 125, e2020JC016095. doi: 10.1029/2020JC016095

Rimet, F., Anneville, O., Barbet, D., Chardon, C., Crépin, L., Domaizon, I., ... Monet, G. (2020). The Observatory on LAkes (OLA) database: Sixty years of environmental data accessible to the public. *J. Limnol.*, 78(2), 164-178. doi: 10.4081/jlimnol.2020.1944

Sadro, S., Melack, J. M., & MacIntyre, S. (2011). Depth-integrated estimates of ecosystem metabolism in a high-elevation lake (Emerald Lake, Sierra Nevada, California). *Limnol. Oceanogr.*, 56(5), 1764-1780. doi: 10.4319/lo.2011.56.5.1764

Schindler, D. W. (1998). Replication versus realism: The need for ecosystem-scale experiments. *Ecosystems*, 1(4), 323-334. doi: 10.1007/s100219900026

Schwefel, R., Gaudard, A., Wüest, A., & Bouffard, D. (2016). Effects of climate change on deepwater oxygen and winter mixing in a deep lake (Lake Geneva): Comparing observational findings and modeling. *Water Resour. Res.*, 52(11), 8811-8826. doi: 10.1002/2016WR019194

Schwefel, R., Steinsberger, T., Bouffard, D., Bryant, L. D., Müller, B., & Wüest, A. (2018). Using small-scale measurements to estimate hypolimnetic oxygen depletion in a deep lake. *Limnol. Oceanogr.*, 63, S54-S67. doi: 10.1002/lno.10723

Soloviev, A., Donelan, M., Graber, H., Haus, B., & Schlüssel, P. (2007). An approach to estimation of near-surface turbulence and CO2 transfer velocity from remote sensing data. *J. Mar. Syst.*, 66(1-4), 182-194. doi: 10.1016/j.jmarsys.2006.03.023

Soomets, T., Kutser, T., Wüest, A., & Bouffard, D. (2019). Spatial and temporal changes of primary production in a deep peri-alpine lake. *Inl. Waters*, 9(1), 49-60. doi: 10.1080/20442041.2018.1530529

Span, D., Arbouille, D., Howa, H., & Vernet, J. P. (1990). Variation of nutrient

<figure><img src="image_5.png" /></figure>
1175

1176

1177

1178

1179

1180

1181

1182

1183

1184

1185

1186

1187

1188

1189

1190

1191

1192

1193

1194

1195

1196

1197

1198

1199

1200

1201

1202

1203

1204

1205

1206

1207

1208

1209

1210

1211

1212

1213

1214

1215

1216

1217

1218

1219

1220

1221

1222

1223

1224

1225

1226

1227

1228

stocks in the superficial sediments of Lake Geneva from 1978 to 1988. *Hydrobiologia*, **207**(1), 161–166. doi: 10.1007/BF00041453

Staehr, P. A., Bade, D., van de Bogert, M. C., Koch, G. R., Williamson, C., Han-
son, P., . . . Kratz, T. (2010). Lake metabolism and the diel oxygen tech-
nique: State of the science. *Limnol. Oceanogr. Methods*, *8*(11), 628–644. doi:
10.4319/lom.2010.8.0628

Staehr, P. A., Christensen, J. P., Batt, R. D., & Read, J. S. (2012). Ecosystem metabolism in a stratified lake. *Limnol. Oceanogr.*, *57*(5), 1317–1330. doi: 10.4319/lo.2012.57.5.1317

Staehr, P. A., & Sand-Jensen, K. (2007). Temporal dynamics and regulation of lake metabolism. *Limnol. Ocean.*, **52**(1), 108–120.

Staehr, P. A., Testa, J. M., Kemp, W. M., Cole, J. J., Sand-Jensen, K., & Smith, S. V. (2012). The metabolism of aquatic ecosystems: History, applications, and future challenges. *Aquat. Sci.*, **74**(1), 15–29. doi: 10.1007/s00027-011-0199-2

Steeman Nielson, E. (1951). The use of radio-active carbon (C14) for measuring organic production in the sea. *Nature*, **167**, 684–685.

Steinsberger, T., Schwefel, R., Wüest, A., & Müller, B. (2020). Hypolimnetic oxygen depletion rates in deep lakes: Effects of trophic state and organic matter accumulation. *Limnol. Oceanogr.*, **65**(12), 3128–3138. doi: 10.1002/lno.11578

Steinsberger, T., Wüest, A., & Müller, B. (2021). Net ecosystem production of lakes estimate from hypolimnetic organic carbon sinks. *Water Resour. Res.*, this issue.

Tadonléké, R. D., Lazzarotto, J., Anneville, O., & Druart, J.-C. (2009). Phytoplankton productivity increased in Lake Geneva despite phosphorus loading reduction. *J. Plankton Res.*, **31**(10), 1179–1194. doi: 10.1093/plankt/fbp063

Tobias, C. R., Böhlke, J. K., & Harvey, J. W. (2007). The oxygen-18 isotope approach for measuring aquatic metabolism in high-productivity waters. *Limnol. Oceanogr.*, **52**(4), 1439–1453. doi: 10.4319/lo.2007.52.4.1439

Vachon, D., & Prairie, Y. T. (2013). The ecosystem size and shape dependence of gas transfer velocity versus wind speed relationships in lakes. *Can. J. Fish. Aquat. Sci.*, **70**(12), 1757–1764. doi: 10.1139/cjfas-2013-0241

Valerio, G., Pilotti, M., Lau, M. P., & Hupfer, M. (2019). Oxycline oscillations induced by internal waves in deep Lake Iseo. *Hydrol. Earth Syst. Sci.*, **23**(3), 1763–1777. doi: 10.5194/hess-23-1763-2019

Wanninkhof, R. H. (1992). Relationship between wind speed and gas exchange. *J. Geophys. Res.*, **97**(92), 7373–7382. doi: 10.1029/92JC00188

Welch, P. D. (1967, jun). The Use of Fast Fourier Transform for the Estimation of Power Spectra: A Method Based on Time Averaging Over Short, Modified Periodograms. *IEEE Trans. Audio Electroacoust.*, **15**(2), 70–73. doi: 10.1109/TAU.1967.1161901

Williams, P. J. L. B., Quay, P. D., Westberry, T. K., & Behrenfeld, M. J. (2013). The oligotrophic ocean is autotrophic. *Ann. Rev. Mar. Sci.*, **5**(1), 535–549. doi: 10.1146/annurev-marine-121211-172335

Woolf, D. K. (1997). Bubbles and their role in gas exchange. In *Sea surf. glob. chang.* (pp. 173–206). Cambridge University Press. doi: 10.1017/CBO9780511525025.007

Woolway, R. I., Jones, I. D., Hamilton, D. P., Maberly, S. C., Muraoka, K., Read, J. S., ... Winslow, L. A. (2015). Automated calculation of surface energy fluxes with high-frequency lake buoy data. *Environ. Model. Softw.*, **70**, 191–198. doi: 10.1016/j.envsoft.2015.04.013

Zeng, X., Zhao, M., & Dickinson, R. E. (1998). Intercomparison of bulk aerodynamic algorithms for the computation of sea surface fluxes using TOGA COARE and TAO data. *J. Clim.*, **11**(10), 2628–2644. doi: 10.1175/1520-0442(1998)011[2628:IOBAAF]2.0.CO;2

Accepted Article
<figure><img src="image_6.png" /><figcaption>**Figure 1.** Hourly time-series of mooring measurements of temperature (A, T), dissolved oxygen (B, DO) and dissolved oxygen anomaly with respect to concentration at saturation (C, DO – DOsat). Hourly mixed layer depth (calculated as the depth where temperature exceeds the value at 0.5 m by 0.2 °C) is shown as a black line. The gray line in panel C represents the zero saturation isoline. Power spectra of the DO as a function of sensor depth are shown in D. The spectrum of the time-integral of air-lake oxygen flux (black line) is shown for comparison. The slanted black solid line represents the characteristic slope of the internal wave continuum spectrum ~ f⁻². The periods of the Kelvin and inertial waves are indicated as vertical lines. The shaded area indicates periods around the diel cycle (18-28 hours).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 2.** Decimal logarithm of vertical turbulent diffusivity ($K_T$) diagnosed from the heat budget (A), Rhône River intrusion inflow (Q) (B) and divergence terms in the dissolved oxygen budget (Eq. 5): concentration rate of change (C), vertical turbulent diffusion plus mixed layer entrainment (D), Rhône intrusion (i.e. horizontal convergence, E), and basin-scale hydraulic uplift induced by the river intrusion (F). The solid black line represents the mixed layer depth and the white lines indicate the zero isolines.</figcaption></figure>
<figure><img src="image_8.png" /></figure>

**Figure 3.** Net ecosystem production (NEP) of oxygen diagnosed from the oxygen budget (A), depth integrated terms (mmol m⁻² d⁻¹) of the oxygen budget in the production-dominated layer ($z < z_c = 8.75$ m, B) and in the respiration layer ($z_c < z < 30$ m, C). Displayed terms are net ecosystem production (NEP, green), rate of change of DO concentration ($\partial DO/\partial t$, black), turbulent diffusion through $z_c$ (Diffusion/Diff. above, blue solid), diffusion through 30 m (Diff. below, blue dashed, only in B), the net effect of Rhône intrusion and uplift (red), mixed layer entrainment (brown), and air-lake exchange (purple). Mixed layer (solid) and photic zone (dashed, 1% incident light) depths are shown in panel A.

Accepted Article
<figure><img src="image_9.png" /><figcaption>**Figure 4.** Gross primary production (GPP), respiration (R) and net ecosystem production (NEP) estimates computed with the scale-separation method (seasonal DO budget plus spectral diel method): A. integrated values in the depth range of 0–30 m, B. depth resolved GPP values, and C. time-mean profiles. Mixed layer (solid) and photic zone (dashed, 1% incident light) depths are shown in panel B. The percentage of segments where the diel cycle was detected according to the prescribed criteria (i.e. "good" segments) is shown as gray steps in C.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 5.** Implementation example of the fitting procedure to remove the contribution of internal displacements to the DO signal. Data is for 22 April at 20 m depth. A. quadratic fit of DO' against T', B. temperature time-series, C. raw (solid) and corrected (dashed) DO time series.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 6.** Determination coefficient ($R^2$) of the daily second-order polynomial fit of DO' against $T'$ (Eq. 4).</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 7.** Time-series of depth-resolved (A) gross primary production (GPP), (B) respiration (R), and (C) net ecosystem production (NEP) estimates computed with the time-domain diel method. Weekly averaged mixed layer and photic zone depths are shown as a black solid and dashed lines, respectively.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 8.** A. Weekly mean of depth-integrated (0–30 m) and B. time-mean profiles of GPP, R and NEP derived from the time-domain diel method. The percentage of segments where the diel cycle was detected according to the prescribed criteria (i.e. "good" segments) is shown as gray steps in B. Results from the scale-separation method are shown as dashed lines for comparison.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 9.** Phase difference (in hours) between the oxygen and light signals at the diel frequency range computed from the co-spectrum of both variables. A positive phase difference means that the oxygen signal lags behind the light signal, as expected if oxygen accumulation is due to photosynthesis.</figcaption></figure>
<table><caption>**Table 1.** Air-lake oxygen fluxes and diagnosed depth-integrated NEP and GPP using different parameterizations for the piston velocity, $k_w$. Mean ± standard deviation for the sampling period (19 April to 4 October 2019) are reported. NEP and GPP are integrated in the production layer, $z < z_c = 8.75$ m, where $z_c$ is the compensation depth, and down to 30 m, respectively. Mean and standard deviation of the different models are shown in the last row. Relative errors were calculated as the ratio of the standard deviation and the mean values for all the different parameterizations.</caption><thead><tr><th>Parameterization</th><th>k<sub>w</sub><br>[cm h<sup>-1</sup>]</th><th>F<sub>air-lake</sub><br>[mmol m<sup>-2</sup> d<sup>-1</sup>]</th><th>NEP<br>[mmol m<sup>-2</sup> d<sup>-1</sup>]</th><th>GPP<br>[mmol m<sup>-2</sup> d<sup>-1</sup>]</th></tr></thead><tbody><tr><td>Soloviev et al. (2007)</td><td>5.8 ± 5.9</td><td>65 ± 22</td><td>69 ± 18</td><td>236 ± 120</td></tr><tr><td>Read et al. (2012)</td><td>4.9 ± 4.5</td><td>55 ± 19</td><td>60 ± 15</td><td>227 ± 116</td></tr><tr><td>Cole and Caraco (1998)</td><td>3.7 ± 2.2</td><td>44 ± 17</td><td>49 ± 16</td><td>217 ± 117</td></tr><tr><td>Crusius and Wanninkhof (2003)</td><td>5.1 ± 9.5</td><td>52 ± 19</td><td>57 ± 18</td><td>224 ± 116</td></tr><tr><td>Wanninkhof (1992)</td><td>3.9 ± 6.8</td><td>41 ± 15</td><td>46 ± 15</td><td>212 ± 180</td></tr><tr><td>Vachon and Prairie (2013)</td><td>9.8 ± 5.8</td><td>114 ± 43</td><td>119 ± 39</td><td>284 ± 131</td></tr><tr><td>Mean</td><td>5.5 ± 2.0</td><td>62 ± 25</td><td>67 ± 25</td><td>233 ± 24</td></tr><tr><td>(Excluding Vachon</td><td>4.7 ± 0.8</td><td>51 ± 8</td><td>56 ± 9</td><td>223.0 ± 8.2)</td></tr><tr><td>Relative error</td><td>36.4%(16.6%)</td><td>39.6%(16.3%)</td><td>36.7%(15.0%)</td><td>10.2%(3.7%)</td></tr></tbody></table>
<table><caption>**Table 2.** Sensitivity of time-mean depth-integrated GPP rates (mmol m⁻² d⁻¹) calculated with the spectral and time-domain methods to the steps performed in order to remove the physical signal (Methods and Supplementary Information). The % reduction with respect to the previous step is given in brackets. The percentage of segments that were discarded according to the different criteria is also indicated. Mean daily GPP values and standard deviation are given in both cases. For the time domain method the standard deviation is calculated from weekly GPP means, in order to be comparable with the spectral method. R²max is threshold for the determination coefficient of the fit between DO and temperature anomalies. If R² exceeded the critical values we considered that physical process dominated DO variations at daily scale and that the diel method could not be applied. The selected methods, for which the main results are shown in the figures and described in the manuscript, are marked with a *. The results indicated with a † correspond to the sensitivity test where the temperature-correlation correction (step i in the frequency domain) was not implemented for the two upper-most sensors located at 0.5 and 2.5 m depth. In this case, the percentage of change indicated in brackets corresponds to a comparison with the base case indicated with a *.</caption><thead><tr><th>Correction</th><th>GPP<br>(0–30 m)</th><th>GPP<br>(0–8.75 m)</th><th>GPP<br>(8.75–30 m)</th><th>%<br>discarded</th></tr></thead><tbody><tr><td colspan="5" style="text-align: center; font-style: italic;">Spectral diel method</td></tr><tr><td>None</td><td>601 ± 242</td><td>304 ± 181</td><td>297 ± 104</td><td>0%</td></tr><tr><td>(i)</td><td>424 ± 174 (–29%)</td><td>237 ± 143 (–22%)</td><td>187 ± 59 (–45%)</td><td>0%</td></tr><tr><td>(i,ii)</td><td>287 ± 130 (–32%)</td><td>196 ± 120 (–18%)</td><td>92 ± 49 (–51%)</td><td>28%</td></tr><tr><td>(i,ii,iii)</td><td>263 ± 124 (–8%)</td><td>190 ± 118 (–3%)</td><td>73 ± 49 (–20%)</td><td>48%</td></tr><tr><td>(i,ii,iii,iv)*</td><td>236 ± 118 (–10%)</td><td>184 ± 114 (–3%)</td><td>52 ± 47 (–29%)</td><td>58%</td></tr><tr><td>(i,ii,iii,iv)†</td><td>250 ± 139 (+6%)</td><td>197 ± 137 (+7%)</td><td>52 ± 47 (0%)</td><td>59%</td></tr><tr><td colspan="5" style="text-align: center; font-style: italic;">Time domain</td></tr><tr><td>No filter</td><td>476 ± 160</td><td>224 ± 123</td><td>251 ± 77</td><td>40%</td></tr><tr><td>Filter, R<sup>2</sup><sub>max</sub> = 1.00</td><td>326 ± 120 (–31%)</td><td>178 ± 49 (–20%)</td><td>148 ± 41 (–41%)</td><td>44%</td></tr><tr><td>Filter, R<sup>2</sup><sub>max</sub> = 0.75*</td><td>234 ± 100 (–28%)</td><td>154 ± 80 (–14%)</td><td>80 ± 35 (–45%)</td><td>68%</td></tr><tr><td>Filter, R<sup>2</sup><sub>max</sub> = 0.50</td><td>193 ± 81 (–17%)</td><td>133 ± 65 (–13%)</td><td>61 ± 29 (–24%)</td><td>80%</td></tr><tr><td>Filter<sup>†</sup>, R<sup>2</sup><sub>max</sub> = 0.75*</td><td>239 ± 97 (+2%)</td><td>159 ± 77 (3 + %)</td><td>80 ± 35 (–45%)</td><td>68%</td></tr></tbody></table>