A Hydrological Model for Reconstructing Past Precipitation from Lake Level Records

Electronic SupplementaryMaterial at

Reconstructing past precipitation from lake levels and inverse modelling for Andean Lake La Cocha

J.H. van Boxel, Z. González-Carranza, H. Hooghiemstra, M. Bierkens and M.I. Vélez.

The paper is about reconstructing past precipitation by inverse modelling using a hydrological model. This supplementary materials contains a brief description of the equations used in the hydrological model, details on the model sensitivity test and description ot the methods used to estimate response times of the lake level.

Supplementary materials 1: Global description of the model

The hydrological model calculates the water level in Lake La Cocha as a function of daily values of precipitation, temperature, wind speed, sunshine duration, and relative humidity. The model includes processes such as precipitation, evaporation, soil water storage, drainage, overland flow and discharge (Fig. S01). Modelling of some of these processes is described in more detail below. The integration time step is one day and the grid cells are 100m x 100m.

Figure S01:Main processes incorporated in the hydrological model.

The water storage capacity of the soil is determined by soil depth, porosity, field capacity and wilting point. These and other parameters of the three different surfaces used in the model are defined in Table S01.

Table S01:Model parameters for the three surface types used in the model.

Supplementary materials 2: Evaporation

The evaporation of the lake is taken equal to Penman’s open water evaporation. Potential evaporation over land surfaces is calculated by multiplying open water evaporation (E0)with a crop factor (Table S01). Actual evaporation is calculated from potential evaporation by multiplying it with a drought factor (Fd) that depends on the soil water content. Fd = 1.0 if the soil water content greater or equal to field capacity; Fd = 0.0 if the soil water content is less or equal to the wilting point and between field capacity and the wilting point Fd is calculated by linear interpolation between those two points.

So actual evaporation over land (Ea) is calculated as:

[kg m-2 s-1] (S.1).

Open water evaporation (E0) is calculated using Penman’s (1948) equation:

[kg m-2 s-1] (S.2),

where s is the slope of the saturation vapour curve [mbar K-1], Q is net radiation [W m-2], γ is the psychrometric constant [mbar K-1], es is the saturation vapour pressure at a height of 2m [mbar], e2 is the real water vapour pressure at 2m elevation [mbar] and L is the latent heat of vaporization for water [J kg-1].

The wind function f(u) is give by (Penman 1948):

[W m-2 mbar-1] (S.3),

with u2 being the wind speed 2 m above the surface [m s-1].

The psychrometric constant (γ) is not constant at all, since it includes air pressure, which in the mountains depends on altitude:

[mbar K-1] (S.4),

Where P is the barometric pressure [mbar], cp is the heat capacity of dry air at constant pressure [J kg-1 K-1], Md is the molecular weight of dry air [kg kmol-1] and Mw is the molecular weight of water [kg kmol-1].

Note: instead of mbar you can also write hPa.

Supplementary materials 3: Radiation modelling

Net radiation (Q) is the global radiation (K↓) minus reflected short wave radiation (K↑) plus atmospheric radiation (L↓) minus the long wave radiation emitted by the surface(L↑):

[W m-2] (S.5).

To calculate global radiation, we first calculate the amount of solar radiation entering the atmosphere (potential radiation) and multiply that with the transmissivity of the atmosphere (using Linke’s turbidity factor). For the clouded hours the transmissivity of the clouds is also taken into account.

The daily average of potential radiation (Ip)is calculated from:

[W m-2] (S.6).

where I0 is the solar constant (1367 W m-2), F is an orbital correction, ω is the angular velocity of the earth (7.27x10-5rad s-1), φ is the latitude, δ is the solar declination and t0 is half the length of the day in seconds (time of sun set).

The orbital parameter (F), correcting for the varying distance between the sun and Earth is approximated by:

[-] (S.7),

where J is the Julian day number.

The declination (δ = the latitude where the sun is exactly in the zenith) varies with the season:

[º](S.8).

The time of sunset (half day length) can be calculated from:

[s](S.9).

For the hours that the sky was clear global radiation was calculated as the sum of direct radiation (unaffected by the atmosphere) and diffuse radiation (scattered by air molecules and aerosols).

[W m-2](s.10).

Direct radiation was calculated from the potential radiation using Linke’s turbidity factor:

[Wm-2] (S.11),

Where TR is Rayleigh transmission and T is Linke’s turbidity factor, which is set at 2.7 (a value for clean air).

Rayleigh transmission (TR) was calculated from the optical mass (M) and approximated by a 4th order polynomial:

[-](S.12).

The Rayleigh transmission of the cloudless sky varies with solar elevation. The transmission at 3:00 hours (solar time) was taking as being representative for the daily average transmission. Therefore the optical mass was calculated for 3:00 solar time from:

[-](S.13),

Where P is the atmospheric pressure in mbar, which is a function of elevation and cos(Z) is the cosine of the zenith angle calculated for 3:00 solar time (t = 10800 s) from:

[-](S.14).

