Rising variance, changing skewness and critical slowing down: Early indicators of a tipping point in Andean Lake La Cocha

(1)

UNIVERSITY OF AMSTERDAM

BACHELOR THESIS

Rising variance, changing skewness and critical

slowing down: Early indicators of a tipping

point in Andean Lake La Cocha

Author: Supervisor:

dhr. K.C. (Koen) van Hattum dhr. dr. ir. J.H. (John) van Boxel

(2)

2

Abstract

Ecological and climatological systems are occasionally subject to an abrupt, non-linear change that is described as a tipping point. These changes are often characterized by great biodiversity and ecosystem services losses. Rising variance, changing skewness and critical slowing down (CSD) have been argued to be early indicators of an approaching tipping point. As recently freshwater lakes and arctic ponds dried up, the aim of this research is to verify that these early warning signals are also applicable to lake-level dynamics and in particular the modelled lake-level dynamics of Lake La Cocha. The research uses a rewritten version of an existing hydrological model to analyse the presence of these early indicators, which occur if the lake stops discharging. To detect the early indicators, the amount of precipitation, used as input for the model, is varied. The existing model yields a fit of 74% and a standard error of estimate 0.15 m with the observed lake levels. Subsequently, the rewritten model produces a 100% fit with the existing model. Both variance and CSD prove to be a valuable indicator as they both appear to rise preceding the tipping point. Skewness does rise prior to the tipping point, however, only with a fairly gentle constant slope. Early warning signals of an approaching tipping point can be more robust if both rising variance and CSD are identified. In addition, future research should focus on the identification of critical values instead of only looking at the changes of early indicators. Lastly, recent research found that flickering could also be a new early indicator. This method does not require a high-resolution time series, which are most often not available. Therefore, it could be a valuable contribution to the identification of tipping points in the future.

(3)

3

1. Introduction

An ecological or climatological system occasionally experiences a large, abrupt change in the structure and functioning of the system. These changes are generally associated with great biodiversity and ecosystem services losses (Scheffer et al., 2001). Holling (1973) was the first to acknowledge the existence of these changes and named it as ‘regime shifts’ to describe the behaviour of the non-linear dynamics of a system. The responsiveness of an ecosystem to disruptions, in the context of resistance and recovery, is what Holling (1973) called the ‘resilience’ of a system. Nowadays, the term ‘tipping point’ is more commonly used to describe a regime shift (IPCC, 2013).

The resilience of a system is the ultimate indication of the sensitivity to change to an alternative regime. A small change ensures that a forward shift takes place at F2 (fig. 1), also described as a catastrophic bifurcation. As a result, the system ends up in an alternative state. This change is an example of hysteresis, namely: when after the regime shift the driver of the system is restored to the initial value, it will not automatically turn back into the original regime (Scheffer et al., 2001).

Figure 1: The shift to an alternative state. Source: Scheffer et al. (2001) 1.1 Scientific and societal relevance

Identifying early indicators can be of great importance for the prediction of climate change, because they could make it possible to foresee a nearing tipping point (Scheffer et al., 2009). In addition, such warning signals can have an enormous influence on the way ecosystems are managed by identifying vulnerable systems and taking action to prevent regimes from shifting.

Smol (2012) recently discovered that ponds in the Canadian Artic have been drying out because of a higher temperature, which led to a higher evaporation. Lake Poopó in Bolivia has also dried up due to the combination of subtraction of water for irrigation and a low precipitation pattern (NASA, 2016; Revollo, 2001). Nevertheless, both water systems experienced a great decline of biodiversity and ecosystem services e.g. loss of an ancient aquatic ecosystem, breeding spots for birds, fishing industry. Therefore, it is important to study tipping points concerning lake-level and pond dynamics. If, in the case of lake La Cocha the lake drops below a certain value and stops discharging into the Guamués River, this will also reduce the flow of the Putumayao River. This will lead to a decrease in navigation opportunities for the local population. Furthermore, there will be less fertile land available at the lower Guamués River because the seasonal flooding of the nutrient-rich river sediment would cease to exist (WWF, 2001).

(5)

5 1.2 Current research and knowledge gap

Current research has concentrated on the study of ecological processes to detect early warnings signals that indicate an approaching regime shift. Examples of these ecological processes are: the eutrophication of lakes and oceans close to shore (Carpenter et al., 2011), changes in land cover in vast areas (Foley et al., 2003; Kéfi et al., 2012), degradation of coral reefs (deYoung et al., 2008) or past climate change (Lenton et al., 2008). Rising variance (Brock & Carpenter; Carpenter & Brock, 2006), changing skewness (Guttal & Jayaprakash, 2008) and the slowing down of the fluctuations, also known as the critical slowing down (CSD) (Dakos et al., 2008) have been argued to be early indicators of a tipping point.

Van Boxel et al. (2014) created a hydrological model that rather accurately simulates the lake level of Andean Lake La Cocha as a function of weather conditions. The model was used to reconstruct the precipitation change during the Holocene and Late Glacial through inverse hydrological modelling from a pollen-based lake level reconstruction. The lake currently discharges into consecutively the Guamués, Putumayo, San Miguel Rivers and finally the Amazon River. Nevertheless, van Boxel et al. (2014) dispute the probability of the lake discharging in the past. If the lake level drops below a certain value, causing the lake to stop discharging, this can be seen as passing a regime shift or a tipping point.

1.3 Research aims and questions

The aim of this research is to verify that early warning signals of tipping points, such as rising variance, changing skewness and CSD, are also applicable to lake-level dynamics and in particular the lake-level dynamics of Lake La Cocha. The research poses the question ‘Does rising variance, changing skewness and CSD indicate tipping points in the hydrological model of Andean Lake La Cocha?’. To answer the research question, four sub-questions were defined:

 To what extent does the rewritten model correlates with the model of van Boxel et al. (2014)?  Does the variance rise?

 Does the skewness change?  Do the fluctuations slow down?

This paper will provide a comprehensive insight in the feasibility of the detection of early indicators of tipping points and will critically discuss the results and practices involved.

(6)

6

Figure 2: Maps of Lake La Cocha. A: Geographical location of the lake. B: Land cover and elevation map. Source: van Boxel et al. (2014)

Figure 3: Main processes incorporated in the hydrological model of Lake La Cocha. Source: Van Boxel et al. (2014)

(7)

7

2. Scientific approach and methods

This section describes the materials, model and methods in this research. First, a description of the site will be provided. Then, the general functioning of the hydrological model is explained, followed by a description of the modern climate and lake-level data. Lastly, the methods to identify the tipping point and to analyse the early indicators will be discussed.

2.1 Site description

