Solute transport is modelled by a diffusion equation, whereas the presence of density instabilities is modelled through homogenizing the unstable zone.

The vertical diffusion of a passive tracer at concentration C in the water column follows the one-dimensional equation:

Diapycnal diffusivity (which means diffusion through layers of the same density) is calculated by employing the method described below. The mass density of water is deduced from temperature and salinity, which was recalculated from conductivity at 25 °C (Aeschbach-Hertig et al., 2002). The Brunt-Väisälä frequency (N) is the frequency at which a vertically displaced parcel oscillates within a statically stable environment (gravity wave). The squared Brunt-Väisälä frequency (N²) equals:

In the vicinity of the mixolimnion–monimolimnion interface, i.e. at depths between 55 and 60 m, the salinity guarantees the water column stability by countering the thermal gradient due to heating by the geothermal flux. N² varies from 10⁻⁴ s⁻² at the chemocline to 10⁻⁸ in the hypolimnetic part of the mixolimnion (depending on the size of the box used to average temperature and salinity gradients and on the inter- and intra-annual variability). Values of N² calculated for two independent surveys in 2006 and 2007 are shown in Fig. 1. The interannual variability of stratification in the mixolimnion appears to be very high.

If Eq. (2) gives a positive result, the vertical turbulent diffusion coefficient can be calculated with Osborn's Law (1980):

In the case where the right member of Eq. (2) is negative or exceeds a threshold corresponding to that of convective instability (i.e. K_z > 10⁻² m²·s⁻¹), the tracer content of the water layers is homogenized respecting mass conservation. In 2007, this threshold was exceeded at several depths in the water column.

In the model, the heat and solute diffusion coefficients are assumed to be equal. In fact, K_z values in the mixolimnion imply that turbulence is the main phenomenon that leads to mass and heat transport. Close to the chemocline, K_z values are close to heat and molecular diffusivity (sometimes even below heat diffusivity). Differences in diffusivity of solute and heat are supposed to lead to detectable double-diffusive staircases in the density profiles, when the turbulent diffusivity is very low. These staircases are the result of mixing zones (isodensity) alternating with stratified zones. The calculation of the ratio between saline and temperature contributions to density confirms the possibility to observe double diffusive staircases close to the chemocline. Moreover, winter mixing increases turbulent mixing close to the chemocline and prevents the development of staircases. Equation (1) is then solved by implementing a finite differences method with Δ_z = 25 cm and Δ_t = 0.3 days. The sublacustrine spring flow modifies the temperature and conductivity gradients; therefore the diapycnal diffusivity coefficient needs to be calculated at each time step.

The assumptions on the initial conditions and boundary conditions are presented together with the simulation results, because the observations determined the values of input parameters for the model.

Temperature drops on the profiles were observed several times during 2007, while no such drops could be detected on profiles acquired throughout 2006. Detection of cooling events varies from one year to another confirming previous observations (Aeschbach-Hertig et al., 2002). Drops are combined with mixing phenomena, as reflected by locally fluctuating temperature measurements between 40 and 50 m (on the profile dated June 27, 2007). Fig. 2 shows the appearance of cold thermal drops on the SCAMP profiles for May and June 2007. Thermal drops of lower amplitude were also detected also in April 2007 which may be attributed to the same external cause. As April water temperatures are lower than those in May water temperature in the mixolimnion, these drops are less pronounced. The anomalies range between 40 and 55 m and may reach 0.05 °C. At this stage, two hypotheses can be formulated: either the injection of cold water could not be detected in April 2007 as the temperature of the anomaly is similar to the mixolimnion one at this time of year, or there is no anomaly at all, meaning that such inflow is occurring intermittently at these depths.

The temperature “anomaly” was estimated by ascribing the deviation on a calculated theoretically “simply-diffusive” temperature profile between monimolimnion and mixolimnion. In April 2007, the temperature “anomaly” was estimated by ascribing the deviation to a large scale sliding averaged temperature profile (Fig. 2). For the following months, the April “simply diffusive” profile was modified by linearly modifying temperature boundary conditions between April and June 2007 (at 30 and 70 m) and applying the heat diffusion equation to the temperature using K_z calculated from Osborn's formula (see Materials and methods). Even if the deviation from the “simply-diffusive profile” cumulates different external influences (especially in April 2007), it is one method of quantifying the perturbation. This drop was −0.027 °C at 46 m at the end of May 2007 and −0.056 °C at 50 m one month later. It is worthwhile to note that no conductivity anomaly was detected by SCAMP recordings at the centre of the lake. The electrochemical composition of this water inflow thus approximates the mixolimnion composition.

