A novel approach for accurate quantification of lake residence time - Lake Kinneret as a case study.

Water residence time, which is affected by increasing water demands and climate change, plays a crucial role in lakes and reservoirs since it influences many natural physical and ecological processes that eventually impact the water quality of the waterbody. Thus, accurate quantification of the water residence time and its distribution is an important tool in lake management. In this study we present a novel approach for assessing the residence time in lakes and reservoirs. The approach is based on the Leslie matrix model that was originally developed for the analysis of age-structured biological population dynamics. In this approach the water in the lake is divided into different age classes each representing the time since the "parcel" of water entered the lake and provides an overall picture of the water age structure. The traditional approach for calculating residence times, which relies only on the lake volume and annual inflow or outflow volumes thereby disregarding any previous information, is very sensitive to large interannual variation. While the proposed approach produces the fraction and volume distribution curves of all age classes within the lake for each simulated timestep. Thus, in addition to mean residence time, the fraction of young water (FYW), quantifying the "young" fraction of water in the lake can be analyzed. The same is true for any other age class of water. The approach was applied to Lake Kinneret (Sea of Galilee) historical data collected over 32 years (1987-2018) and for prediction of long-term time series based on several future scenarios (inflows and outflows). It offers a more accurate quantification of the mean residence time of water in a lake and can easily be adapted to other waterbodies. Comparison of simulation results may serve as basis for determining the lake's management policy, by controlling the inflows and outflows, that will affect both the mean residence time and the fraction of "young/old" age classes of water.

Introduction

Increasing water demands from lakes and reservoirs in arid and semi-arid regions, due to population growth, agricultural irrigation, and economic activities, in addition to changes in climate patterns, result in significant fluctuations in the volume and depth of waterbodies. These affect the physical, hydrological and ecological processes in the aquatic ecosystem (Li et al., 2010). Climate change models for the Mediterranean region predict a 25–30% decrease in precipitation and higher evaporation by the end of the twenty-first century accompanied by even stronger reduction in annual runoff volumes in rural areas (Erol and Randyr, 2012). In addition, occurrence of extreme events is expected to increase (e.g., flooding, extended droughts), which will magnify the seasonal and multiannual water level fluctuations, creating hydrological stresses, such as prolonged hydraulic retention time in lakes (Jeppesen et al., 2015). Thus, understanding lakes and reservoirs’ hydrodynamic processes and characteristics under varying water levels is essential for determination of optimal management policy. However, simple monitoring and flows and fluxes analyses may produce difficulties in gaining valuable information on the ecosystem's functioning (Deleersnijder et al., 2001). Therefore, estimating a specific time-scale analysis variable is often used, such as water age and residence time, both of which have been widely applied to numerous waterbodies (e.g.: Qi et al., 2016; Viero and Defina, 2016; Liu et al., 2012; Li et al., 2010; Monsen et al., 2002; Deleersnijder et al., 2001; Delhez et al., 1999; etc.).

Water age and residence time are conducive variables used for estimating water mass exchange and can also be treated as effective factors for investigating patterns of water quality in lakes and reservoirs. Zimmerman (1976) defines residence time as ‘‘the time it takes for any water parcel of the sample to leave the waterbody through its outlet”, and water age as ‘‘the time it has spent since entering the waterbody through one of the boundaries.’’ Water age is, therefore, defined by spatial heterogeneity; particles at different locations within a water body will have different ages (Monsen et al., 2002). Water age and residence time naturally complement each other and while age traces water according to when it entered the domain of interest, residence time follows the water until it leaves it (Delhez et al., 2014).

Water time-scale variables play a key parameter in limnology affecting biochemical processes, nutrient concentrations and water quality in general. For example, more than 40 years ago, Vollenweider (1976) developed a regression relationship for average planktonic agal chlorophyll concentration as a function of the annual phosphorusload, normalized by waterbody area, mean depth, and hydraulic residence time, for a group of European waterbodies. Several studies have reported on a significant positive correlation between chlorophyll-a concentration in the surface water layer and the time-scale variables, residence time or water age (Zou et al., 2020; Gao et al., 2018; Qi et al., 2016; etc.). Zhang et al. (2010) demonstrated that an increase in time-scale values yielded a decrease in dissolved organic matter, and an increase in water salinity. Maavara et al. (2018) who developed a mechanistic model to predict N2O emmission through nitrification and denitrification in various water bodies, reported that a water residence time kinetically limits the extent of these processes in waterbodies. Other studies investigated the relationship between residence time and the availability of nutrients, such as the significant impact of residence time on the total nitrogen change rates in waterbodies (Tong et al., 2019). Furthermore, increase in cyanobacteria biomass was positively associated with retention time for various lake types, such as humic and clear lakes (Richardson et al., 2018). Recently, Zhao et al. (2022), who studied the effects of water residence time on nitrogen fixation, reported that short lake water residence time can severely inhibit nitrogen fixation capacity by inhibiting the growth of nitrogen fixing cyanobacteria, and hence may serve as management tool for controlling cyanobacteria growth and eutrophication development. These publications demonstrated the effect of time-scale variables on lake ecology and as consequence on lake water quality and may therefore provide a quantitative tool to support the management policy of lakes.