Diffuse radiation was calculated from an approximation of Erbs formula (Erbs et al. 1982), derived by Van Boxel (2002):

(S.15),

With:[-](S.16).

When the sky is clouded all radiation reaching the earth surface is diffuse radiation. In this case global radiation was calculated by multiplying the potential radiation with a mean transparency of the clouds, which was estimated at 0.40, following Van Boxel (2002).

[W m-2](S.17).

The reflected short wave radiation is calculated from global radiation and surface albedo (α):

[W m-2](S.18).

Long wave radiation emitted by the surface was calculated from the temperature (T) using the Stefan-Boltzmann equation and an emissivity of 0.98:

[W m-2](S.19),

where T is the absolute temperature of the surface and σ is the Stefan-Boltzmann constant (5.67x10-8 W m-2 K-4).

The atmospheric radiation was also calculated with the Stefan-Boltzmann equation, but now using an apparent emissivity of the atmosphere (εa), which is a function of cloud cover (CC) and water vapour pressure (e), using Brunt’s (1932) approximation for the clear sky:

[W m-2](S.20),

with:

[-](S.21),

Supplementary materials 4: Sunshine relative humidity in paleoconditions

When calibrating the model under modern conditions daily mean/total values for relative humidity and sunshine duration were available. However these are not available for the past 14000 years and we do not have a suitable proxy for these parameters.

The model is not sensitive to sunshine duration and relative humidity, but nevertheless these parameters were estimated as a function of mean annual precipitation (MAP) in mm.

We assumed the relative humidity increases to 100% and the sky is completely clouded (i.e. a sunshine duration of zero) if mean MAP is extremely high. If, on the other hand, MAP is zero it is assumed that relative humidity is 0% and sunshine duration is 12 hours (no clouds). Between these values sunshine duration (S) and relative humidity (R) are interpolated using an exponential relation that also satisfies the modern climatological conditions (MAP=1384 mm, S=2.5 hour day-1 and R=86.6%):

[%] (S.22),

[hr day-1] (S.23).

These relations are graphically depicted in figure S02.

Figure S02:Sunshine duration and relative humidity estimated by precipitation. The blue and red circles indicate modern conditions.

Supplementary materials 5: The power in the river discharge formula

In order to correctly model lake-level dynamics it is necessary to calibrate how the discharge responds to the lake level. The discharge was described as:

[m3s-1](S.24),

where D is the discharge, Hl is the lake level, b is the level of the bottom of the river and a is a proportionality constant (in m0.5 s-1) and p is the power of the relation.

A power of 1.5 would correspond to a rectangular cross section of the river (the width of the river does not increase if the level rises) and a power of 2.5 would correspond to a more or less triangular cross section of the river, comparable to a V-notch (ISO 1980; ASTM 1993).

We tried to find the best value for the power by ranging the power (p) from 1.5 to 2.5 and optimizing the parameters a and b for each value of the power. The standard error of estimate of the simulated lake levels was used as an indicator for the goodness of fit.

The standard error of estimate ranged from 0.146 m (p=1.5) to 0.149 m (p=2.5). This difference is so small that on the basis of the fit it could it cannot be decided which power is correct. Since it is more likely that the width of the river changes if the lake level rises we decided to use the power of 2.5.

Supplementary materials 6: Model sensitivity

Model sensitivity was tested for model parameters (precipitation correction PC and drainage time constant τd) and meteorological variables (precipitation, temperature, relative humidity, sunshine durarion and wind speed).

Figure S03:Standard error of estimate and correlation between modelled and observed lake levels for a range of values of the precipitation correction and drainage time constant (τd).

Because of the rather humid conditions in the area the soil water content was rarely limiting evaporation. This resulted in the model being insensitive to model parameters such as soil depth, porosity, field capacity and wilting point. The model was further tested for its sensitivity to the precipitation correction (PC = ratio between areal mean precipitation and precipitation observed at the meteorological station) and the drainage time constant (Fig. S03).

If PC increases the standard error of estimate also increases, and the correlation decreases, indicating that the fit gets worse as PC gets larger. The fit was relatively insensitive to the drainage time constant, so it was decided to keep this at its initial estimate of 5 days. For the range of values tested a PC of 1.0 yielded the best fit with a standard error of estimate of 0.15 m for the modelled water level and a correlation between modelled and observed water levels of 0.87. The optimal discharge parameters for PC=1.0 were a=1.364 m0.5s-1 and b=-0.050 m, yielding a mean annual discharge of 3.57 m3s-1.

The meteorological variables for which model sensitivity was tested are: precipitation, temperature, relative humidity, sunshine duration and wind speed. The model was run until equilibrium, changing one meteorological variable over a wide range, keeping everything else constant (Fig.S04).

Modelled lake levels and discharge were most sensitive to changes in precipitation and temperature (Fig. S04).