In addition, the LDS temperature analysis between April 25 and July 19, 2007 (Fig. 3A) indicates a bottom-to-top propagation in the water column of a colder water flow, detectable from the 53 m sensor and capable of propagating intermittently and extending to the 25 m sensor. Three major temperature drops are also highlighted (May 14–16, 2007, June 12–20, 2007 and July 11–13, 2007) throughout the period. The beginning and end time of the cooling events are determined through different processing steps. First, the sliding mean of the time series is computed at a large time-scale (10 days) to obtain the trend (T)–see Fig. 3A. Second a sliding mean (called SLIM) is calculated at a small time-scale (with a 1 hour window) to remove the effects of small-scale turbulence. If the difference between the trend (T) and the sliding mean (SLIM) exceeds twice of the standard deviation σ, then it is the starting point of a cooling event. When the SLIM reaches the general trend T, it is the end date of the event. Vertical velocities related to the motion of the parcel of water are computed using the different beginning time of the cooling events at the different depths.

These events represent thermal drops that nearly reach −0.1 °C, the first two of which generate an impact to a depth of 25 m. The most significant temperature decreases were recorded on thermometers placed 50 and 53 m, with a time-lag indicating the bottom-to-top cold water mass propagation.

Fig. 3B shows the second event occurring between June 12 and 20, 2007. Based on these observations, we are in a position to calculate the apparent velocity of the rise of the thermal anomaly within the water column for this particular event. The apparent velocity is calculated by considering the time-lag between the detection of cooling at two different depths. Between 53 and 45 m, this velocity equals 5.7 × 10⁻⁴ m·s⁻¹, whereas between 45 and 25 m, it drops to 2.9 × 10⁻⁴ ms⁻¹. During the third event (July 11–13, 2007), the cold anomaly lasts for a shorter length of time and only extends to a depth of 45 m.

For both of these events, the temperature at 35 m, which represents the minimum profile value, does not exhibit any significant variation as the water mass passes by. This water mass displacement may be ascribed either to an advective effect related to the sublacustrine spring kinetic energy or to a convective effect related to a density difference.

The duration of these cooling cycles serves to distinguish between local mixing phenomena and “overturns” that often appear in slightly-stratified hypolimnions of freshwater lakes. In this case, the resulting temperature variations would affect the water column over a height on the order of a few meters and would last up to several hours. The length T of a static or convective or shear instability event within a stratified medium is inversely proportional to the Brunt-Väisälä frequency, which lies in the range of 10⁻³·s⁻¹ in the Lake Pavin mixolimnion (Assayag et al., 2008). The period of an instability can be approximated by: T ≈ 24 N⁻¹ (Thorpe and Jiang, 1998), i.e. at most about 6 hours in the lower part of the mixolimnion. Since these events were observed over several days, it may be deduced that local instabilities are not the cause of the recorded thermal anomalies.

The sudden temperature drops noticed in the mixolimnion cannot be ascribed to heat fluxes on the lake surface or bottom given that temperature recordings below the thermocline (at 25 m) or above the chemocline (at 58 m) do not reveal any temperature variation. This layer of the lake is thus behaving like an isolated system, at the episode scale and at both its upper and lower boundaries, between 25 and 58 m. The temperature decreases recorded in the water column between these two depths during the studied periods indicate the necessary presence of a cold-water inflow originating from outside the system. Moreover, the minimum temperatures reached during these temperature drops were not recorded elsewhere in the water column and cannot be correlated with internal waves that would have caused intrusions at these depths. The sole hypothesis that could actually be envisaged therefore entails a sublacustrine spring.

The three cooling events correspond to temperature variations around a mean of −0.05 °C, which on average require half a day to form (Fig. 3B). The heating and cooling dissymmetry in the first event deserves our attention, as it demonstrates that cooling is caused by a sudden inflow of cold water whereas heating occurs gradually (Fig. 3B).

An evaluation of the daily heat flux was made between 50 and 55 m, i.e. in the zone of highest recorded temperature variations:

According to the bathymetry established by Delebecque in 1898, the volume of water between 50 and 55 m is 10⁶ m³.

On the first event (Fig. 3B), a temperature drop of 0.07 °C was noticed in one single day. This observation over such a short time enables us to eliminate convective transfers towards the upper levels which occur with a time lag. Accordingly, we have estimated that the spring has a thermal impact of the same magnitude, between 50 and 55 m, over approximately 1/6 of the lake surface. This estimate of the horizontal plume extent was derived from spatial information collected a few years ago during a field campaign of thermal anomaly detection (Assayag et al., 2008). Equation (4) then leads to:

The variation in heat quantity may be ascribed to the cold spring, which has an annual average flow rate, estimated by various methods, at 20 L/s (Aeschbach-Hertig et al., 2002; Assayag et al., 2008), i.e. 1700 m³ per day. The intermittent nature of the observed cooling events reflects the intermittence of the sublacustrine spring. A number of flow rate hypotheses were tested; temperature deviations between the spring and the mixolimnion vary depending on the selected hypothesis from 2.72 °C for a flow rate of 50 L/s, 1.36 °C for a flow rate of 100L/s to 0.68 °C for a 200 L/s flow rate.

Instantaneous flow rates in excess of 200 L/s are not very realistic since they would lead to greater turbulence and high horizontal velocities in the lake.

Finally, an instantaneous flow rate of 100 L/s was selected for the following modelling approach due to its compatibility with both the heat balance and realistic spring temperatures.

In the vicinity of Lake Pavin, the hole “Creux du Soucy”, located at about 2 km from the lake, is a good candidate to explain cool water entries in the lake. The surface water inside the hole has an altitude of 1275 m (25 m higher than the elevation of the lake). The water temperature and conductivity were measured in July 2007 and were respectively 2.1 °C and 23 μS/cm (conductivity at 25 °C is 40 μS/cm). Old legends report that this hole may be directly connected to Lake Pavin by a kind of siphon (Bakalowicz, 1971). This assumption is confirmed by recent diving explorations of the hole, which showed that the hole is connected to a 50 m long quasi-vertical tunnel. The presence of fractured rocks may explain why the tunnel behaves like a siphon. Therefore, the low conductivity of Creux de Soucy water and its location are consistent with the assumption of a sublacustrine input of this water into Lake Pavin. Moreover, the temperatures of surface springs falling into Lake Pavin lie above 5 °C, far above the temperature of the entering water.

The general annual trend in oxygen concentration within the Lake Pavin oxycline, outside of the icy period, is to decrease with time at a given depth (LGE data for years 2006 and 2007). This decrease is explained in particular by the oxidation of various dissolved compounds, such as reduced iron and manganese, in addition to ammonium and methane diffusing slowly from the mixolimnion–monimolimnion interface, consuming oxygen once its concentration begins to increase in the mixolimnion, and precipitating. Precipitates increase the particle content: a rise in particle concentration was observed by means of Doppler acoustic measurements (LADCP). Metal oxides are formed at different oxygen concentration levels, therefore at different depths, which results in observable independent layers with higher particle concentrations between 50 and 60 m (data not provided).

In 2007, periods of oxygen recharge reaching 50 μmol/l were observed at 56 m, even though no known biogeochemical mechanism is able to offer any explanation. This phenomenon was not noticed in 2006.

Fig. 4 A, B shows that these oxygen recharge periods coincide with temperature drops detected between 45 and 56 m. From May 14 to 16, 2007, oxygen recharge appeared to lag behind the cold thermal anomaly. For the second period, June 12–20, 2007, the increase in oxygen also coincided with a slight time lag compared to the temperature anomaly, which might be attributed to anomaly propagation time between temperature and oxygen sensors spaced 150 m apart horizontally. Even if the velocity of anomaly displacement lies on the order of 10⁻³ ms⁻¹ and is consistent with the effect of horizontal turbulent diffusion that can be estimated at 1.5 × 10⁻¹ m²s⁻¹ from the dimensions of Lake Pavin (Stevens et al., 2004), these measurements might also reflect the spatial (horizontal) intermittency or inhomogeneity of the spring water distribution in the hypolimnion. A temperature sensor located close to the oxygen sensor would have been useful to definitely conclude on a correlation between water drop and oxygen recharge events.

At a depth of 56 m, the oxygen recharge periods correspond to oxygen concentration increases of between 40 and 50 μM for an average value of 200 μM at 56 m. The other optode positioned at a depth of 50 m deep displays coincidental yet smaller fluctuations, amounting to 15 μM for an average value of 245 μM over the study period. The sensor at 56 m would thus seem to lie closer to the spring inflow depth.

Three processes may be involved in explaining temperature evolution: turbulent diapycnal diffusion, the initial kinetic energy of the spring (advection), and buoyancy (convection). Based on the velocity of the thermal anomaly rise in the water column, the first and the third hypotheses can be rejected.