Lake La Cocha (fig. 2) is a freshwater lake in the south-west of Colombia, just north of the Ecuadorian border (106’N, 7709’W) (González-Carrenza et al., 2012). The lake is located on the eastern side of the Andes and lies at an elevation of 2780 m. The surface of the lake covers 42 km2 with 16 km in length, a maximum width of 6 km and depth of 75 m. The water volume is 1.55 x109 m3 _{(Matabanchoy & Lopez, 1999). The La Cocha catchment lies between 2780 m and 3650 m} elevation and covers an area of 214 km2_{(Moreno Díaz, 2004). The land cover of the basin consists of} namely wetlands, cloud forest, agricultural land and high montane grassland (páramo). (WWF, 2001; Gonzalez-Carrenza et al., 2012). The soils in the area can best be characterized as thick andosols with high porosity (Poulenard et al., 2001; Buytaert et al. 2005; Tonneijck et al. 2010).

2.2 The hydrological model

This research focuses on the Andean Lake La Cocha in Colombia. The hydrological model of Lake La Cocha (fig. 3) was created by Van Boxel et al. (2014) and was written in Delphi. This research implements the model in MATLAB. The model uses an integrated time step of 1 day and grid cells of 100 x 100 m2_{. The input variables of the model are daily values of precipitation, temperature,} sunshine hours, relative humidity and wind speed. The output variables of the model are discharge and evaporation. The lake level regulates the discharge of the river, whereas the evaporation is strongly determined by the temperature and thus radiation.

The precipitation is uniformly distributed over the catchment. Precipitation that falls directly into the lake is added to the lake volume and precipitation that falls on the land is added to the soil water content. Soil water beneath field capacity is stored in the soil whereas soil water above field capacity drains through the soil into the lake with a drainage time constant of 5 days. Overland flow appears if the soil water content exceeds the porosity. The water will run consequently into the river and the lake (van Boxel et al., 2014).

Evaporation from the lake is directly subtracted from the water volume. On land, the evapotranspiration is calculated by multiplying the open water evaporation with a crop factor. The evaporation is determined by a linear interpolation between potential evaporation at the field capacity and zero evaporation at the wilting point (van Boxel et al., 2014).

To determine the evaporation, the net radiation is required. The global radiation is calculated using the mean latitude of the catchment, Julian day number and cloudiness. Subsequently, the albedo of the surface types is used to calculate the reflected shortwave radiation. Stefan-Boltzmann’s law is applied to calculate the emitted long wave radiation and clouds are presumed to radiate as black bodies (van Boxel et al., 2014).

Van Boxel et al. (2014) assumed that the discharge (D) of the river fluctuates with a power of 2.5 of the lake level (H1) where the discharge coefficient (a) is 2.814 and the bottom of the river (b) is -0.050 m.

(8)

8 The Digital Terrain Model (DTM) obtained from van Boxel et al. (2014) holds 21370 grid cells of which 4187 cells currently fall within the lake. However, the size of the lake is subject to change due to lake-level variations. The surface type differs with elevation between water (lake), forest (cloud forest and agricultural land) and páramo (high montane grassland) with all three having different hydrological properties (table 1). The grid cells containing forest or páramo are divided into 50 m elevation classes. The calculations of radiation, evaporation and water content are preformed per elevation class as well as the determination of temperature, relatively humidity and atmospheric pressure, which also vary with elevation. A more detailed explanation of the hydrological model can be found in the supplementary material of van Boxel et al. (2014).

Table 1: Hydrological properties for different surface types. Source: Van Boxel et al. (2014)

Surface type Water Forest Páramo

Albedo 0.08 0.15 0.20 Soil depth [m] 1.5 1.5

Porosity 0.7 0.7

Water content at field capacity 0.5 0.5 Water content at wilting point 0.2 0.2 Crop factor 1.0 0.8 0.7 Drainage time constant [days] 5.0 5.0 Infiltration capacity [mm] 5.0 5.0

Table 2: Statistics of daily average climatological parameters and lake level for the period 1986-2001. Source: Van Boxel et al. (2014)

Parameter (units) Mean Standard deviation