Methods for estimating these time-scale variables (water age and residence time) rely on in situ tracer studies or on different type of models. For example, Kratzer and Biagtan (1997) used dye tracer (Rhodamine) and followed its movement downstream and indicated that the average tracer travel time at a downstream location was higher than the average tracer travel time at an upstream station. In other words, at each site the dye tracer arrived at different times resulting in longer or shorter travel times compared to average time. Other studies have used computerized tools, mainly numerical hydrodynamic models for modeling time-scales, most of them three-dimensional models that included fluid dynamics modules (e.g.: Qi et al., 2016; Li et al., 2010; Zhang et al., 2010; Monsen et al., 2002; Deleersnijder et al., 2001; Delhez et al., 1999). These tools are resource intensive and require preliminary preparation of the model settings in addition to its calibration and verification, skilled users and sufficient computational resources.

In this study we present a novel approach for estimating the water age distribution and residence time in lakes. The approach is based on the Leslie matrix model that was originally developed for the analysis of age-structured biological population dynamics (Leslie, 1945). The approach has been modified and applied to non-steady state reactors and treated wastewater storage reservoirs (Juanico and Friedler, 1994), and can provide information on the hydraulic age distribution of the lake or reservoir.

The Leslie matrix model (Leslie, 1945) is a deterministic model widely used to assess the age structure of populations (e.g., Iwasaki et al., 2010; Bergek et al., 2012 and the references therein; Monte, 2018 and the references therein; Sun et al., 2021). Briefly, this model classifies individuals in a population into age classes. Population density is observed at regularly spaced times and is partitioned into several age classes. The number of individuals in each age class is calculated at each timestep based on: (1) survival rate - the fraction of individuals in the current timestep that survive until the next timestep. These, advance from age class x to age class x+1, (2) effective fertility - the number of offspring born to an individual parent in age class x that survive until the next census, (3) harvesting mechanism that specifies the harvesting policy of the population which can be defined independently for each age class)e.g., how many trees from each age class the forestry industry can fell in a given tree stock(. These can be expressed by two matrices the Leslie Matrix representing the first two transformation laws and the harvesting matrix representing the harvesting policy. Following these steps, the age structure of the population can be represented by a vector at any time.

Water in lakes can be considered as a population divided into different age classes each representing the time since a conceptual “parcel” of water entered the lake. This concept of dividing lake water to different parcels having different transport time, has been addressed in the literature (e.g.: Qi et al., 2016; Monsen et al., 2002; Zimmerman,1976). In this study, we modified the Leslie Matrix model to represent the population dynamics as a water age structure within a non-steady-state flow lake ecosystem. The water age distribution allows us to estimate the mean residence time of water in the lake and to analyze the water age structure at given times. We applied the developed approach to the sub-tropical Lake Kinneret (Sea of Galilee), with data collected over 32 years (1987-2018) and then used the method to evaluate the effect of long-term future management scenarios on the age of the water in the lake.

Methodology

Model description

The main concept

When applying the Leslie Matrix model to a lake, the water volume in the lake is divided into groups based on age classes, where at each timestep the volume and the fraction of each age class is quantified. In this case, survival in the Leslie Matrix model represents water that remains in the lake, i.e., the water that does not evaporate, or flow out from the lake (analogous to natural mortality). In the modified model, water which enters the lake at each timestep plays the role of fertility (“producing” young water in the lake) while water withdrawals are analogous to harvesting.

Description of the modified model

In the Leslie Matrix model, at each timestep the number of living individuals is quantified for each age class. Similarly, the modified model estimates the water volume of each age class by using two indexes; i (1…n) which counts the number of timesteps, and j (1…k) which counts the number of age classes, and the model estimates the water volume in each age class (j) during the time interval (i).

In the original model, in parallel to shifting forward in timesteps, individuals that survive move from one age class to the next, older, age class (“getting older”). Individuals that do not survive, leave the class, and are subtracted from the population. In the modified model, when shifting forward in timesteps a certain amount of water leaves the lake (e.g.: withdrawals, evaporation, river outflow etc.) and the water volume is reduced (analogous to natural mortality and harvesting), the remaining water (“the survivors”) progresses in the age classes and moves to the next class, thereby “getting older" which in essence is equivalent to water residing for a longer time in the lake.

Similar to the original Leslie Matrix model where offspring are added to the population as the youngest age class, in the modified model, inflows to the lake enter the youngest age class. This is because inflows are the newest water added to the lake and are thus considered the youngest class at each timestep. As an example, if the model timestep is one year, then the inflow is defined, at the end of the interval, as age class 1, which is one-year-old; that is to say the water stays in the lake for one year. Assuming completely mixed conditions, water outflow (mortality and harvesting, see above) from each age class is proportional to its fractional volume. Completely mixed conditions can be assumed if the timestep is similar or longer than the lake mixing cycle, i.e., for a monomictic lake, such as Lake Kinneret (Berman et al., 2014), the timestep should be at least one year.

It should be noted that we are fully aware of some possible short-circuiting in the water flow, i.e. water may enter and leave the epilimnion during the stratified period. However, accounting this will require a much more detailed and complex calculation and finer time-steps with much richer data requirements such inflow water density, density of the epilimnion, etc.

The required input data, for the modified model, includes the initial lake volume (Vinitial, scalar; analogous to the total initial population) and water flow data, i.e., water entering (offspring) the lake during the determined timestep and leaving (mortality and harvesting) the lake during this timestep. Therefore, if the timestep is one year the input data for inflows (Qiin (i=1...n)) and outflows (Qiout (i=1...n)), are annual flows.

Model development

The age distribution of the water is constructed in steps, assuming that the lake water is divided into k-age classes. If vi,j = volume of water in age class j at time i; and Fi,j = fractional volume of age class j at time i, then Eq. (1) defines the relationship between the two, where Vitotal is the total volume of the lake at time i Eq. (2).

Combining Eqs. (1) and (2) yields the following expression:

If completely mixed conditions are assumed, outflow water (withdrawals, natural stream outflows and evaporation) from each age class is proportional to its fractional volume.

Where is the outflow volume that leaves age class j at the end of timestep i and Qiout is the total outflow from the lake at time interval i.

Inflows only enter the youngest age class. In our case, it is assumed the inflows enter the lake at the beginning of each time interval and part of it leaves the lake by the end of the time interval. The volume of the lake at the end of each time interval is given by:

Where is the total volume of the lake at time i-1.

Let the mid-interval total lake volume () be defined as the arithmetic mean of the volume at the beginning and at the end of each time interval:

Note: is the total lake volume at the beginning of time interval i, which is equal to the total volume at the end of time interval i-1.

The age classes of water in the lake may be classified into the following two types:

Youngest age class (j = 1) which contains water that entered the lake at time interval i.

Age classes j>1

Youngest age class (j = 1)

Since we assumed that inflows enter the youngest class (j=1) at the beginning of each time interval and part of it leaves the lake by the end of the time interval, the volume of this class is modified over the course of the time interval:

Combination of Eqs. (1),(4) and (7) when using the mean lake volume yields an expression for the fractional volume of age class 1 at the end of time interval i:

Age class j>1

The volume of the j-th class by the end of time interval i can be expressed as:

Where is the volume of age class j-1 at time interval i-1.

Note that moving from time interval i-1 to i, age class j-1 gets older and becomes age class j. Now, similar to (8), combination of (1), (4) and (9) when using the mean lake volume yields an expression for the fractional volume of age class j:

Fig. 1 illustrates the framework of the model, water from age class 1 at timestep 1, which remained in the lake (“survived”), aged by one year and is present at timestep 2 in age class 2. It should be noted that the volume and fraction of this class are smaller than the volume at timestep 1 since some of the volume of the age class left the lake and the total lake volume has changed.

Fig. 1

Simplified scheme of the modified model framework, input and output data

At the end of the simulation two matrixes are obtained, one describes the volume of each age class and the other expresses the age classes’ fractional volumes, for each timestep (Fig. 1). Each vertical vector (column) of these two matrixes represents a timestep while each horizontal vector (row) defines the age of the class. Namely, if the simulation timestep is one year, then each vertical vector is a different year where row 1 describes age class of year one, row 2 represents age class of 2 years and so on.

The fraction matrix obtained, can be used to calculate the mean residence time (τi) based on the age structure for each timestep (each simulation year), as follows:

Where is the fraction of each age class j at timestep i and k is the maximum age class which contains the oldest water in the lake, that is “k years-old” (or older). For example: assuming that three-timesteps are modeled and the fractions of age class 1, 2 and 3 are 0.3, 0.6 and 0.1 respectively. It should be noted that age class 3 contains the oldest water in the lake, which entered the lake three years earlier, and age class 1 contains the youngest water in the lake. Hence, the average age of the water (or the mean residence time) is 1.8 years (0.3 × 1+0.6 × 2+0.1 × 3).

Study site

Lake Kinneret (Sea of Galilee) is a meso-eutrophic lake with mean annual primary production of 650 gC m−2 (Berman et al., 1995), located at an elevation of about -210 m,i.e. 210 m below mean sea level, in the northern part of the Dead Sea Rift Valley (part of the Afro-Syrian Rift Series). The lake is 22 km long and 12 km wide at maximum width; when full, the maximum and mean depths are 44 and 24 m, respectively, and its surface area is 170 km2. Lake Kinneret is a warm, monomictic, lake and due to its importance, its inflows and outflows are continuously monitored by Israel Water Authority (IWA) and by Mekorot, Israel's National Water Company. The lake receives inflows (average of 600 mcm/y; million cubic meters per year), during the years 1987-2018) from 6 main sources including the Jordan River with an average of 390 mcm/y and approximately 70% of the inflows (Gal et al. 2003), direct rain with an average of 66 mcm/y and direct runoff averaging 62 mcm/y. The lake outflow is composed of five sinks, where the major ones are The National Water Carrier, NWC, with an average of 250 mcm/y and evaporation with an average of 240 mcm/y, while only about 10 mcm/y exits the lake through its natural outlet to the Southern Jordan River. Since the mid-1990s the lake ecosystem has undergone a number of significant limnological changes (Zohary, 2004; Gal and Williamson, 2010; Zohary and Ostrovsky, 2011) including five consecutive drought years between the winter of 2013/14 to 2017/8 that exhibited extremely low inflows into the lake and as a result, low outputs, mainly low intake to the NWC, and a declining water level (Tal, 2019; Lachmani, 2019).

Model application

Initial conditions and baseline simulation

As Lake Kinneret is monomictic we assess the lake water age distribution (residence time) using a one-year timestep. The calculations are made for hydrological years which start on 1st October and ends on 30th September the following year. Data used in this study were collected as part of the on-going Lake Kinneret monitoring program (Sukenik et al., 2014) and included annual inflow and outflow data for 32 years over from 1987 through 2018. A lake water volume of 3,755 mcm estimated for October 1, 1987 served as the initial volume.

At the beginning of the simulation the structure of the water age classes in the lake is not known, hence, all water receives an age of one year. During the simulation of the first years, age classes are constructed during each timestep, the age class structure is not stable and contains only a limited number of age classes. Therefore, the model requires a spin-up period. This was achieved by repeated simulation of the year 1987 for 100 times at equilibrium conditions (inflow=outflow), so the lake volume on October 1st did not change between years. The structure of the water age classes in the lake obtained after 100 timesteps was used as an initial condition for the simulation, and then calculated for each year over a period of 32 years (1987-2018).

Future scenarios

Israel Water Authority (IWA), the body that manages the lake, developed a strategic plan for protecting the lake ecosystem in coming decades (2018-2050), in order to cope with extended drought periods (IWA, 2018; Tal, 2019). The strategic plan followed an unprecedented five-year drought in northern Israel, which resulted in decreasing flows in the streams feeding the lake, to the lowest ever recorded (Lachmani et al., 2019) and a lake level closest to the lowest ever recorded. Under the assumption that climate change will continue (Gal et al., 2020) and intensify in the lake region (Tal, 2019), the strategic plan includes importing excess desalinated water, from the national water supply system, into the lake. To evaluate the potential impact of climate change on inflows into the lake and water demand from the lake, thirty different long-term time series, predicting the inflow and outflows entering and leaving the lake for the years 2020-2050, were developed by IWA (based on RCP 4.5 and RCP 8.5; IWA – Israel Water Authority, 2018a) . The long-term time series realizations referred to the predicted regional water demand, in addition to climate change effects expressed by elevated air temperature, reduced precipitation and lower inflows to and outflows from the lake. While all realizations had the same characteristics, they were shifted by one-year resulting in a displacement in the timing of peaks and minima. In addition, there was an annual inflow decrease of 10 mcm/y over the course of the simulation period. The average predicted overall inflows and outflows for 2020-2060 were 515 mcm/y for both inflow and outflow which is approximetly 100 mcm/y less in comparison to 2000s. The strategic plan prescribed that an additional annual volume of imported water that will be added to the lake during a certain year will be utilized during the same year. Importing water to the lake will be executed gradually, starting with an addition of about 10 mcm/y in 2023 and reaching its maximum in 2050. This was planned under two possible scenarios (Table 1) that included, (1) importing the excess water directly into the lake (LAKE scenario), reaching a maximum 304 mcm/y in 2050, and averaging 184 mcm/y over the entire period (2030-2050), and, (2) adding part of the imported water directly into the lake and part to streams in the upper part of the lake's watershed (WATERSHED scenario) a portion of which will drain through the natural drainage system into the lake. In this case, water addition to the lake will reach a maximum 245 mcm/y by 2050 and average 129 mcm/y over the whole period. The two scenarios as well as “business as usual” (BAU) scenario, in which no water is imported into the lake from external sources, were compared by assessing the age distribution of the water within the lake over the period of the realizations.

Table 1

Amount of imported water entering the lake at both possible scenarios during the examined period of the strategic plan (2020-2060).

	Amount of imported water (mcm/y)
Year	BAU scenarioNo water is imported to the lake	LAKE scenario Importing water directly into the lake	WATERSHED scenario Importing water via the lake's watershed
2025	0	21	19
2030	0	67	47
2035	0	96	74
2040	0	122	94
2045	0	153	118
2050	0	184	129
2051-60	0	184	129

Thirty long-term time series forecasting inflow and outflow to the lake during the years 2020-2050 were established by IWA for each scenario. The long- term time series forecasting were extended to year 2060 in order to study a longer period and to analyze the effect of water import at the end of the period. These data served as input for the three scenarios, creating a total of 90 model simulations.

Model output analysis

The developed model was initially applied to data collected from Lake Kinneret during the years 1987-2018 than it was used for analyzing future scenarios. We then will focus on two very different years including a high (2013) and low inflow year (2018). While 2013 was a rainy year that experienced large inflows into the lake (776 mcm/y; Lachmani et al., 2019), 2018 was a dry year following 5 consecutive drought years (2014-2018), with inflows between 312-390 mcm/y (Lachmani et al., 2019).

Results

Model implementation for the years 1987-2018

Analysis of the age structure in the lake 1987-2018

The age structure obtained by the water age distribution approach, shows, as expected, that for each examined year, some of the water had resided in the lake for a long time, while some was younger. The fraction of young water (FYW) is a quantification of the “young” fraction of water in the lake (modified from Juanico and Friedler, 1994). For example: FYW10 is the fractional volume of water that has resided 1-10 years in the lake, expressed as a proportion of total lake volume (Fig. 2a).

Fig. 2

Age distribution of water in the lake: (a) fraction of water volume in the year 2000, the marked area is the FTW10 for that year, and, (b) the cumulative fractions of young water (FYWn) over the simulated period (1987-2018).

Analysis of FYWs for varying ages highlights the differences in response to dry and wet years. As expected, when inflows were higher, the FYWs increased and vice versa during dry years (Fig. 2b). During the 1992 hydrological year, for example, inflows and outflows were very high (1,526 and 980 mcm/y, respectively), resulting in FYW2 =36%, and FYW15=∼100% (i.e., practically all water in the lake was younger than 15 years). On the other hand, in 2001 inflows and outflows were much lower with values of 558 mcm/y and 468 mcm/y representing approximetly 30% and 50% of the outflows in 1992, respectively. As a consequence, in 2001 only 15% of the water was ≤2 years old, i.e. 85% older than 2 years old, and as little as 18% of the water was older than 15 years.

Mean residence time 1987-2018

For each simulated period, the mean residence time was calculated based on the age structure approach (Eq. (11)). The mean residence time of the water in the lake increased by 24% over the period 1987-2018 from 6.6 years in 1987 to 8.6 years in 2018 (Fig. 3), with linear slope of 0.05 y−1 (R2=0.67; p<0.005). Over the course of this period, inflows to- and outflows from- the lake significantly changed and affected the mean residence time which increased and decreased inversely to the flows. For example, in 1992 the inflow and the outflow were very high, therefore the mean residence time was low (6.8 years), on the other hand, during the dry 2001 the mean residence time was higher (7.3 years). From 2013 onwards the mean residence time steadily increased as a result of five consecutive dry years due to low inflows and as a result minimal water extraction. As a consequence, the linear increase in mean residence time calculated for the 2013-2018 drought period was 0.3 y−1 (R2=0.98; p<0.005).

Fig. 3

Annual inflow (Qin), outflow (Qout), mean residence time and for comparison two estimates of the traditional residence time based on the ratio between lake volume and inflows (V/Qin) and lake volume to outflows (V/Qout).

Calculation of the traditional residence time (lake volume divided by the total inflow or outflow; Wetzel, 2001) for the simulated years inversely followed the inflows and yielded a much more notable variation between the simulated years, compared to mean residence time calculated according to the age distribution approach. The increase in residence time during the 1987-2018 period, based on the traditional calculation, was from 4.2 (1987) to 10.9 y (2018) or from 6.2 (1987) to 9.5 y (2018), depending if inflows or outflows were used, respectively. In these cases, the change represents at least a 259% and 152% increase in residence time over that period based on inflows and outflows, respectively.

2013 and 2018

The high inflows of the year 2013 entered age class 1 and after subtracting the relative outflows leaving the class during the time interval, age class 1 constituted about 15% of the total volume of water in the lake (Fig. 4a). For 2013, age class 1 is the largest one followed by an almost exponential reduction of the water volumes of the following age classes. In the case of steady inflows and outflows, an exponential trend is expected, since during every timestep of the simulation, the volumes of the age classes decrease (due to outflows), therefore, an “older” class would become smaller and smaller (for more details see Juanico and Friedler, 1994). On the other hand, in 2018, the age class volumes distribution was not as orderly as in 2013 as it followed five consecutive drought years (2014-2018), with low inflows into the lake. These low inflows are well noticed (Fig. 4b) and are expressed in age classes 1-5, 2018 to 2014, respectively, especially in relation to the relative fraction of the same age classes in 2013 (Fig. 4a). The high inflow of 2013 is clearly pronounced as age class 6 in 2018. Though, by 2018 the class 6 containing the 2013 inflows was smaller than when it first entered the lake (a total decline of 249 mcm) it was nevertheless the largest class representing 10% of the volume of water in the lake, as opposed to only 7% of age class 1, 2018 inflows.

Fig. 4

Model output for water entering the lake in 2013 (a) and 2018 (b). Age class volume (mcm; left Y-axes) and Age class fraction (right Y-axes) vs. age class

Model prediction for the years 2020-2060

As aforementioned (Section 2.4 above) the water age distribution approach, developed in this study was applied to proposed water importing scenarios, developed by the IWA, using 30 long-term time series forecasting of inflow and outflow to, and from, the lake during the projection period. An example for the input and output of one time series is shown in Fig. 5, for all three scenarios. The results show that the BAU scenario is substantially different than the two other scenarios: starting from 2030 inflows and outflows are typically much lower than the ones projected for the two water import scenarios, with an average difference of about 100 MCM/y; consequently, a higher mean residence time is predicted for the BAU scenario. For the BAU scenario, the water exchange is low, and the mean residence time considerably increases with time, reaching 12.6 y in 2060 compared to 8.0 y in 2000. Moreover, the gap between its value and the mean residence time for the two other scenarios increases over time. The maximum difference between the mean residence time predicted for the BAU and the LAKE and WATERSHED scenarios, in 2060, was 3.4 and 2.7 y, respectively. The difference between the LAKE and WATRERSHED scenarios is less pronounced, but consistently increases over time from no difference in 2030 to a difference of 0.25 y (3% difference) in 2040 and a maximum difference of 0.8 y (10% difference) in 2060.

Fig. 5

Simulated annual inflow, Qin, (A) outflow, Qout, (B) and mean residence time for the three scenarios (C).

For each water import scenario (BAU, LAKE, WATERSHED), the model was executed over 30 time series and thus predicted 30 solutions for each simulated year, creating three ranges of possible results, one for each scenario (Fig. 6a). The increase in the mean residence time, over the simulated period, for the BAU scenario is obvious in all 30 predicted time series,rangingfrom 8.0-8.3 y in 2020 to 10.3-13.0 y in 2060 (the end of the forecast period). As water is planned to be imported starting from 2030, the mean residence times of all three scenarios overlap until 2031 and then deviate from each other. The maximum difference in mean residence times, between the examined scenarios occurred at the end of the predicted period (2060), with 5 and 4 year differences between BAU and LAKE and BAU and WATERSHED scenarios, respectively. These gaps over the period 2032-2060 were found to be statistically significant (ANCOVA, p<0.05).

Fig. 6

Range of all 30 average residence time solutions (for 30 time series) for each of the three scenarios (A), and FYW10, the fraction of water residing in the lake 10 years or less (B) obtained from the various scenarios. The median results for each scenario are depicted by the solid lines. The shadowed areas describe overlap of the solutions obtained for the three scenarios

The FYW10 in 2020 for all three alternatives ranged between 67% to 86%, meaning that about 13%-34% of the water resided in the lake longer than 10 years (Fig. 6b). Since the forecast predicts a long-term decrease in inflows and outflows throughout the period, the obtained FWY10 also decreased, signifying a larger volume of water residing in the lake longer than 10 years. For example, the FYW10 in 2040 ranged between 54 to 70% for the BAU scenario, and being higher for the other two scenarios, with 63-76% and 65-79% for WATERSHED and LAKE scenarios, respectively. In 2060 the FYW10 decreased by 5%-6% for the BAU alternative, while a more moderate reduction was obtained for the WATERSHED (3%-4%) and LAKE (up to 3%) scenarios (compared to 2040).

Analysis of the solutions for the LAKE and WATERSHED scenarios indicated similarity in the mean residence times from 2020 to 2038. From 2039 a statistically significant difference is observed between the mean residence times of these two scenarios (ANCOVA, p<0.05). Despite the relatively small differences in annual inflows in the two scenarios (average difference of about 40 mcm/y), the cumulative impact of the differences over time results is a statistically significant difference of up to 2 y in the mean residence time in 2060 (∼25% difference).

Discussion

The novel approach presented here for estimating the lake mean residence time, accounts for the relative annual contribution of all inputs and outputs of water into a lake. The results show that flows changed significantly over the simulated period, while changes to the mean residence time were much less pronounced and as expected its value slightly increased during dry years with low flows, and vice versa. The traditional residence time, calculated by the ration of lake volume to total inflow or outflow, is a well-known method in limnology (Wetzel, 2001). It is, however, a simplification of the complex calculation of the real residence time, since it re-calculates the residence time each year without considering previous years, i.e., the calculation does not account for the residence time at the beginning of the year. Thus, it can be suitable only for lakes and reservoirs operating under steady or quasi-steady state, where inflows and outflows do not differ significantly from one year to another. So, with the traditional calculation, an extremely dry year would provide an estimated very long residence time for a lake though the previous years may have been extremely wet with very high inflows. In contrast, the developed mean residence time calculation includes a “memory” component as it depends on past data of the examined period, and therefore provides a more realistic and accurate solution. For example, calculating the traditional residence time, based on inflow, for the consecutive years 1991 and 1992 yielded: 9 y and 3 y respectively (Fig. 2), a gap of 6 years caused by the high inflows during 1992 that had inflows and outflows of 1,526 and 980 mcm/y, respectively, compared to 1991 which had inflows and outflows of 430 mcm/y. The traditional-theoretical residence time in 1992 does not consider the effect of the preceding three-year drought, and relies only on the high inflow during 1992 and therefore it dramatically decreases in 1992. On the other hand, the water age distribution approach which accounts for conditions through previous years yielded quite similar values of 7.6 and 7.7 y for 1991 and 1992, respectively.

The water age distribution approach yielded the definition of the fraction of young water, FYW, where FYWs increased with increasing flows and vice versa, decreased during dry years. The significant change in inflows during 1991-1992 is well noticed in the calculated FYWs (Fig. 2), therefore for a better understanding of the overall picture, both parameters should be analyzed; mean residence time and FYWs.

Our approach provides an overall picture of the water age classes distribution obtained for each year of the simulation, as was illustrated for 2013 and 2018 (Fig. 4). These data allowed us to compare between years and to draw conclusions regarding age classes and the residence time of the water in the lake. For example, in 2013, 50% of the water resided in the lake five years or less (obtained by summing the fraction of age classes 1-5). In contrast, in 2018 only 28% of the water resided in the lake 5 years or less, while 72% had been in the lake longer than 5 years, namely, the water in the lake was much “younger” during 2013 compared to 2018. When inflows and outflows are low the water mass exchange is poor and the water resides in the lake much longer than during years with high inflows and outflows, when higher water mass exchange generates a high proportion of younger water age classes (Fig. 2b).

The mean residence time calculated by the new approach for 1987-2018 was less sensitive to fluctuations in the annual inflows and outflows than the FYW. This is due to the fact that the mean residence time depends on the complete hydraulic age distribution of the water in the lake and thus has inertia that smooths annual variability and generally better reflects long-term trends. In contrast, FYW of n years depends only on changes in inflows and outflows during the last n years and is therefore more sensitive to large interannual variation resulting in a larger variability than the mean residence time. The smaller the n of FYW the higher the variability. Hence, considering both terms, mean residence time and FYWn, is essential for obtaining an overall assessment of the state of water in lakes.

Dividing the waterbodies to different parcels of water having different transport time, is well established in many hydrological models (e.g.: McGuire and McDonnell, 2006 and ref. therein) some of them are very complex and interpreting their results may be difficult. Our developed approach may seem like a simplistic method which is based on a numer of assumptions such as assuming completely mixed conditions and disregarding stratification. However, as we use an annual time-step, we assume that stratification may occur during part of the year but complete mixing does prevail throughout the rest of the year. Furthermore, the approach combines a number of time-scale parameters: average residence time, age group distribution and FYW, and it relies on past data of the entire examined period and therefore may provide a broad picture of the lake water exchange. In addition, the developted approach would not only allow to assess water residence time through water age distribution, but it may also be adjusted and applyed for estimating nutrients residence time, by calculating annual loads and corresponding dynamic annual budgets. This informationis is expected to attract the attention of water managers (e.g.: in the context of eutrophication).

It should be noted that for lakes the fraction of young water is significant and provides an indication of lake water exchange which influences many natural physical and ecological processes, hence indirectly expresses the lake water quality (Li et al., 2010). In contrast, in treated wastewater storage reservoirs, which are a seasonal storage, and hence the timescale is days unlike years in lakes, the fraction of old water is much more important as it indicates longer time for natural treatment processes to occur which improve effluent quality. In other words, for treated wastewater storage reservoir a large fraction of old water is desirable (Juanico and Friedler, 1994).

A similar approach of calculating both mean residence time and the fraction of young water, was previously applied for treated wastewater storage reservoirs. These variables were used as a design criteria and operation rules for treated wastewater storage reservoirs in Eastern Sicily, Italy (Consoli et al., 2011; Cirelli et al., 2008). The authors indicated that the percentage of fresh effluents correlates with poor performance of the reactor and enables estimation of the removal of pollutants in non-steady-state systems, in addition to forecasting the quality of the effluent released from the reservoir. Other studies at wastewater stabilization reservoirs identified the percentage of fresh effluents as the main parameter effecting effluent quality which was sensitive to operational changes (Friedler et al., 2003), thus, serve for optimization of wastewater stabilization reservoirs (Mancini et al., 2017). The results describe the significance of both variables for design and operation of reservoirs. However, implementation of this approach for natural waterbodies, as lake, has not been performed yet.

Analysis of the results of the 90 model simulations provides insight into the likely difference in mean residence time under the three projected futures. They indicate that the predicted increase in mean residence time in 2060 under the BAU scenario is expected to be 19-50% higher than the mean residence time calculated for 2018 following a unprecedent five-year drought period. In addition, higher values of FYW10 in 2060 were obtained for the WATERSHED and LAKE scenarios compared to BAU scenario. This means that the water in the lake will be much younger for these scenarios. These results indicated that for the BAU scenario, with lower inflows and limited withdrawal, water exchange in the lake will decrease, as expected, and as a consequence the water residence time in the lake will increase, compared to the other scenarios. Since the forecast predicts long-term decreases in inflows and outflows throughout the period, at the end of the predicted period, i.e. 2060, is expected to be mostly affected by these changes, and therefore the highest mean residence time is expected then. Indeed, the decrease in inflows and outflows throughout the period is well noticed in the predicted range of the mean residence time of BAU scenario. However, importing water into the lake moderated the decline in inflows and therefore moderated the influence on the mean residence time, as predicted for the WATERSHED and LAKE scenarios. These results strengthen the need for importing water into the lake in order to mitigate the reduction of natural inflows due to climate change. In addition, the obtained differences between WATERSHED and LAKE scenarios indicate that according to our approach, even when the differences in annual flows are quite small between the two scenarios, the cumulative difference affect the water age distribution, and the mean residence time and as a result may impact the lake ecosystem.

Water time-scale variables play a crucial role in waterbodies since they affect many natural biochemical and physical processes (e.g., downstream output, water exchange, sedimentation, degradation and nutrients exchange). These significantly impact the water quality of the waterbodies. Hence, the study of lakes and reservoirs time-scale variables was found relevant for management purposes in numerous studies. For example, Camacho and Martin (2013) suggested to use time-scales to support the development of minimum flow criteria along with the design of watershed management programs of St. Louis Bay estuary, Mississippi. Qi et al. (2016) quantified a time-scale value, for lake ecosystem management of Lake Poyang (China), using data of a series of scenarios based on inflows, water levels, and meteorological data for different representative periods. These authors concluded that time-scale indicators may provide environmental monitoring and could help to determine an appropriate management policy. Gao et al. (2016) referred time-scale values as a practical water management measure, indicating persistence of pollution in the lake, and signaling the need for an appropriate response by authorities to address lake water quality problems. These examples emphasize the importance of accurate estimation of the water time-scales which is a substantial task not only from the management point of view but also for basic research.

Conclusions

We developed a novel approach for estimating the mean residence time of water within lake ecosystems under non-steady-state flow conditions. This approach based on the water age distribution, was adapted from the Leslie matrix model. Major advantages associated with estimating the water mean residence time in lakes were demonstrated by our unique approach of calculating the water age distribution. Opposed to the traditional residence time calculation, our water age distribution approach incorporates a “memory” component thereby providing a realistic solution, better representing the aquatic ecosystem in question.

The water age distribution approach produces the fraction and volume distribution curves of all age classes within the lake for each simulated timestep. Thus, in addition to mean residence time, the fraction of young water (FYW), which quantifies the “young” fraction of water in the lake is analyzed. Mean residence time relies on the complete hydraulic age distribution, hence it is less sensitive to changes in the annual inflows and outflows but succeedes in representing the underlying long-term trends. In contast, the FYWn relies on changes in inflows and outflows during the last n years and is much more sensitive to large interannual variation resulting in a larger variability than the mean residence time.

The developed approach was used for comparing three scenarios of importing water into Lake Kinneret, assuming climate change in the lake region. The approach indicated that even when the variation in flows is quite small the cumulative changes affect the water age distribution and thus the difference in mean residence time can be quite significant. This observation would not have been obvious based on the traditional-theoretical residence time calculation, which relies on current conditions and ignores the lake's conditions at previous timesteps.

The current approach offers a more accurate estimation of water mean residence time in lakes and can easily be adapted to other water bodies such as reservoirs, wetlands, and estuaries. It provides an overall picture of the water age structure, expressed by fraction and volume of age classes. Simulation results of future scenarios may serve as basis for determining lakes management policy, by controlling either inflows or outflows or both, that will affect both the mean residence time and the fraction of “young/old” water age classes.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.