If MAP was varied from 1000 to 2000 mm the lake level varied from 0.80 m to 1.91 m, so the range in modelled lake levels is 1.11 m (=1.91m–0.80m). The response is not linear (R2=0.971 for a linear fit) because the lake surface increases at increasing lake levels. A MAP of 1000 mm yielded a discharge of only 0.9 m3s-1, whereas 2000 mm gave a discharge of 7.5 m3s-1; a range of 6.6 m3s-1. The response of the discharge to precipitation is linear (R2=0.99999).

The main effect of higher MAT (mean annual temperature) is that it increases evaporation and the lake level and discharge decrease as a consequence. Keeping MAP constant at 1380 mm and varying MAT from 5° to 15 ºC yielded lake levels from 1.59 m till 1.28 m and discharge values from 4.8 m3s-1 till 2.8 m3s-1. Both lake level and discharge respond almost linearly to temperature change (R2=0.97 and 0.98).

Although the other parameters (relative humidity, sunshine duration and wind speed) were varied in a fairly wide range, changes in lake level and discharge were less sensitive to these parameters and in all cases the response was linear (R20.99).

Supplementary materials 7: Calculating the response time of the lake system

The response time of the lake level to changes in precipitation depends on whether the lake level is above or below 0.0 m. (Fig. 10b). At a lake level above 0.0 m precipitation estimates are mainly determined by the lake level; at lower lake levels precipitation estimates by inverse modelling are mainly sensitive to temperature.

During the last 6 months of 1997 the lake level dropped from ~2.0 m to ~1.0 m (Fig. 6), equivalent to a discharge of 42 x 106 m3 of water. At a lake level of 2.0 m the discharge is 8.2 m3s-1. The calculated response time (the water volume to be removed divided by discharge) is 59 days. Since discharge is not linearly related to lake level, the response time will depend on the lake level (with shorter response times at higher lake levels), but two months is a good estimate for the order of magnitude of the response time.

The response time was also tested with the hydrologicalmodel. The model was fed with a constant precipitation if 1380 mm yr-1, yielding a lake level of 1.38 m. Then the lake level was increased or decreased by 0.5 m keeping precipitation constant. The modelled lake level went back to 1.38 m and the dynamics could be described by a simple exponential decay (first order system). Time constants were 95 days for the decreased lake level and 75 days for the increased lake level.

We can also estimate the response time (τ) from Fig. 10b as the ratio of equilibrium lake-level change (dL) over change in required MAP (dP), corrected with the ratio of the surface of the catchment (AC) and the lake (AL):

[s] (S.25)

At 11.0 °C the required MAP for an equilibrium lake level of 2.0 m is 2386 mmyr-1 and for an equilibrium lake level of 1.0 m a MAP of 1175 mm yr-1is needed. Using these numbers we arrive at an estimated response time of 0.16 year, which also equals 2 months.

In the modern situation, with a river discharging the lake, there is a strong feedback between lake level and discharge, so the lake can respond quickly to changes in the environmental conditions. Response times are in the order of 2 to 3 months, also depending on the lake level. This response time is a measure for the time the lake needs to return to its equilibrium level. However a rain storm can increase the lake level almost instantly.

If the lake is not discharging strong feedback between lake level and discharge disappears and the response becomes slower. If we apply eq. S.25 to estimate the response time for lake levels below 0 m it is longer. For example, at 11 °C we read 930 mm yr-1at a lake level of 0 m, and 915 mm yr-1at a lake level of -15 m. This yields a response time in the order of 173 years. If we do the response test with the model around a lake level of -2 m we obtain a response time of 65 years if we increase the lake level by 0.5 m, and 135 years if we decrease it by 0.5 m. So under constant climatological conditions, after disturbance it will take a long time before the lake level returns to equilibrium. However if climatological conditions are not constant the lake level can respond quickly. A wet year or period will immediately result in a higher lake level and a dry year will result in a lower lake level.

The pollen-based environmental reconstruction showed that quasi-stable conditions lasted for 200 to 600 years (González-Carranza et al. 2012). This means that even with a response time in the order of 100 years the lake can be assumed to be in equilibrium.

Supplementary materials references:

ASTM (1993)Standard method for open-channel flow measurement of water with thin-plate weirs. American Society for Testing and Materials, document nr. ASTM D5242.

Brunt D (1932) Notes on radiation in the atmosphere. Quarterly Journal of the Royal Meteorological Society 58: 389-420.

Erbs DG, Klein SA, Duffie JA (1982)Estimation of the diffuse radiation fraction for hourly, daily and monthly-average global radiation. Solar Energy28: 293–302.

ISO (1980) Water flow measurement in open channels using weirs and venturi flumes - Part 1: Thin plate weirs.International Organization of Standards, ISO 1438.

Penman HL (1948) Natural evaporation from open water, bare soil and grass. Proceedings Royal Society London Series A193: 120-146.

Van Boxel JH (2002) Modelling global radiation for the Portofino area in Italy. Report of the Institute for Biodiversity and Ecosystem Dynamics (IBED), University of Amsterdam, The Netherlands, 21 pp. (

- 1 -