For a vertical diffusion K_z equal to 10⁻⁴ m²s⁻¹ (a high value for the Lake Pavin hypolimnion) (Assayag et al., 2008; Bonhomme et al., in prep), the layer thickness influenced by cooling specific to the spring may be estimated by: H=KzT, which yields 3 m in a single day. This velocity is a full order of magnitude lower than that observed for the thermal anomaly rise (∼25 m/day).

We can now compare the advective with the convective hypothesis. The horizontal plume velocity is supposed to be much higher than the vertical plume motion. According to this hypothesis, the homogenisation is very rapid in the horizontal dimension compared with the vertical one and the study of horizontal and vertical motions can be decoupled. After the maximum spatial extent of the plume is reached, the horizontal velocity is no longer taken into account. The situation is simplified in the following way: the maximum spatial extent of the plume is assimilated to a cylindrical jet, whose diameter is the horizontal distance at which temperature differences are detectable. The injection of a fluid is supposed to occur at the lowest depth at which the temperature anomaly is detectable (55 m) and its initial vertical injection velocity is accordingly seen as the lowest vertical measurable velocity of cooling event propagation. On the vertical near the point of fluid ejection, the kinetic effect dominates the transport effect. Floatability becomes dominant beyond a characteristic distance l_m, which is the ratio of movement quantity (M₀) to floatability (Q₀) (Eq. (6)) (Fischer et al., 1979):

In considering a plume of circular cross-section, M₀ is defined as: πD2w02/4, where D is the radius of the plume created by the spring, w₀ the vertical velocity of the fluid measured at the plume origin and centre (assumed located between 50 and 60 m), and Q₀ defined as πD2Δb0w0/4, with Δb₀ representing the floatability flow whose expression is:

Since no salinity difference with respect to the spring can be measured, Δb₀ is to be determined using Morton's formula (Eq. (7)) (Morton et al., 1956), by which one can calculate Δb₀ from remote measurements.

Velocity w was estimated by setting a high value for the velocity of rising thermal drops observed on the LDS recordings, without any prior knowledge of the precise plume center location: w ∼ 5 × 10⁻⁴ m·s⁻¹. This value has to be linked to measured w (see section 3.1). By using a high estimate of the flow rate q ∼ 20 Ls⁻¹ (which is the most likely hypothesis according to (Aeschbach-Hertig et al., 2002)), an entrainment coefficient α equal to 0.1 (usual value in the case of a plume), a plume height H ∼ 25 m (as the plume rises about 25 m in the water column), a value of Δb₀ ∼ 2.3 × 10⁻⁹ m⁴s⁻³ is obtained.

Application of the definition of M₀ with D ∼ 150 m (because the influence of the plume is detectable within a radius of 150 m - see section 3.4), w₀ ∼ 5 × 10⁻⁴ m/s and Δb₀ ∼ 2.3 × 10⁻⁹ m⁴s⁻³ yields a value of l_m ∼ 0.3 m. The accurate estimate of D is not of great importance since the calculation of l_m is not very sensitive to D. This very low value means that floatability dominates vertical displacement in the water column after the first meter.

Note that such a low value for Δb₀ is consistent with non-measurable salinity anomalies correlated with the Lake Pavin mixolimnion spring, since:

The following values of a₀, b₀ and β (a₀ = 0.601, b₀ = 8.45 × 10⁻⁴ and β_s = 0.778 × 10⁻³) were extracted from (Aeschbach-Hertig et al., 2002). Expected conductivity differences are thus on the order of 10⁻² μS/cm, making them non-measurable, which is consistent with observations.

Only the convective transport hypothesis can be used to explain the observed anomaly propagation. The objective here is to focus on the impacts of the sublacustrine spring on mass transfer at the mixolimnion–monimolimnion interface. The modelling tool described above was used to investigate the role of the sublacustrine spring role in maintaining meromixis.

Changes at the mixolimnion–monimolimnion interface were simulated over a 300-day period, which corresponds to the mean unfrozen period on Lake Pavin, in two scenarios:

• • the absence of any spring;
• • the presence of a sublacustrine spring, represented by temperature and conductivity intermittently forced to their initial values on the water column between 45 and 55 m. A consequence of this forcing is to reduce the temperature and conductivity at these depths by countering the diffusive tendency of the monimolimnion towards the mixolimnion.

The initial temperature and conductivity profiles were taken from SCAMP measurements recorded on April 4, 2007. The modelled zone lies between depths of 30 and 80 m. At the boundaries of this zone, the conductivity and temperature readings are forced to their initial values at each time step.