Minimum Maximum Missing days (%) Precipitation (mm day-1₎ _3.7 _5.7* _0.0 _52.3 _2.4 Temperature (C) 11.6 1.0 6.2 15.0 17.8 Relative humidity (%) 86.7 3.8 62.0 100.0 20.3 Sunshine (h day-1₎ _2.5 _2.3 _0.0 _11.4 _3.4 Wind speed (m s-1₎ _1.3 _0.6 _0.1 _7.7 _4.9 Lake level (m) 1.38 0.30 0.77 2.79 0.0

*Distribution is strongly skewed

Figure 4: Mean monthly precipitation and temperature with minima and maxima. Source: Van Boxel et al. (2014)

(9)

9 2.3 Modern climate and lake-level data

The daily climate data (1986-2001) used as input for the hydrological model was retrieved from van Boxel et al. (2014). The original data is collected from the Encano meteorological station of the Instituto de Hidrología, Meteorología y Estudios Ambientales de Colombia (IDEAM) placed at 1 km from the lake shore at 2830 m elevation (109’N, 7711’W). The data consists of precipitation, temperature, relative humidity, sunshine and wind speed (table 2). The mean monthly precipitation (fig. 4) is 100 mm on average, with values of 150 mm in the wet period (April until July). The mean monthly temperature is around 12 C, but drops to almost 10.5 C in the coldest month (August). The performance of the model will be tested via cross correlation with the simulated lake-level values and the estimated lake-level values from van Boxel et al. (2014). Missing data in the observed climate and lake-level data was replaced by the monthly average of the month at issue.

2.4 Tipping point identification

The change from a state where the lake is discharging to an alternative state where the lake is not discharging is stated as the tipping point in this system. As there is a strong feedback between the precipitation and the discharge, this research uses a precipitation coefficient (P) to gradually reduce the input of water into the system (van Boxel et al., 2014). At a tipping point, the model will enter the forward loop in which the state of the lake will change from discharging to not discharging. Test runs with the model of van Boxel et al. (2014) indicate that the precipitation coefficients which affect the tipping point lie between P = 0.45 and 0.7. Therefore, the simulation will be run between these precipitation coefficients with step size P = 0.01. The tipping point will be indicated as the point of entry of the forward loop.

2.5 Early indicators analysis

The lake-level series that will be produced are first detrented to effectively apply and examine the statistical analysis of the early indicators. For each precipitation coefficient the lake-level series will run until the start level and end level differ less than 1 mm, which can be considered as the equilibrium.

After the equilibrium is reached, the variance and the skewness of the lake-level time series will be taken to detect the rising variance and changing skewness. The CSD will be tested via an autocorrelation spectrum with a varying lag of 1 year. The CSD will be reviewed at the point where the R2_{of the autocorrelation falls to 0.5, which indicates that 50% of the time series is described by its} lagged series.

(10)

10

3. Results

This section describes the results of the research. First, the correlation between the lake levels of the rewritten model in MATLAB and the model from van Boxel et al. (2014) will be presented. Then, the identification of the tipping point will be reported, followed by the outcome of the detection of the early indicators.

3.1 Model evaluation

The rewritten model in MATLAB reports a correlation coefficient (R) of 1.00 with the lake-level values of van Boxel et al. (2014). In addition, all the other modelled variables such as radiation, evaporation, runoff, drainage and discharge have a R of 1.00. However, this research has made an alteration to the model to improve the correlation between the modelled and observed lake levels. It was found that the reference height used to estimate the temperature, relative humidity and atmospheric pressure in the original model was set equal to the lake level (2780 m) instead of the reported height of the meteorological station (2830 m). After changing the reference height, the model yielded a coefficient of determination (R2_{) of 73.9% with a standard error of estimate (SEE) of 0.154} m whereas the modelled lake levels of van Boxel et al. (2014) report a R2_{of 73.5% with a SEE of} 0.155 m. The observed and modelled lake-level values are visualized in fig. 9. The fit of the model improved with 0.40% and the SEE dropped with 1 mm. Nevertheless, the alteration did not substantially increase the performance of the model.

3.2 Tipping point identification

The change from a state where the lake is discharging to an alternative state where the lake is not discharging is stated as the tipping point in this system. The system shows that this change of states occurs between P = 0.60 and 0.56 (fig. 5) as the discharge shifts from being active 100% to 0% of the days. The decreasing input of precipitation causes the lake level to drop after passing the tipping point, starting at P = 0.57. However, between P = 0.50 and 0.48 the lake level seems to stabilize at a lake-level depth of -30 m. This can be explained by the bathometry of the lake.

After reaching the equilibrium of the lake level per precipitation coefficient, the time series of the lake-level data were analysed by calculating the variance, skewness and an autocorrelation spectrum.

As P moves from 0.70 to 0.45, the variance decreases first (fig. 6), but starts rising at P = 0.62 and continues until the end of the regime shift after which it starts to drop. Subsequently at P = 0.51 the variance tends to rise again, followed by a drop at P = 0.50 until P = 0.49 and a small rise between

P = 0.49 and P = 0.48 again.

The skewness of the lake level changes in a moderate manner while approaching the tipping point (fig. 7). As the P moves towards 0.60 the frequency distribution becomes right-skewed. Furthermore, a small peak in the skewness occurs between P = 0.50 and P = 0.48. Both the variance and the skewness show the small period of lake-level stabilization. Therefore, it can also be described as a change of states.

The autocorrelation is assessed at the point where R2_{= 0.500 (indicated with the dotted line in} fig. 8). As P moves from 0.64 to 0.60, the lag of the autocorrelation increases with 1 month (from 2.5 to 3.5 months). As the regime shift occurs, the lag continues to increase with 2.5 months up to 6.0 months, after which it stabilizes around 5.5 months.

(11)

11

Figure 5 (above): Discharge (% days) and mean lake level as function of precipitation coefficient (P). Figure 6 (below): Variance as function of Precipitation Coefficient (P).

Figure 7: Frequency distributions of the lake level at different Precipitation Coefficients (P) and skewness as a function of Precipitation Coefficients (P).

(12)

12

Figure 8: Autocorrelation spectrum with a varying lag up to 1 year for different Precipitation Coefficients

(13)

13

4. Discussion

This research examines the possibilities to detect early indicators of a tipping point in the lake-level dynamics of Andean Lake La Cocha using a hydrological model. First, the model, its assumptions and the response time of the lake level to changing environmental conditions will be discussed. Then, the methods to identify the tipping point will be considered, followed by a discussion about the early indicators analysis. Finally, the applicability of early indicators will be reviewed.

4.1 Model evaluation

The model provides a good fit with the observed lake levels and therefore offers a valuable representation of the reality. Van Boxel et al. (2014) optimized only two model parameters (a and b in Eq. 1) and reached an R2_{of 0.735 with a standardized error of estimate 0.155 m. The alteration of the} reference height used to estimate the temperature, relatively humidity and atmospheric pressure improved the R2_{of the model with only 0.40% and the SEE with 1 mm, which did not substantially} upgrade the performance of the model. The difference between de modelled and the observed lake levels can partly be explained by the use of the precipitation input. The precipitation data from only one meteorological station were used to determine the precipitation for the whole catchment. If a strong rain storm passed over the meteorological station, the amount of precipitation was generalized for the whole catchment, resulting in an overestimation of the lake level. This can be noticed in fig. 9 where some peaks in the modelled lake levels cannot be seen in the observed lake levels.

Parameters involving the water retention capacity of the soil and the pace at which soil water drains into the lake were estimated instead of calibrated. Van Boxel et al. (2014) argue that it would not strongly influence the model because the calibration of the discharge parameters a and b in Eq. 1 appeared to be an efficient mechanism to increase the performance of the model.

After the system has passed the tipping point, the equilibrium value is determined by the feedback between the lake level and evaporation. In this situation, the return time to the equilibrium is much longer because there is no discharge feedback. In addition, although the model recounts the surface types for every lake level, the surface area of the lake will not change severely because the slopes around the lake are rather steep. Therefore, the difference between evaporation from the vegetation or evaporation directly from the lake surface will not vary strongly, except if there are arid or semi-arid conditions (van Boxel et al., 2014).

There is a big difference in the number of cells per elevation between -29 m and -30 m (table 1A in Appendix A), which indicates that from one depth to another, a rather large area changes from surface type. As the crop factor for grid cells outside the lake is lower than for grid cells inside the lake (table 1), less water evaporates after the depth changes from -29 to -30, resulting in a short period of lake-level stabilization. The difference between the number of grid cells per surface type between a depth of -29 m and -30 m, which caused the lake level to stabilize for a short period, can be subscribed to estimation of the bathometry by van Boxel et al. (2014). The grid cells that represent the south corner of the lake (fig. 2b) almost all contain an elevation of 2750 m, which concur with the transition from depth -29 m to -30 m. Van Boxel et al. (2014) developed the model to estimate precipitation changes from the paleo lake levels. Since the length of the core was only 12 m, lake levels beyond a depth of -12 m could not be produced. Consequently less attention was paid to the bathymetry below -12 m.

(14)

14 4.2 Tipping point identification

To examine the behaviour of the lake level around the tipping point, this research only varies one climatological variable, namely precipitation. In reality, changing systems endure a greater collection of changing variables which can have different influences. To assess the influence of other climatological variables, model simulations were run with varying temperatures. All the other variables were kept constant. The tipping point appears at a temperature coefficient (T) of 4.10, which implies a temperature rise around 35 °C, as the new mean temperature at the meteorological station (2830 m) becomes 47 °C. The lake is discharge free at T = 4.80, which implies a mean temperature of 56 °C. This fairly high temperature rise is not likely to happen as even the most dramatic predictions about global warming does not come close to this value (IPCC, 2013). Therefore, the influence of temperature on the lake-level dynamics of Lake La Cocha can be considered as non-substantial. However, for the grid cells that are located outside the lake, the hydrological model first performs the calculations for the infiltration, drainage and run off, followed by the evaporation. Thus, all soil water above field capacity first drains or flows into the lake, after which it can be subject to evaporation. As a rise in temperature directly leads to an increase in evaporation, the hydrological model underestimates the actual evaporation of the soil grid cells, eventually leading to a higher lake level. Therefore, while only varying the temperature, the tipping point can only be approached if an extreme temperature rise is intended.

In addition, not only climatological and ecological factors need to be reviewed. Many ecological disasters involving lake-level dynamics have been majorly influenced by anthropogenic factors, such as the over extraction of water for agriculture and industry in lake systems. This has already happened to Lake Poopó, Bolivia (Revollo, 2001), the Aral Sea, Kazakhstan (Edelstein, 2012), Lake Urmia, Iran (AghaKouchak et al., 2015), the Dead Sea (Yechieli et al., 2006) and Lake Chad (Coe & Foley, 2001). In the case of Lake La Cocha, water extraction for agriculture is a minor factor as the catchment endures rather wet conditions all year round. Nevertheless, future models studying lake-level dynamics should therefore, if applicable, also incorporate anthropogenic factors to identify nearing tipping points.

The variance appears to be a good indicator of a tipping point in this system as it clearly starts to rise prior to the tipping point. These findings are in line with recent research (Carpenter & Brock, 2006; Carpenter et al., 2011; Lenton, 2011) Nonetheless, the period in which it changes is rather short. Therefore, it could be difficult to adequately deal with an approaching tipping point.

The skewness also tends to rise preceding the tipping point. However, the slope is fairly gentle and does not change while approaching the tipping point. Only after passing the tipping point, the slope changes. Thus, the skewness does not appear be a good indicator of an approaching tipping point in this system. This stands in contrast with other research where the skewness occurred to be a valuable warning signal (Guttal & Jayaprakash, 2008; Carpenter & Brock, 2006).

The increase in autocorrelation is properly noticeable and evens seems to increase while approaching the tipping point. The CSD is therefore a good instrument to use as an early indicator of a tipping point in this system, which is comparable with other research (Dakos et al., 2008; Lenton 2011; Kéfi et al., 2012).

(15)

15 4.4 Applicability of early indicators

Ditlevsen & Johnsen (2010) argue that an early warning signal of an approaching tipping point can only be obtained if both rising variance and CSD are identified in a climatological or any other dynamical system, as it has been proven difficult to identify statistical significant changes in the variance and autocorrelation. Moreover, Significant increases of an indicator seem only to appear if a tipping point has been passed and is already changing to another state, often too late to stop the regime from shifting and to adequately respond to the sudden changes. To improve the ability to avoid regimes to shift Biggs et al. (2009) argue that instead of detecting changes in the early indicators, research should focus on the definition of critical values. Nevertheless, the application of these statistical instruments can be difficult to examine in real life time series. There have been reports of a class of ecological models where tipping points were passed without any warning signals, which would also be likely to happen in nature. Many ecological systems do not produce continuous potentials and thus will not show the described early indicators prior to a tipping point (Hastings & Wysham, 2010).

In addition to the early indicators in this research, Wang et al. (2012) found that flickering, which induces a system to keep switching between different states in reaction to rather large influences, appears to be a possible new early warning signal. It seems to be more robust and easier to detect, as it does not require a high resolution time series. As high resolution time series are most often not available, this new method could be a valuable contribution to the identification of tipping points in the future.

(16)

16

5. Conclusion

This research posed the question: ‘Does rising variance, changing skewness and CSD indicate tipping points in the hydrological model of Andean Lake La Cocha?’. To asses this question, this research used a rewritten MATLAB model of an existing hydrological model of Andean Lake La Cocha from van Boxel et al. (2014). The rewritten model in MATLAB provided a good fit with the modelled lake levels from van Boxel et al. (2014). After implementing a small alteration, the model provided an improved fit, but this did not substantially change the results.

The tipping point of the system was approached by gradually decreasing the precipitation using a precipitation coefficient. The tipping point was identified as the value of the precipitation coefficient at which the lake changed from a discharging to a non-discharging state. The tipping point was located at a precipitation coefficient of 0.60 at which the regime started to shift. The system was found to be discharge free at a precipitation coefficient of 0.56.

To detect the early indicators, the detrended time series of the lake level around the tipping point were analysed by calculating the variance, skewness and CSD. The variance and CSD tended to be a valuable indicator of a nearing tipping point as they both appeared to rise preceding the tipping point. The skewness also rose prior to the tipping point, however, with a fairly gentle constant slope, which did not change. The skewness has therefore not been proven to be a valuable indicator of an approaching tipping point.

The applicability of early indicators in ecological and climatological models was discussed. Early warning signals of an approaching tipping point can be more robust if both rising variance and CSD are identified. In addition, future research should focus on the identification of critical values instead of only looking at the changes of early indicators. Lastly, recent research detected a new early indicator that does not require a high resolution time series, which are most often not available. This new finding could be a great benefit to field of research involving tipping points.

(17)

17

6. Acknowledgements

Sincere gratitude goes out to dr. ir. J.H. (John) van Boxel for supervising this research and providing all the materials of the hydrological model and its explanation.

(18)

18

7. References

AghaKouchak, A., Norouzi, H., Madani, K., Mirchi, A., Azarderakhsh, M., Nazemi, A., Nasrollahi, N., Farahmand, A., Mehran, A. and Hasanzadeh, E. (2015). Aral Sea syndrome desiccates Lake Urmia: call for action. Journal of Great Lakes Research, 41(1), 307-311.

Biggs, R., Carpenter, S. R., & Brock, W. A. (2009). Turning back from the brink: detecting an impending regime shift in time to avert it. Proceedings of the National academy of

Sciences, 106(3), 826-831.

Brock, W. A., & Carpenter, S. R. (2006). Variance as a leading indicator of regime shift in ecosystem services. Ecology and Society, 11(2).

Buytaert, W., Wyseure, G., De Bièvre, B., Deckers, J. (2005) The effect of land-use changes on the hydrological behaviour of Histic Andosols in south Ecuador. Hydrol Process 19:3985–3997 Carpenter, S.R., & Brock, W.A. (2006). Rising variance: a leading indicator of ecological

transition. Ecology Letters, 9(3), 311-318. http://dx.doi.org/10.1111/j.1461-0248.2005.00877.x Carpenter, S. R., Cole, J. J., Pace, M. L., Batt, R., Brock, W. A., Cline, T., T., Coloso, J., Hodgson,

J.R., Kitchell, J.F., Seekell, D.A. & Smith, L. (2011). Early warnings of regime shifts: a whole-ecosystem experiment. Science, 332(6033), 1079-1082.

Coe, M. T., & Foley, J. A. (2001). Human and natural impacts on the water resources of the Lake Chad basin. Journal of Geophysical Research: Atmospheres, 106(D4), 3349-3356.

deYoung, B., Barange, M., Beaugrand, G., Harris, R., Perry, R. I., Scheffer, M., & Werner, F. (2008). Regime shifts in marine ecosystems: detection, prediction and management. Trends in

Ecology & Evolution, 23(7), 402-409.

Dakos, V., Scheffer, M., van Nes, E.H., Brovkin, V., Petoukhov, V., & Held, H. (2008). Slowing down as an early warning signal for abrupt climate change. Proceedings Of The National

Academy Of Sciences, 105(38), 14308-14312. http://dx.doi.org/10.1073/pnas.0802430105

Ditlevsen, P. D., & Johnsen, S. J. (2010). Tipping points: early warning and wishful thinking. Geophysical Research Letters, 37(19).

Edelstein, M. R. (2012). Aral Sea Demise as a Dry Run for Climate Change: From Cumulative to Cascading Impacts. In Disaster by Design: The Aral Sea and its Lessons for

Sustainability (pp. 415-441). Emerald Group Publishing Limited.

Foley, J. A., Costa, M. H., Delire, C., Ramankutty, N., & Snyder, P. (2003). Green surprise? How terrestrial ecosystems could affect earth’s climate. Frontiers in Ecology and the

Environment, 1(1), 38-44.

González-Carranza Z., Hooghiemstra H., Vélez M.I. (2012) Major altitudinal shifts in Andean vegetation on the Amazonian flank show temporary loss of biota in the Holocene. Holocene 22:1227–1241

Guttal, V. and Jayaprakash, C. (2008), Changing skewness: an early warning signal of regime shifts in ecosystems. Ecology Letters, 11: 450-460. doi:10.1111/j.1461-0248.2008.01160.x

(19)

19 Hastings, A., & Wysham, D. B. (2010). Regime shifts in ecological systems can occur with no

warning. Ecology letters, 13(4), 464-472.

Holling, C.S. (1973). Resilience and Stability of Ecological Systems. Annual Review of Ecology and

Systematics. 4 (2), pp.1-23

IPCC (2013): Summary for Policymakers. In: Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change [Stocker, T.F., D. Qin, G.K. Plattner, M. Tignor, S.K. Allen, J. Boschung, A. Nauels, Y. Xia, V. Bex and P.M. Midgley (eds.)]. Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA. http://www.ipcc.ch/pdf/assessment-report/ar5/wg1/WG1AR5_SPM_FINAL.pdf

Kéfi, S., Dakos, V., Scheffer, M., Van Nes, E. H., & Rietkerk, M. (2013). Early warning signals also precede non‐catastrophic transitions. Oikos, 122(5), 641-648.

Lenton, T. M., Held, H., Kriegler, E., Hall, J. W., Lucht, W., Rahmstorf, S., & Schellnhuber, H. J. (2008). Tipping elements in the Earth's climate system. Proceedings of the national Academy

of Sciences, 105(6), 1786-1793.

Lenton, T. M. (2011). Early warning of climate tipping points. Nature Climate Change, 1(4), 201. Matabanchoy, J.C., López, O.D. (1999) Pasto entre la sed y el agua. Associacio´n para el Desarollo

Campesino, Medellín, Colombia

Moreno Díaz, C.A. (2004) Proyecto de incentivos para la laguna de La Cocha como sitio Ramsar; Informe final primera fase. Instituto Alexander Von Humboldt, Bogotá, Colombia, p 89 NASA. (2016). Bolivia’s Lake Poopó Disappears : Natural Hazards. Retrieved from

https://earthobservatory.nasa.gov/NaturalHazards/view.php?id=87363

Poulenard, J., Podwojewski, P., Janeau, J.L., Collinet, J. (2001) Runoff and soil erosion under rainfall simulation of Andisols from the Ecuadorian páramo: effect of tillage and burning. Catena 45:185–207

Revollo, M.M. (2001). Management issues in the Lake Titicaca and Lake Poopó system: Importance of developing a water budget. Lakes And Reservoirs: Research And Management, 6(3), 225-229. http://dx.doi.org/10.1046/j.1440-1770.2001.00151.x

Scheffer, M., Carpenter, S.R, Foley, J.A., Folke, C., & Walker, B.H. (2001). Catastrophic shifts in ecosystems. Nature, 413(6856), 591-596. http://dx.doi.org/10.1038/35098000

Scheffer, M., Bascompte, J., Brock, W.A., Brovkin, V., Carpenter, S.R., Dakos, V., Held, H., Nes, E.H. van, Rietkerk, M., Sugihara, G. (2009) Early-warning signals for critical transition.

Nature, 461, 53-59.

Smol, J. (2012). Climate Change: A planet in flux. Nature, 483(7387), S12-S15. http://dx.doi.org/10.1038/483s12a

(20)

20 Tonneijck F.H., Jansen B., Nierop K.G.J., Verstraten J.M., Sevink J., De Lange L. (2010) Towards

understanding of carbon stocks and stabilization in volcanic ash soils in natural Andean ecosystems of northern Ecuador. Eur J Soil Sci 61:392–405

Van Boxel, J.H., González-Carranza, Z., Hooghiemstra, H., Bierkens, M., & Vélez, M.I. (2014). Reconstructing past precipitation from lake levels and inverse modelling for Andean Lake La Cocha. Journal Of Paleolimnology, 51(1), 63-77. http://dx.doi.org/10.1007/s10933-013-97551 Wang, R., Dearing, J. A., Langdon, P. G., Zhang, E., Yang, X., Dakos, V., & Scheffer, M. (2012).

Flickering gives early warning signals of a critical transition to a eutrophic lake state. Nature, 492(7429), 419.

WWF. (2001). Managing Rivers Wisely: La Cocha case study. Cali: WWF. Retrieved from http://d2ouvy59p0dg6k.cloudfront.net/downloads/mrwlacochacasestudy.pdstudy.pdf

Yechieli, Y., Abelson, M., Bein, A., Crouvi, O., & Shtivelman, V. (2006). Sinkhole “swarms” along the Dead Sea coast: reflection of disturbance of lake and adjacent groundwater

(21)

21

8. Appendices

8.1 Appendix A – Extra figures and tables

Table A1: Differences in number of cells per elevation.

Depth (Lake

Level) Elevation Count Difference -40 2740 1866 0 -39 2741 1922 56 -38 2742 1976 54 -37 2743 2030 54 -36 2744 2089 59 -35 2745 2144 55 -34 2746 2202 58 -33 2747 2257 55 -32 2748 2308 51 -31 2749 2373 65 -30 2750 2429 56 -29 2751 2702 273 -28 2752 2742 40 -27 2753 2788 46 -26 2754 2834 46 -25 2755 2883 49 -24 2756 2920 37 -23 2757 2962 42 -22 2758 3005 43 -21 2759 3044 39 -20 2760 3099 55

Figure A1: Lake-level time series for different value of precipitation coefficient (P) during the regime shift

(22)

22 8.2 Appendix B – Digital material

A folder containing all digital material is available on request by sending an email to kc.vanhattum@gmail.com. The folder consists of:

 Digital Terrain Model (DTM)

 Daily values of precipitation, temperature, sunshine hours, relative humidity and wind speed from the climatological station.

 Hydrological model in MATLAB script.

 Output file of the lake-level and discharge time series after simulation used to perform the statistical analysis and to detect the regime shift.

(23)

23 8.3 Appendix C – Script

% Hydrological Model Lake La Cocha, Colombia % University of Amsterdam

% Bachelorproject Earth Science 2018 % Student: Koen van Hattum

% Supervisor: John van Boxel

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % I N I T I A L I S A T I O N % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %

clear close all

% Loading observed data from meteo station

data = load('DailySeries_K1.txt'); Date = data(:,1); ObsLakeLevel = data(:,2); Precip = data(:,3); ObsWind = data(:,4); ObsTemp = data(:,5); ObsSun = data(:,6); ObsHumid = data(:,7); ObsEvap = data(:,8);

% Loading standaard output daily series van Boxel et al. (2014)

data2 = load('StandaardOutput-DailySeries.txt'); Date2 = data2(:,1); Temp2 = data2(:,2); NetR2 = data2(:,3); Precip2 = data2(:,4); Evap2 = data2(:,5); Runof2 = data2(:,6); Drain2 = data2(:,7); Discharge2 = data2(:,8); Level2 = data2(:,9);

% transform Date into year, month, day

ConvDate = num2str(Date);

Year = str2num(ConvDate(:,1:4)); Month = str2num(ConvDate(:,5:6)); Day = str2num(ConvDate(:,7:8));

% Transform to Julian day numbers

Dates = datetime(Year,Month,Day); J = day(Dates,'dayofyear');

% Load DTM and intilialize Matrix

[Longitude, Latitude, DEM, OPT] = surfer_read('Catchment100x100-New.grd'); [ny, nx] = size(DEM); % Grid size dx = 100; % Pixel length in X-direction [m] dy = 100; % Pixel length in Y-direction [m]

(24)

24 % Time constants

StartTime = 1; % Start Time for simulation [day]

EndTime = 5844; % Time at which simulation Ends [day]

dt = 1; % Calculation time step

Time = StartTime; % Start simulation at first day

% Hydrological properties [Outside, Lake, Forest, Paramo]

Albedos = [0,0.08,0.15,0.20]; % Fraction of radation reflected for water, forest and paramo

SoilDepths = [0.00, 9.99, 1.50, 1.50]; % Soil depth for vegetation below and paramo [m]

Porosities = [0,1,0.7,0.7]; % Porosity of the soil [m3/m3]

FieldCaps = [0,1,0.5,0.5]; % Water content at field capacity [m3/m3]

Wiltings = [0,0,0.2,0.2]; % Water content at wilting point [m3/m3]

CropFactors = [0,1.0,0.8,0.7]; % Crop Factor water, forest, paramo

Infiltrations = [0,99,5.0,3.0]; % Infiltration Capacity [mm]

RunofFractions = [0,0,0.4,0.6]; % Runoff fractions

DrainConstants = [1,1,5.0,5.0]; % Drainage time constant [days]

% System Constants

meanlat = degtorad(mean(Latitude)); % Lattitude in Rad

CellKm2 = dx * dy / 10^6; % Cell area [km2]

LATENT = 2.5*10^6; % latent heat of vaporization for water [J/kg]

CP = 1006; % Heat capacity of dry air at constant pressure [J/kg/K]

EPSILON = 18/29; % molecair weight of water [kg/kmol]

P = 1013.25; % barometric pressure [mbar]

Omega = 2*pi / 24; % units hr^-1

Boltzmann = 5.67 * 10 ^ -8; % stefan Boltzmann constant

AngV = 7.27*10^-5; % Angular velocity of the earth [rad/s]

I0 = 1367; % Solar constant [W/m2]

Tlinke = 2.7; % Linke's turbidity factor

LapseRate = 0.6/100; % 0.6 K per 100 m

Temp_forestLine = 9.5; % Temperature till which trees grow [Celsius]

TRef = 11.62; % First temperature at height 2780 [Celsius]

ZRef = 2830; % Initial lake level and reference point temp. [m]

LakeElev = 2780; % Elevation of lake [m]

LakeLevel = 0.972; % Initialize lake level [m]

Bottom = - 0.050; % Bottom of the river [m]

DischargeCoeff = 2.814; % Discharge Coefficient [m^0.5/s]

(25)

25

DrainConstant = 5.0; % Drainage time constant [days]

%%% I N I T I A L I Z E A R R A Y S A N D M A T R I C E S %%%

% Create matrix for landcover en heightclasses

LandCover = zeros(ny,nx); Number = zeros(ny,nx); % Create arrays PrecipArray = ones(30,1); Runoff = zeros(1,30); Drainage = zeros(1,30); ExtraRunoff = zeros(1,30); Infiltration = zeros(1,30); Sunshine = zeros(1,30); SkyCover = zeros(1,30); PotEvap = zeros(1,30); Evap = zeros(1,30); WaterContent = zeros(1,30); Elevation = zeros(1,30); Temperature = zeros(1,30); SurfType = zeros(1,30); Albedo = zeros(1,30); SoilDepth = zeros(1,30); Porosity = zeros(1,30); FieldCap = zeros(1,30); Wilting = zeros(1,30); CropFactor = zeros(1,30); InfilCap = zeros(1,30); RunofFract = zeros(1,30); DrainConst = zeros(1,30); PotRad = zeros(1,30); TransR = zeros(1,30); Direct = zeros(1,30); CT = zeros(1,30); Diffuse = zeros(1,30); Global0 = zeros(1,30); Global1 = zeros(1,30); Global = zeros(1,30); Reflect = zeros(1,30); Kelvin = zeros(1,30); Emission = zeros(1,30); Esaturation = zeros(1,30); Eact = zeros(1,30); EpsAir0 = zeros(1,30); EpsAir = zeros(1,30); Atmosph = zeros(1,30); NetRad = zeros(1,30); Psych = zeros(1,30); S = zeros(1,30); M = zeros(1,30); Count = zeros(1,30);

% Create first water content array at field capacity

for i = 1:30

if i == 1

WaterContent(i) = 1; WaterContent(i+1) = 1; else ,WaterContent(i) = 0.5;

(26)

26 end end %%% I N I T I A L I Z E S T O R I N G OF V A R I A B L E S %%% StorePrecip = zeros(EndTime,1); StoreEvap = zeros(EndTime,1); StoreNetRad = zeros(EndTime,1); StoreRunoff = zeros(EndTime,1); StoreDrainage = zeros(EndTime,1); StoreDischarge = zeros(EndTime,1); StoreLevel = zeros(EndTime,1); % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % D Y N A M I C P A R T % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %

while Time <= EndTime

%%% D E F I N E T H E D E M %%%

% Locate the different surface types

ZStep = 50;

ZStart = ZRef - 2*ZStep;

ZForestLine = ZRef + (TRef - Temp_forestLine) / LapseRate;

for i = 1:ny for j = 1: nx

if DEM(i,j) > 9000 LandCover(i,j) = 1; elseif DEM(i,j) <= LakeElev LandCover(i,j) = 2; elseif DEM(i,j) > ZForestLine LandCover(i,j) = 4; else ,LandCover(i,j) = 3; end if LandCover(i,j) == 1 Number(i,j) = 1; elseif LandCover(i,j) == 2 Number(i,j) = 2;

else ,Number(i,j) = 3 + round((DEM(i,j) - ZStart) / ZStep); end

end end

for i = 1:2

Number(1:ny,1:nx) == i;

Count(i) = sum(sum( Number == i)); Elevation(i) =2780; Temperature(i) = ObsTemp(Time); SurfType(i) = i; Albedo(i) = Albedos(SurfType(i)); SoilDepth(i) = SoilDepths(SurfType(i)); Porosity(i) = Porosities(SurfType(i)); FieldCap(i) = FieldCaps(SurfType(i)); Wilting(i) = Wiltings(SurfType(i)); CropFactor(i) = CropFactors(SurfType(i));

(27)

27 InfilCap(i) = Infiltrations(SurfType(i)); RunofFract(i) = RunofFractions(SurfType(i)); DrainConst(i) = DrainConstants(SurfType(i)); end Z = ZStart; for i = 3:30 Elevation(i) = Z; Z = Z + ZStep;

Temperature(i) = TRef - (Elevation(i)-ZRef) * LapseRate; if Temperature(i) >= Temp_forestLine

SurfType(i) = 3;

Count(i) = sum(sum( Number == i));

Temperature(i) = ObsTemp(Time)- (Elevation(i)-ZRef) * ...

LapseRate; Albedo(i) = Albedos(SurfType(i)); SoilDepth(i) = SoilDepths(SurfType(i)); Porosity(i) = Porosities(SurfType(i)); FieldCap(i) = FieldCaps(SurfType(i)); Wilting(i) = Wiltings(SurfType(i)); CropFactor(i) = CropFactors(SurfType(i)); InfilCap(i) = Infiltrations(SurfType(i)); RunofFract(i) = RunofFractions(SurfType(i)); DrainConst(i) = DrainConstants(SurfType(i)); else ,SurfType(i) = 4;

Count(i) = sum(sum( Number == i));

Temperature(i) = ObsTemp(Time)- (Elevation(i)-ZRef) * ...

LapseRate; Albedo(i) = Albedos(SurfType(i)); SoilDepth(i) = SoilDepths(SurfType(i)); Porosity(i) = Porosities(SurfType(i)); FieldCap(i) = FieldCaps(SurfType(i)); Wilting(i) = Wiltings(SurfType(i)); CropFactor(i) = CropFactors(SurfType(i)); InfilCap(i) = Infiltrations(SurfType(i)); RunofFract(i) = RunofFractions(SurfType(i)); DrainConst(i) = DrainConstants(SurfType(i)); end end

% Locate forest & paramo (Soil) and Lake

Soil = find(SurfType == 3 | SurfType == 4); Lake = find(SurfType == 2);

%%% T E M P E R A T U R E - P R E S S U R E - S K Y C O V E R %%%

% set all values to zero

SkyCover = zeros(1,30); PotRad= zeros(1,30); Global= zeros(1,30); NetRad= zeros(1,30); PotEvap = zeros(1,30); Evap = zeros(1,30); % define pressure

Pressure = CalcPressure (Elevation);

(28)

28

% define actual evaporation

RHref = ObsHumid(Time);

Eact = RHref /100 * CalcEsat_mbar (ObsTemp(Time));

% define relative humidity

Esat = 6.1078*exp(17.2694*Temperature ./...

(237.3+Temperature));

RelHum = 100 * Eact ./ Esat;

% for diagnostics

RelHum = min(RelHum, 100);

% Calc. Sky cover using relative humidity

DayLength = CalcDayLength(meanlat,J(Time)); if RelHum >= 100 SunShine(1:30) = 0; SkyCover(1:30) = 1; elseif ObsHumid(Time) >= 100 SunShine(1:30) = 0; SkyCover(1:30) = 1;

else ,SunShine(1:30) = ObsSun(Time) .* (100-RelHum(1:30))./...

(100-ObsHumid(Time));

SkyCover(1:30) = 1 - SunShine(1:30) / DayLength; end %%% C R E A T E P R O P E R T Y T A B L E %%% Properties = table(Count',SurfType',Elevation',Temperature',... RelHum',Pressure',SkyCover',Albedo',SoilDepth',Porosity',... FieldCap',Wilting',CropFactor',InfilCap',RunofFract',DrainConst'); % % % R A D I A T I O N % % %

% Calc. orbital correction paramater (F) [-]

OrbitCor = 1 + 0.033 * cos(2*pi*J(Time)/366);

DeclRad = degtorad(-23.45)*cos(2*pi*(J(Time)+10)/366);

HalfDay = DayLength/2;

PotRadDay = 2/Omega * I0 * OrbitCor * (Omega* HalfDay *...

sin(meanlat) * sin(DeclRad) + cos(meanlat) *...

cos(DeclRad) * sin(Omega*HalfDay)) / 24;

% W/m2

% Calc. Cosine of the zenith angle (cosZ) [-]

CosZenit3 = sin(meanlat) * sin(DeclRad) + ...

cos(meanlat) * cos(DeclRad) * cos(Omega*3);

% Zenit angle 3:00 WZT

for i = 1:length(Count)

PotRad(i) = PotRadDay; % Potential radiation

M(i) = (Pressure(i) ./ 1013.25) * (1 / CosZenit3 );

% M = optical mass

TransR(i) = 0.959-0.0636*M(i)+0.00282*(M(i)^2)-6*(10^-5)*...

(M(i)^3)+ 4*(10^-7)*(M(i)^4);

(29)

29

Direct(i) = PotRad(i) * power(TransR(i),Tlinke);

% direct radation

CT(i) = power(Direct(i)/PotRad(i),0.363);

% (CT) van Boxel 2002 (16)

Diffuse(i) = Direct(i)*(-0.035 + 2.236 * CT(i) - 4.000 *... % diffuse radiation S.15

(CT(i)^2) + 2.134 * (CT(i)^3) - 0.272 * (CT(i)^4)); Global0(i) = Direct(i) + Diffuse(i);

% Global rad under a clear sky

Global1(i) = 0.40 * Direct(i);

% Global rad under a clouded sky

Global(i) = SkyCover(i) .* Global1(i) + (1-SkyCover(i))...

* Global0(i);

Reflect(i) = Albedo(i) .* Global(i);

% Reflected shortwave

Kelvin(i) = Temperature(i) + 273.15;

Emission(i) = 0.98 * Boltzmann .* (Kelvin(i).^4);

% Emission long wave

Esaturation(i) = CalcEsat_mbar(Temperature(i));

% Saturation vapour pr.

Eact(i) = Esaturation(i) .* RelHum(i) ./ 100;

EpsAir0(i) = 0.520 + 0.065 * sqrt(Eact(i));

% Ref: Brunt 1932

EpsAir(i) = SkyCover(i) .* 1.0 + (1-SkyCover(i)) .* EpsAir0(i); Atmosph(i) = EpsAir(i) .* Boltzmann .* (Kelvin(i).^4);

% Atmospheric LW rad

NetRad(i) = Global(i) - Reflect(i) + Atmosph(i) - Emission(i);

% Net Radiation

end

%%% I N F I L T R A T I O N A N D R U N O F F %%%

% Define precipitation for current day

ObsPrecip = PrecipArray .* Precip(Time);

% Calculate infiltration and runoff

if ObsPrecip(i) <= InfilCap(i) Infiltration(i) = ObsPrecip(i); Runoff(i) = 0;

else ,Infiltration(i) = InfilCap(i) + ((1-RunofFract(i)) .* ...

(ObsPrecip(i)-InfilCap(i)));

Runoff(i) = ObsPrecip(i) - Infiltration(i); Runoff(i) = max(Runoff(i),0);

end

WaterContent(i) = WaterContent(i) + mmToWC(Infiltration(i),... % Convert to fraction

SoilDepth(i));

if WaterContent(i) > Porosity(i)

ExtraRunoff(i) = (WaterContent(i) - Porosity(i)) .* ... % Convert to [mm]

SoilDepth(i) .* 1000;

Runoff(i) = Runoff(i) + ExtraRunoff(i);

Infiltration(i) = Infiltration(i) - ExtraRunoff(i); WaterContent(i) = Porosity(i);

end

(30)

30

%%% E V A P O R A T I O N %%%

% Calc. wind function (Fu) Penmann's wind function [W/m2/mbar]

Fu = 7.4*(0.5+0.54*ObsWind(Time));

Psych(i) = CP * Pressure(i) / (LATENT * EPSILON);

% (Cp * P) / (L* Eps)

S(i) = CalcEsat_mbar(Temperature(i)+0.5) - CalcEsat_mbar...

(Temperature(i)-0.5);

PotEvap(i) = (S(i)*NetRad(i) + Psych(i)*Fu*Esat(i)*(1-...

RelHum(i)/100)) / (S(i)+Psych(i));

PotEvap(i) = PotEvap(i) / LATENT * 24 * 3600;

% Convert to mm/day

PotEvap(i) = PotEvap(i) * CropFactor(i);

% Include Crop factor

if SurfType(i) == 1 Evap(i) = 0; elseif SurfType(i) == 2 Evap(i) = PotEvap(i);

else ,Evap(i) = PotEvap(i) * (WaterContent(i)-Wilting(i)) /...

(FieldCap(i)-Wilting(i));

Evap(i) = min(PotEvap(i),max(Evap(i),0));

WaterContent(i) = WaterContent(i) - mmToWC(Evap(i),... % Convert to fraction SoilDepth(i)); end end %%% D R A I N A G E %%%

% Calculates the drainage of soil water through the soil matrix

% if the soil water content exceeds the field capacity

for i = 1:length(Count) if i <= Lake

Drainage(i) = 0;

elseif WaterContent(i) < FieldCap (i)

Drainage(i) = 0;

else ,Drainage(i) = (WaterContent(i) - FieldCap(i)) / ...

DrainConst(i);

Drainage(i) = WCTomm(Drainage(i),SoilDepth(i)); end

WaterContent(i) = WaterContent(i) - mmToWC(Drainage(i),...

SoilDepth(i));

end

%%% W A T E R B A L A N C E %%%

% Calculates the water balance of the lake and new lake level

TotalDrainage = 0; TotalRunOff = 0;

for i = Soil

TotalDrainage = TotalDrainage + Count(i) * Drainage (i); TotalRunOff = TotalRunOff + Count(i) * Runoff(i);

(31)

31

Cellm2 = CellKm2 * 1.0e6;

% Convert km2 --> m2

LakeSurface = Cellm2 * Count(Lake);

% Surface area of the lake [m2]

TotalDrainage = TotalDrainage * Cellm2 / 1000;

% Now units are [m3]

TotalRunOff = TotalRunOff * Cellm2 / 1000;

% Now units are [m3]

Drainage_mm = TotalDrainage / LakeSurface * 1000;

% And now it is [mm]

Runoff_mm = TotalRunOff / LakeSurface * 1000;

% And now it is [mm]

Change_mm= ObsPrecip(Lake) - Evap(Lake) + Drainage_mm +...

Runoff_mm;

LakeLevel = LakeLevel + Change_mm / 1000;

% Lake level in [m]

LakeElev = LakeElev + Change_mm / 1000;

% Lake level in [m]

%%% D I S C H A R G E %%%

Discharge = 0.0;

Depth = LakeLevel - Bottom; if Depth > 0.0

Discharge = DischargeCoeff * power(Depth, 2.5);

% Discharge in [mm]

LakeLevel = LakeLevel - Discharge / 1000;

% Convert [mm] to [m]

LakeElev = LakeElev - Discharge / 1000; end %%% S T O R E R E S U L T S %%% StorePrecip(Time) = ObsPrecip(Lake); StoreEvap(Time) = Evap(Lake); StoreNetRad(Time) = NetRad(Lake); StoreRunoff(Time) = Runoff_mm; StoreDrainage(Time) = Drainage_mm; StoreDischarge(Time) = Discharge; StoreLevel(Time) = LakeLevel; Time = Time + dt; end % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % E N D O F D Y N A M I C P A R T % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %%% C O R R E L A T I O N %%% corDischarge = corr(StoreDischarge(1:EndTime),Discharge2(1:EndTime)) corRunoff = corr(StoreRunoff(1:EndTime),Runof2(1:EndTime)) corDrainage = corr(StoreDrainage(1:EndTime),Drain2(1:EndTime)) corLevel = corr(StoreLevel(1:EndTime),Level2(1:EndTime)) corEvap = corr(StoreEvap(1:EndTime),Evap2(1:EndTime)) corPrecip = corr(StorePrecip(1:EndTime),Precip2(1:EndTime))

(32)

32

corNetRad = corr(round(StoreNetRad(1:EndTime)),NetR2(1:EndTime))

ObsLevel = ObsLakeLevel ./ 100;

% Convert from cm to m

R_Obs_Mod = corr(StoreLevel(1:EndTime),ObsLevel(1:EndTime))

% Cross correlation of observerd and modelled lake levels

R2_SEE_Obs_Mod = fitlm(StoreLevel(1:EndTime),ObsLevel(1:EndTime))

% R-squared and SSE of observed and modelled lake levels

%%% V I S U A L I Z A T I O N %%% figure(1) plot(StoreLevel(1:EndTime),'k') hold on plot(ObsLevel(1:EndTime),'r') figure(2) subplot(4,1,1) plot(StoreLevel(1:EndTime),'k') hold on plot(Level2(1:EndTime),'r') title('Lake Level')

subplot(4,1,2) plot(StoreDrainage(1:EndTime),'k') hold on plot(Drain2(1:EndTime),'r') title('Drainage') subplot(4,1,3) plot(StoreDischarge(1:EndTime),'k') hold on plot(Discharge2(1:EndTime),'r') title('Discharge') subplot(4,1,4) plot(StoreEvap(1:EndTime),'k') hold on plot(Evap2(1:EndTime),'r') title('Evaporation')

Rising variance, changing skewness and critical slowing down: Early indicators of a tipping point in Andean Lake La Cocha

UNIVERSITY OF AMSTERDAM

BACHELOR THESIS