Fluxes at the mixo–monimolimnion interface in the direction of the mixolimnion were evaluated from calculated concentration variations of a fictitious independent passive tracer, with a profile analogous to that of the conductivity and varying between 0 and 1. This independent passive tracer is only transported, does not react with any other tracer in the environment and its concentration does not influence the water density. Matter transfer was determined by calculating variations in the area delimited by the profile between 30 and 60 m, while still incorporating matter losses at zone boundaries.

Heat and dissolved matter both diffuse in the monimolimnion-to-mixolimnion direction. In turn, this diffusion changes K_z. In the scenario of no contribution from sublacustrine water around a 50-m depth, the decrease in heat and salinity gradients through diffusion causes a stability decline at the interface. The presence of water inflow able to maintain a low conductivity near the interface raises the interface stability, while stimulating mixing with the remainder of the mixolimnion shallower than 50 m. This additional water therefore helps insulate the mixolimnion from the monimolimnion by reducing vertical diffusion near the interface. Fig. 5 shows the evolution of a passive tracer concentration after 300 days of simulation in the absence or presence of the spring. The presence of the spring seems to play an important role in the mass transfer between the two compartments of the lake.

The flux evaluation performed with a passive tracer over a one-year period is summarized in Table 1. The presence of a spring would thus help reduce fluxes between the monimo- and mixolimnion by some 20%, with this effect tending to become more pronounced over time. This flux decrease proves significant at the annual scale and serves to limit diffusion above the monimolimnion. On the other hand, the absence of a spring causes a gradual stability decline by diffusion.

The 2005–2006 winter was especially harsh and snowy, accompanied by a very long freezing-over of the lake. Minimum temperatures in the water column during 2006 remained below T_md (the maximum density temperature) calculated (Eklund, 1965) at a depth of 25 m from April to September 2006. In Fig. 6, no thermal anomaly is detected between 50 and 55 m in 2006. Moreover, a comparison of thermal and conductivity profiles between 2006 and 2007 (Fig. 6) indicates that the mixolimnion bottom in 2006 exhibited a higher temperature in 2007. Moreover, the hypolimnion was more stratified in 2006. Observation of LDS data over the months of July and August 2006 confirms this finding (Fig. 7). In comparison with Fig. 3, no element reveals the presence of sublacustrine water intrusion or rising motions at a depth of 50 m during all of 2006 near the lake centre. The sublacustrine spring may have been absent (or had a very low flow) in 2006 or the position of the LDS at the lake center may not have detected the thermal impact of the spring. For example, if the density difference of the spring with mixolimnion environment is stronger in 2006 than in 2007, the rising plume moved faster up the water column. In this case, the minimum observed temperature at a depth of around 30 m throughout all of 2006 might be due to the horizontal plume extension.

During 2007, the mixolimnion was much more homogeneous. The spring periodically added a sizable volume of low-conductivity water. In assuming overflows comparable to what has been estimated in this study, the spring contributed a volume of 4300 m³ in half a day, which is equivalent to nearly 0.1% of the volume existing from the lake bottom to around 60 m.

A seasonal overturn was observed during the 2006–2007 winter thanks to a chain of thermistors in place all winter long underneath the ice. The seasonal overturn of the lake mixolimnion occurred in two stages: in the initial stage, water adjacent to the density maximum at these depths (3.9 °C) (Eklund, 1965) dipped at the end of January to a depth in the range of 50 m. The 10 m layer above the chemocline (between 50 and 60 m) was therefore not mixed (Fig. 8). An inverse stratification then occurred between the surface and 50 m from the bottom in the mixolimnion (data not provided).

The sudden inflow of a large quantity of low-salinity water around February 25, 2007 (Fig. 8) eventually flattened the salinity gradient, which enabled the seasonal overturn to touch the chemocline (at 60 m) during a second stage. A heat balance identical to that depicted in Section 3.4 between 40 and 60 m makes it possible to attribute the temperature drop at these depths to a cold external water inflow and thus undeniably to the presence of a spring.

During the 2006–2007 winter, it thus seems likely that the sublacustrine spring helped fix the mixing depth at the time of spring mixing as well as diffusion at the mixolimnion–monimolimnion interface.

The difference between the 2006 and 2007 profiles (Fig. 6) may be explained by seasonal overturn differences due not only to changing meteorological conditions, but to varying hydrological processes, favouring either the occurrence of water inflow from the spring or an absence of such inflows.