Quantifying energy and water fluxes in dry dune ecosystems of the Netherlands

Abstract

Coastal and inland dunes provide various ecosystem services that are related to groundwater, such as drinking water production and biodiversity. To manage groundwater in a sustainable manner, knowledge of actual evapotranspiration (ET${}_{\mathrm{a}})$ for the various land covers in dunes is essential. Aiming at improving the parameterization of dune vegetation in hydrometeorological models, this study explores the magnitude of energy and water fluxes in an inland dune ecosystem in the Netherlands. Hydrometeorological measurements were used to parameterize the Penman–Monteith evapotranspiration model for four different surfaces: bare sand, moss, grass and heather. We found that the net longwave radiation ($R_{\text{nl}}$) was the largest energy flux for most surfaces during daytime. However, modeling this flux by a calibrated FAO-56 $R_{\text{nl}}$ model for each surface and for hourly time steps was unsuccessful. Our $R_{\text{nl}}$ model, with a novel submodel using solar elevation angle and air temperature to describe the diurnal pattern in radiative surface temperature, improved $R_{\text{nl}}$ simulations considerably. Model simulations of evaporation from moss surfaces showed that the modulating effect of mosses on the water balance is species-dependent. We demonstrate that dense moss carpets (Campylopus introflexus) evaporate more (5 %, $+$14 $\mathrm{mm}$) than bare sand (total of 258 $\mathrm{mm}$ in 2013), while more open-structured mosses (Hypnum cupressiforme) evaporate less ($-$30 %, $-$76 $\mathrm{mm}$) than bare sand. Additionally, we found that a drought event in the summer of 2013 showed a pronounced delayed signal on lysimeter measurements of ET${}_{\mathrm{a}}$ for the grass and heather surfaces, respectively. Due to the desiccation of leaves after the drought event, and their feedback on the surface resistance, the potential evapotranspiration in the year 2013 dropped by 9 % ($-$37 $\mathrm{mm}$) and 10 % ($-$61 $\mathrm{mm}$) for the grass and heather surfaces, respectively, which subsequently led to lowered ET${}_{\mathrm{a}}$ of 8 % ($-$29 $\mathrm{mm}$) and 7 % ($-$29 $\mathrm{mm}$). These feedbacks are of importance for water resources, especially during a changing climate with an increasing number of drought days. Therefore, such feedbacks need to be integrated into a coupled plant physiological and hydrometeorological model to accurately simulate ET${}_{\mathrm{a}}$. In addition, our study showed that groundwater recharge in dunes can be increased considerably by promoting moss vegetation, especially of open-structured moss species.

Introduction

Coastal and inland sand dunes are major drinking water production sites in the Netherlands. Approximately 23 % of Dutch drinking water originates from aquifers in these dunes, which are replenished by both natural groundwater recharge and artificial infiltration of surface waters. Another ecosystem service of groundwater in dune systems is that shallow groundwater tables sustain nature targets with a very high conservation value. Such targets, like wet dune slacks and oligotrophic pools, are often legally enforced, e.g., by the European Habitat Directive and by the Water Framework Directive. Furthermore, a deep layer of fresh groundwater in coastal dunes protects the hinterland from the inflow of saline groundwater.

Under a warming climate, summers are expected to become dryer and the water quality of surface waters may degrade (Delpla et al., 2009), especially during dry periods with low river discharge rates (Zwolsman and van Bokhoven, 2007; van Vliet and Zwolsman, 2008). To maintain current drinking water quality and production costs, water production in the future may have to rely more on natural groundwater recharge. This implies that drinking water companies need to search for new water production sites or intensify current groundwater extractions, while protecting groundwater-dependent nature targets.

For sustainable management of renewable groundwater resources, groundwater extractions should be balanced with the amount of precipitation that percolates to the saturated zone, the groundwater recharge. Knowledge of actual evapotranspiration (ET${}_{\mathrm{a}}$, here defined as the sum of plant transpiration, soil evaporation and evaporation from canopy interception) for the various land covers is essential to quantify the amount of recharge. Inland dune systems are predominantly covered with deciduous and pine forest. Well-developed hydrometeorological models are available to simulate ET${}_{\mathrm{a}}$ for these forest ecosystems (Dolman, 1987; Moors, 2012). Other ecosystems, such as heathland and bare sand colonized by algae, mosses, tussock-forming grasses or lichens, received less attention. However, heathland and drift sand ecosystems have a higher conservation value than forest plantations, in particular those of coniferous trees. Nature managers are therefore often obligated to protect and develop certain heathland and drift sand ecosystems at the expense of forest ecosystems (The European Natura 2000 policy). A better parameterization of heathland and drift sand ecosystems in hydrometeorological models would aid in the sustainable management of important groundwater resources and would allow quantifying the cost and benefit of nature conservation in terms of groundwater recharge.

To this end, this study explores diurnal patterns in energy and water fluxes in a dry dune ecosystem on an elevated sandy soil in the Netherlands. Our study aims to improve the parameterization of dune vegetation in hydrometeorological models based on field measurements, focusing on four different surfaces: bare sand, moss (Campylopus introflexus), grass (Agrostis vinealis) and heather (Calluna vulgaris). A second objective is to quantify the effect of moss species on the water balance. Mosses and lichens are present in most successional stages in dry dune ecosystems, either as pioneer species or as understory vegetation. Voortman et al. (2014) hypothesized that moss-covered soils could evaporate less than a bare soil, since the unsaturated hydraulic properties of moss layers reduce evaporation under relatively moist conditions. Such hydraulic behavior could have large implications on the ecological interactions between vascular and nonvascular plants in water-limited ecosystems, as the presence of a moss cover could facilitate the water availability for rooting plants. Such interactions are of importance to groundwater resources, as the resilience of plant communities to drought determines the succession rate and biomass, which subsequently feedback on evapotranspiration.

Organization of the research from measurements to model simulations.

[Figure omitted. See PDF]

A third objective is to gain insight into the delayed effect of dry spells on potential and actual evapotranspiration for heathlands and grasslands. To quantify the evapotranspiration loss term, many hydrological modeling frameworks use the concept of potential evapotranspiration ET${}_{\mathrm{p}}$ (Federer et al., 1996; Kay et al., 2013; Zhou et al., 2006), defined as the maximum rate of evapotranspiration from a surface where water is not a limiting factor (Shuttleworth, 2007). ET${}_{\mathrm{p}}$ is input to modeling frameworks and reduces to ET${}_{\mathrm{a}}$ in cases of water stress. However, if dry spells result in a vegetation dieback, the simulated ET${}_{\mathrm{p}}$ should be adjusted to account for the smaller transpiring leaf area after the dry spell. The model simulations presented in this paper give some guidance on the magnitude of errors in simulated ET${}_{\mathrm{a}}$ if feedbacks of dry spells on ET${}_{\mathrm{p}}$ are neglected.

The knowledge presented in this paper will help to improve and interpret the simulations of water recharge in sand dunes by hydrological models, and will sustain rainwater harvesting in dunes by vegetation management.

Measurements and methods

General setup

A field campaign started in August 2012 to measure energy and water fluxes in the drinking water supply area “Soestduinen”, situated on an elevated sandy soil (an ice-pushed ridge) in the center of the Netherlands (52.14${}^{\circ }$ latitude, 5.31${}^{\circ }$ longitude). Due to deep groundwater levels, the vegetation in this region is groundwater-independent, i.e., relying solely on rainwater (on average 822 mm rain per year, 40 % falling in the first 6 months of the year and 60 % falling in the last 6 months of the year). The reference evapotranspiration according to Makkink (1957) is on average 561 $\mathrm{mm}$ per year. The field data were used to parameterize the Penman–Monteith equation, to calculate ET${}_{\mathrm{p}}$ and to perform hydrological model simulations of ET${}_{\mathrm{a}}$, based on the actual availability of soil moisture. The Penman–Monteith equation is given by $${\text{ET}}_{\mathrm{p}}=\frac{\mathrm{\Delta }\left(R_{\mathrm{n}},-,G\right)+{\mathit{\rho }}_{\mathrm{a}}{c}_{\mathrm{p}}\left(e_{\mathrm{s}},-,e_{\mathrm{a}}\right)/r_{\mathrm{a}}}{\left(\mathrm{\Delta },+,\mathit{\gamma },\left(\mathrm{1},+,\frac{r_{\mathrm{s}}}{r_{\mathrm{a}}}\right)\right)\mathit{\lambda }{\mathit{\rho }}_{\mathrm{w}}},$$ where ET${}_{\mathrm{p}}$ is the potential evapotranspiration ($\mathrm{mm}{\mathrm{s}}^{-\mathrm{1}}$), $\mathrm{\Delta }$ is the slope of the saturation vapor pressure vs. temperature curve ($\mathrm{kPa}{}^{\circ }{\mathrm{C}}^{-\mathrm{1}}$), $R_{\mathrm{n}}$ is the net radiation ($\mathrm{J}{\mathrm{m}}^{-\mathrm{2}}$), $G$ is the soil heat flux ($\mathrm{J}{\mathrm{m}}^{-\mathrm{2}}$), ${\mathit{\rho }}_{\mathrm{a}}$ is the air density ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$), $c_{\mathrm{p}}$ is specific heat of moist air ($\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{}^{\circ }{\mathrm{C}}^{-\mathrm{1}}$), $e_{\mathrm{s}}$ is the saturation vapor pressure of the air (kPa), $e_{\mathrm{a}}$ is the actual vapor pressure of the air ($\mathrm{kPa}$), $r_{\mathrm{a}}$ is aerodynamic resistance to turbulent heat and vapor transfer ($\mathrm{s}{\mathrm{m}}^{-\mathrm{1}}$), $\mathit{\gamma }$ is the psychrometric constant ($\mathrm{kPa}{}^{\circ }{\mathrm{C}}^{-\mathrm{1}}$), $\mathit{\lambda }$ is the latent heat of vaporization ($\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}$) and ${\mathit{\rho }}_{\mathrm{w}}$ is the density of liquid water ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$). Results of Irmak et al. (2005) suggest that estimates of ET${}_{\mathrm{p}}$ on hourly time steps are more accurate than estimates on a daily timescale. Furthermore, Liu et al. (2005) showed that the use of daily input values leads to a systematic overestimation of ET${}_{\mathrm{a}}$, especially for sandy soils. Hence, energy fluxes in the Penman–Monteith equation are preferably simulated at subdiurnal timescales. Furthermore, understanding and simulation of plant physiological processes requires knowledge of the diurnal variation of environmental variables (Nozue and Maloof, 2006). Therefore, field data were aggregated to hourly time steps to maintain the diurnal pattern and to analyze our field results at the same time interval as commonly available climate data.

In this paper evapotranspiration is defined as the sum of transpiration, soil evaporation and evaporation from canopy interception, expressed in mm per time unit. Radiative and soil heat fluxes are expressed in $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$. Figure 1 shows the procedures followed to translate field data (Sect. 2.1) to submodels of the Penman–Monteith equation (Sect. 2.2) and to subsequently calculate ET${}_{\mathrm{p}}$ and simulate ET${}_{\mathrm{a}}$ (Sect. 2.3).

Hydrometeorological measurements

Four homogeneous sites of bare sand, moss (Campylopus introflexus), grass (Agrostis vinealis) and heather (Calluna vulgaris) (Fig. 2) were selected to measure actual evapotranspiration (ET${}_{\mathrm{a}})$, the net radiation ($R_{\mathrm{n}}$), the soil heat flux ($G$) and the albedo. Other meteorological variables such as wind speed ($u$, at 2 $\mathrm{m}$ above the surface), relative humidity (RH, 1.5 $\mathrm{m}$ above the surface), air temperature ($T_{\mathrm{a}}$, 1.5 $\mathrm{m}$ above the surface) and rain ($P$) were measured at a weather station, installed in-between the measurement plots at a maximum distance of 40 $\mathrm{m}$ from each plot. Measurements were collected with data loggers (CR1000, Campbell Scientific Inc.) at a 10 s interval and aggregated to minutely values. Field measurements of bare sand, moss and grass were collected between August 2012 and November 2013. The field measurements in the heather vegetation were collected between June 2013 and November 2013.

The vegetation types studied in this paper, (a) the moss surface with an approximately 2 cm thick layer of Campylopus introflexus (inset), (b) the grass surface, primarily Agrostis vinealis and (c) the heather surface, Calluna vulgaris.

[Figure omitted. See PDF]

Measured incoming solar radiation $R_{\mathrm{s}\downarrow }$ at the four different surfaces. Periods with snow cover or sensor maintenance were omitted.

[Figure omitted. See PDF]

The net radiation was measured with net radiometers (NRLite2 Kip & Zonen B.V.). The net radiometers were installed at a relatively low height of 32, 40, 40 and 50 $\mathrm{cm}$ above the bare sand, moss, grass and heather surfaces, respectively (relative to the average vegetation height), to limit the field of view to a homogenous surface. The incoming solar radiation ($R_{\mathrm{s}\downarrow }$) and reflected solar radiation ($R_{\mathrm{s}\uparrow }$) were measured with an albedo meter (CMA6, Kip & Zonen B.V.) that was rotated between the four surfaces. It was installed next to each $R_{\mathrm{n}}$ sensor. Due to a snow cover (winter months) or sensor maintenance (October 2012, May 2013), some periods were omitted (Fig. 3).

Eight self-calibrating heat flux plates (HFP01SC, Hukseflux B.V.) (two for each site) were installed 8 $\mathrm{cm}$ below the soil surface near the net radiometers. These heat flux plates were programmed to calibrate themselves for 15 $\min$ at 6 h time intervals, based on a known heat flux supplied by an integrated heater. Besides each soil heat flux plate, an averaging thermocouple (TCAV, Campbell Scientific, Inc.) was installed at 2 and 6 cm depth and a soil moisture probe (CS616, Campbell Scientific, Inc) was installed at 4 cm depth to estimate the change in heat storage ($S$) above the heat flux plates. The sum of the measured soil heat flux at 8 cm depth and $S$ represents the heat flux at the soil surface. Sensor installation and procedures to calculate $S$ were followed according to the Campbell Scientific Inc. (2014) HFP01SC instruction manual.

Lysimeter design.

[Figure omitted. See PDF]

Within each surface, one weighing lysimeter was installed. The lysimeters (Fig. 4) had a 47.5 cm inner diameter and were 50 cm deep. Intact soil monoliths were sampled by hammering the PVC tube into the soil, alternated with excavating the surrounding soil to offset soil pressures. The lysimeters were turned upside down, to level the soil underneath and to close this surface with a PVC end cap. To allow water to drain out of the lysimeter bottom plate, a 2.5 cm diameter hole was made in the base plate. A 15 cm long fiberglass wick (Pepperell $\mathrm{2}\times \mathrm{1}/\mathrm{2}$ inch) was installed in the PVC end cap to guide drainage water through the hole into a tipping bucket (Davis 7852) below the lysimeter. The wick, together with two sheets of filter cloth (140–150 $\mathrm{\mu }\mathrm{m}$, Eijkelkamp Agrisearch Equipment), placed at the bottom of the lysimeter tank, prevented soil particles from flushing out of the lysimeter. The tipping bucket below the lysimeter had a resolution of 0.2 $\mathrm{mm}$ for the intercepting area of the tipping bucket, which was equal to 0.024 $\mathrm{mm}$ for the cross-sectional area of the lysimeter. Drainage water was collected in a reservoir installed below the lysimeter.

The lysimeters were weighted with temperature compensated single point load cells (Utilcell 190i, max 200 $\mathrm{kg}$). These load cells were initially connected to the full bridge data ports of the data loggers. However, the measurement resolution of the data loggers was too coarse to fully compensate for temperature effects on weight measurements. Fluctuations of 0.333 $\mathrm{\mu }\mathrm{V}$ due to temperature effects were within the data logger measurement resolution, which equals 36 $\mathrm{g}$ in weight change, i.e., 0.2 $\mathrm{mm}$ of evaporation. To increase the lysimeter precision, digitizers (Flintec LDU 68.1) were installed in May 2013 to process and digitize the load cell signals without interference of the data logger. In this setup, a measurement resolution of 10 g was achieved, i.e., 0.06 mm equivalent water depth, which is adequate for measuring ET${}_{\mathrm{a}}$ for daily time periods (subtracting two values would lead to a maximum error of 0.06 $\mathrm{mm}$ caused by the measurement resolution). Analysis of measured ET${}_{\mathrm{a}}$ was therefore limited to the period after installation of the digitizers.

After a rain event on 7 September 2013, the tipping buckets below the grass and heather lysimeters became partly clogged with beetles nesting underneath the lysimeters. This led to a continuous drainage signal which was out of phase with the weight measurements. Without accurate drainage measurements, lysimeter weight signals cannot be transferred to evapotranspiration. Therefore, ET${}_{\mathrm{a}}$ data on days with a poor drainage signal after 7 September 2013 were disregarded in the analyses for the grass and heather lysimeters.

Parameterization of the Penman–Monteith equation

Net radiation ($R_{\mathrm{n}}$)

The net radiation ($R_{\mathrm{n}})$ is defined as $$R_{\mathrm{n}}=R_{\text{ns}}+R_{\text{nl}}=\left(\mathrm{1},-,\text{albedo}\right)R_{\mathrm{s}\downarrow }+\left({\mathit{\epsilon }}_{\mathrm{s}},R_{\mathrm{l}\downarrow },-,R_{\mathrm{l}\uparrow }\right),$$ where $R_{\text{ns}}$ is the net shortwave radiation, $R_{\text{nl}}$ is the net longwave radiation, $R_{\mathrm{s}\downarrow }$ is the incoming solar radiation, $R_{\mathrm{l}\downarrow }$ is the downwelling longwave radiation from the atmosphere to the surface, $R_{\mathrm{l}\uparrow }$ is the emitted longwave radiation by the surface into the atmosphere and ${\mathit{\epsilon }}_{\mathrm{s}}$ is the surface emissivity representing the reflected downwelling longwave radiation. The albedo in Eq. (2) was determined by linear regression between measured $R_{\mathrm{s}\downarrow }$ and $R_{\mathrm{s}\uparrow }$. Based on the albedo obtained this way, $R_{\text{nl}}$ follows from measurements of $R_{\mathrm{n}}$ by subtracting calculated $R_{\text{ns}}$ from measured $R_{\mathrm{n}}$. Throughout this paper, this back-calculated $R_{\text{nl}}$ is referred to as the measured $R_{\text{nl}}$.

In hydrometeorological models, $R_{\text{nl}}$ is commonly estimated under clear sky conditions and multiplied by a factor to correct for clouds (Irmak et al., 2010; Gubler et al., 2012; Blonquist Jr. et al., 2010; Temesgen et al., 2007). A similar approach was followed in this study in which the Stefan–Boltzmann law is substituted into Eq. (2) for $R_{\mathrm{l}\downarrow }$ and $R_{\mathrm{l}\uparrow }$ under clear sky conditions (Saito and Šimůnek, 2009; Van Bavel and Hillel, 1976) and multiplied by a cloudiness function to obtain $R_{\text{nl}}$: $$R_{\text{nl}}=\left({\mathit{\epsilon }}_{\mathrm{s}},{\mathit{\epsilon }}_{\mathrm{a}},\mathit{\sigma },T_{\mathrm{a}}^{\mathrm{4}},-,{\mathit{\epsilon }}_{\mathrm{s}},\mathit{\sigma },T_{\mathrm{s}}^{\mathrm{4}}\right)f_{\text{cd}},$$ where ${\mathit{\epsilon }}_{\mathrm{a}}$ is the clear sky emissivity of the atmosphere (–), ${\mathit{\epsilon }}_{\mathrm{s}}$ is the surface emissivity (–), $\mathit{\sigma }$ is the Stefan–Boltzmann constant ($5.67\times {10}^{-\mathrm{8}}$ ${\mathrm{Wm}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{1}}$), $T_{\mathrm{a}}$ is the air temperature (K), $T_{\mathrm{s}}$ is the surface temperature (K) and $f_{\text{cd}}$ is a cloudiness function (–; described later). For vegetated surfaces, ${\mathit{\epsilon }}_{\mathrm{s}}=0.95$ was used (based on Jones, 2004), and ${\mathit{\epsilon }}_{\mathrm{s}}=0.925$ was used for bare sand (based on Fuchs and Tanner, 1968). Estimating ${\mathit{\epsilon }}_{\mathrm{a}}$ has a long history and numerous parameterizations are available. In this study, the empirical relationship found by Brunt (1932) was used: $${\mathit{\epsilon }}_{\mathrm{a}}=0.52+0.065\sqrt{e_{\mathrm{a}}},$$ where $e_{\mathrm{a}}$ is the water vapor pressure measured at screen level (hPa). The cloudiness function $f_{\text{cd}}$ in Eq. (3) is limited to $0.05\le f_{\text{cd}}\le \mathrm{1}$ and equal to $$f_{\text{cd}}=\frac{R_{\mathrm{s}\downarrow }}{R_{\text{s0}}},$$ where $R_{\text{s0}}$ is the estimated clear sky solar radiation. We estimated $R_{\text{s0}}$ following the FAO irrigation and drainage paper No. 56 (Allen et al., 1998). Since $f_{\text{cd}}$ is undefined during the night, an interpolation of $f_{\text{cd}}$ between sunset and sunrise is required. According to Gubler et al. (2012) $f_{\text{cd}}$ can be best linearly interpolated between the 4 to 6 h average before sunset and after sunrise. We adopted this approach, applying a 5 h average.

An estimate of $T_{\mathrm{s}}$ is required to fully parameterize Eq. (3). We developed a new approach to simulate the diurnal pattern in $T_{\mathrm{s}}$. Using Eq. (3), we back-calculated $T_{\mathrm{s}}-T_{\mathrm{a}}$ based on measured $R_{\text{nl}}$ for clear hours ($f_{\text{cd}}>$ 0.9). Generally, $T_{\mathrm{s}}-T_{\mathrm{a}}$ will be negative during nighttime (when solar elevation $\mathit{\beta }$ (radians) $<\mathrm{0}$), and will gradually increase to positive values during daytime ($\mathit{\beta }>\mathrm{0}$ 0). We describe this pattern by (Fig. 5): $$T_{\mathrm{s}}-T_{\mathrm{a}}=f_{\text{cum}}\left(\mathit{\beta },,,{\mathit{\mu }}_{\mathit{\beta }},,,{\mathit{\sigma }}_{\mathit{\beta }}\right)\left(T_{\text{s, amp}},+,\mathit{\beta },T_{\text{s, slope}}\right)+T_{\text{s, offset}},$$ where $f_{\text{cum}}$ is a cumulative normal distribution function with mean ${\mathit{\mu }}_{\mathit{\beta }}$ and standard deviation ${\mathit{\sigma }}_{\mathit{\beta }}$, describing the moment at which the surface becomes warmer than the air temperature (${\mathit{\mu }}_{\mathit{\beta }}$) and the speed at which the surface warms up or cools down (${\mathit{\sigma }}_{\mathit{\beta }}$) as a function of solar elevation angle ($\mathit{\beta }$). $T_{\text{s,amp}}$ is the amplitude of $T_{\mathrm{s}}$ (K), $T_{\text{s,slope}}$ is the slope between $\mathit{\beta }$ and $T_{\mathrm{s}}-T_{\mathrm{a}}$ during daytime ($\mathrm{K}{\mathrm{radians}}^{-\mathrm{1}}$) and $T_{\text{s,offset}}$ is the average value of $T_{\mathrm{s}}-T_{\mathrm{a}}$ during nighttime (K). The parameters of Eq. (6), except $T_{\text{s,offset}}$, were fitted to the data by minimizing the root mean square error (RMSE) by generalized reduced gradient nonlinear optimization. The $T_{\text{s,offset}}$ was determined as the average nighttime $T_{\mathrm{s}}-T_{\mathrm{a}}$ to limit the number of parameters during the optimization. Equation (6) was substituted for $T_{\mathrm{s}}$ in Eq. (3) to estimate $R_{\text{nl}}$. This novel approach to derive $R_{\text{nl}}$ was compared to the $R_{\text{nl}}$ model of the FAO-56 approach (Allen et al., 1998), originally derived to obtain daily estimates of $R_{\text{nl}}$ (using minimum and maximum daily $T_{\mathrm{a}}$ divided by 2 instead of $T_{\mathrm{a}}$ in Eq. 7) but commonly applied at hourly timescales (ASCE-EWRI, 2005; Perera et al., 2015; Gavilán et al., 2008; López-Urrea et al., 2006): $$R_{\text{nl}}=-\mathit{\sigma }{T}_{\mathrm{a}}^{\mathrm{4}}(a-b\sqrt{e_{\mathrm{a}}})\left(1.35,\frac{R_{\mathrm{s}}}{R_{\mathrm{s}\mathrm{0}}},-,0.35\right),$$ where the first term between brackets represents the net emittance, which should compensate for the fact that $T_{\mathrm{s}}$ is not measured. The empirical parameters $a$ and $b$ can be calibrated for a specific climate and/or vegetation. The second term between brackets is a cloudiness function. The default parameter values for $a$ and $b$ are 0.34 and 0.14, respectively (Allen et al., 1998). We calibrated these parameters for every site by linear least squares regression for clear days ($R_{\mathrm{s}}/R_{\text{s0}}>0.9$) and compared the performance of both $R_{\text{nl}}$ models (Eqs. 3, 7).

Equation (6) and associated parameters to describe the surface–air temperature difference, substituted for $T_{\mathrm{s}}$ in $R_{\text{nl}}$ (Eq. 3).

[Figure omitted. See PDF]

Soil heat flux ($G$)

The soil heat flux is commonly expressed as a fraction of $R_{\mathrm{n}}$, particularly on large scales using remote sensing (Su, 2002; Bastiaanssen et al., 1998; Kustas et al., 1998; Kustas and Daughtry, 1990; Friedl, 1996). We adopted the same approach, making a distinction between daytime ($F_{\text{day}}$) and nighttime ($F_{\text{night}}$) fractions, determined by linear least squares regression between $R_{\mathrm{n}}$ and the average of the two sets of soil heat flux measurements.

Aerodynamic resistance ($r_{\mathrm{a}}$)

The aerodynamic resistance under neutral stability conditions can be estimated by (Monteith and Unsworth, 1990) $$r_{\mathrm{a}}=\frac{\ln \left[\frac{z_{\mathrm{m}}-d}{z_{\text{om}}}\right]\ln \left[\frac{z_{\mathrm{h}}-d}{z_{\text{oh}}}\right]}{k^{\mathrm{2}}{u}_z},$$ where $z_{\mathrm{m}}$ is the height of wind speed measurements (m), $d$ is the zero plane displacement height (m), $z_{\text{om}}$ is the roughness length governing momentum transfer (m), $z_{\mathrm{h}}$ is the height of the humidity measurements (m), $z_{\text{oh}}$ is the roughness length governing transfer of heat and vapor (m), $k$ is the von Karman's constant (0.41 (–)) and $u_z$ is the wind speed at height $z_{\mathrm{m}}$ ($\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$). For grass, empirical equations are developed (FAO-56 approach) to estimate $d$, $z_{\text{om}}$ and z${}_{oh}$: $$\begin{array}{}& d=& 0.66V& z_{\text{om}}=& 0.123V& z_{\text{oh}}=& 0.1z_{\text{om}},\end{array}$$ where $V$ is the vegetation height. Wallace et al. (1984) found comparable coefficients for heather: $d=0.63V$ and $z_{\text{om}}=0.13V$ and therefore Eqs. (9)–(11) were applied for both surfaces using a constant vegetation height of 7 and 31 $\mathrm{cm}$ for the grass and heather surfaces, respectively. For the moss surface, we used a vegetation height of 2 $\mathrm{cm}$, which is equal to the thickness of the moss mat. For the bare sand surface we assumed $d=\mathrm{0}$ $\mathrm{m}$, and used typical surface roughness values published by Oke (1978): $z_{\text{oh}}=0.001$ $\mathrm{m}$ and $z_{\text{om}}=z_{\text{oh}}$.

Surface resistance ($r_{\mathrm{s}}$) and canopy interception

Canopy interception was simulated as a water storage which needs to be filled before rainwater reaches the soil surface. A maximum storage capacity of 0.50 $\mathrm{mm}$ was defined for heather following the study of Ladekarl et al. (2005). To our knowledge no literature value of the interception capacity of the specific grass species (Agrostis vinealis) is published. Considering the relatively low vegetation height, we assumed a maximum interception capacity of 0.25 $\mathrm{mm}$.

We distinguished wet ($r_{\text{swet}}$) and dry canopy surface resistance ($r_{\mathrm{s}}$), since interception water evaporates without the interference of leaf stomata. During canopy interception (i.e., if the interception store is fully or partly filled), we used a surface resistance of 0 $\mathrm{s}{\mathrm{m}}^{-\mathrm{1}}$, reducing Eq. (1) to the Penman equation (Penman, 1948; Monteith and Unsworth, 1990). After the canopy storage is emptied, the surface resistance switches to $r_{\mathrm{s}}$. The $r_{\mathrm{s}}$ was back-calculated for daytime periods for the heather and grass lysimeters by substituting measured $R_{\mathrm{n}}$, $G$, ET${}_{\mathrm{a}}$, $e_{\mathrm{s}}$ and $e_{\mathrm{a}}$ and simulated $r_{\mathrm{a}}$ into Eq. (1) under nonstressed conditions (i.e., ET${}_{\mathrm{p}}=$ ET${}_{\mathrm{a}}$). Nighttime evaporation was assumed to be equal to 0 $\mathrm{mm}$. To make sure that the back-calculated $r_{\mathrm{s}}$ was based on days at which evapotranspiration occurred at a potential rate, it was back-calculated for every two consecutive days after precipitation events and after emptying of the (calculated) interception store. The surface resistance ($r_{\mathrm{s}}$) of bare sand and moss was assumed to be equal to 10 $\mathrm{s}{\mathrm{m}}^{-\mathrm{1}}$, i.e., similar to the surface resistance under well-watered conditions of bare soil found by Van de Griend and Owe (1994).

During the summer of 2013, a dry spell (from 4 until 25 July 2013) resulted in a vegetation dieback of grass and heather. Surface resistances were back-calculated for periods before and after the drought event. The drought event had 22 consecutive dry days with a cumulative reference evapotranspiration of 85 $\mathrm{mm}$ according to Makkink (1957). Drought events of similar magnitude have been recorded 12 times during the past 57 years (from 1958 until 2014) at climate station “de Bilt” located in the center of the Netherlands (52.1${}^{\circ }$ latitude, 5.18${}^{\circ }$ longitude), 10 $\mathrm{km}$ from the measurement site. The measurements in the heather vegetation started a week before the drought event. During this week, there were 2 days (30 June and 1 July 2013) for which $r_{\mathrm{s}}$ could be back-calculated. The estimated $r_{\mathrm{s}}$ for these days was 35 and 107 $\mathrm{s}{\mathrm{m}}^{-\mathrm{1}}$ respectively. We selected the $r_{\mathrm{s}}$ value of the second day to use in our model simulations (107 $\mathrm{s}{\mathrm{m}}^{-\mathrm{1}}$) because it was in close agreement with the median surface resistance found by Miranda et al. (1984) of 110 $\mathrm{s}{\mathrm{m}}^{-\mathrm{1}}$ in a comparable heather vegetation. After the drought event, $r_{\mathrm{s}}$ increased to 331 $\mathrm{s}{\mathrm{m}}^{-\mathrm{1}}$ ($N=14$, standard error $=102$ $\mathrm{s}{\mathrm{m}}^{-\mathrm{1}}$). For the grass vegetation, the surface resistance before the drought event was 181 $\mathrm{s}{\mathrm{m}}^{-\mathrm{1}}$ ($N=\mathrm{9}$, standard error $=68$ $\mathrm{s}{\mathrm{m}}^{-\mathrm{1}}$). After the drought event, the surface resistance increased to 351 $\mathrm{s}{\mathrm{m}}^{-\mathrm{1}}$ ($N=\mathrm{4}$, standard error $=47$ $\mathrm{s}{\mathrm{m}}^{-\mathrm{1}}$). Since mosses of these habitats are desiccation-tolerant and quickly rehydrate after drought (Proctor et al., 2007), we did not assess the effect of the dry spell on the surface resistance of the moss surface.

The parameters thus obtained were used to parameterize the Penman–Monteith equation and to calculate hourly ET${}_{\mathrm{p}}$ values for each surface.

Model simulations of ET${}_{\mathrm{a}}$

Using hourly ET${}_{\mathrm{p}}$ for the year 2013 (876 mm precipitation), we used Hydrus 1D (Šimůnek et al., 2008) to simulate ET${}_{\mathrm{a}}$. If meteorological data of the local weather station were missing due to snow cover or sensor maintenance, the meteorological data of weather station “de Bilt” were used for the calculation of ET${}_{\mathrm{p}}$.

First, we simulated ET${}_{\mathrm{a}}$ for the lysimeter surfaces and compared our results with the lysimeter measurements of ET${}_{\mathrm{a}}$. The lower boundary condition in the model was a seepage face with hydraulic pressure equal to 0 at a depth of 65 $\mathrm{cm}$ below the surface (50 cm soil and 15 cm wick). This boundary condition assumes that the boundary flux will remain zero as long as the pressure head is negative. When the lower end of the soil profile becomes saturated, a zero pressure head is imposed at the lower boundary and outflow calculated accordingly. Second, we simulated ET${}_{\mathrm{a}}$ for the groundwater-independent surroundings. We expected that the availability of soil moisture in the lysimeter tanks to be larger than in the groundwater-independent surroundings, because the lowest sections of the lysimeters need to be saturated before drainage occurs. To estimate the yearly ET${}_{\mathrm{a}}$ of dune vegetation in environments with deep groundwater levels, we used a free drainage boundary condition (i.e., a pressure head gradient of 0 and an elevation head of 1) located 2.5 $\mathrm{m}$ below the surface. Third, we investigated the magnitude of the vegetation dieback in the summer of 2013 on both ET${}_{\mathrm{p}}$ and ET${}_{\mathrm{a}}$, by using two different surface resistances: one derived from the period before, and one for the period after the vegetation dieback.

Soil hydraulic properties in the hydrological model were described by the Van Genuchten relationships (Van Genuchten, 1980). Soil samples (100 ${\mathrm{cm}}^{\mathrm{3}}$) collected next to each lysimeter at 5 and 15 cm depth were used to derive the drying retention function. The average drying retention parameters (of the two samples collected next to each lysimeter) were used in the hydrological model, taking hysteresis into account by assuming the wetting retention curve parameter (${\mathit{\alpha }}_{\text{wet}}$) to be twice as large as the drying retention curve parameter (${\mathit{\alpha }}_{\text{dry}}$) (Šimůnek et al., 1999). The unsaturated hydraulic properties (parameters $l$ and $K_{\mathrm{0}}$) were estimated using the Rosetta database and pedotransfer functions, providing the fitted drying retention curve parameters as input (Schaap et al., 2001). The hydraulic properties of the 15 cm long wick, guiding drainage water below the lysimeter into the tipping bucket, were taken from Knutson and Selker (1994).

Since mosses have neither leaf stomata nor roots, ET${}_{\mathrm{a}}$ from the moss surface is limited by the capacity of the moss material to conduct water to the surface. This passive evaporation process is similar to the process of soil evaporation, i.e., evaporation becomes limited if the surface becomes too dry to deliver the potential rate. The unsaturated hydraulic properties of the dense Campylopus introflexes moss mat covering the lysimeter soil were based on the hydraulic properties derived by Voortman et al. (2014) and used in the first 2 cm of the model domain. Macro pores in the moss mat were neglected by Voortman et al. (2014), which implies that direct implementation of these hydraulic properties would result in large amounts of surface runoff generation or ponding, since the unsaturated hydraulic conductivity ($K_{\mathrm{0}}$) of the moss mat is lower than 0.28 $\mathrm{cm}{\mathrm{d}}^{-\mathrm{1}}$. Therefore, the dual porosity model of Durner (1994) was used to add 1000 $\mathrm{cm}{\mathrm{d}}^{-\mathrm{1}}$ to the hydraulic conductivity curve of Voortman et al. (2014) between $-$1 and 0 cm pressure head (Appendix A). This permits the infiltration of rainwater at high intensity rain showers without affecting the unsaturated hydraulic behavior at negative pressure heads. Because of the complex shape of the retention function of the moss mat, hysteresis in the soil hydraulic functions in the underlying soil was neglected for the simulation of evaporation from moss surfaces. The sensitivity of this simplification on the model outcomes was investigated by adjusting the soil hydraulic function of the soil from the drying to the wetting curve. This had a negligible effect ($<$ 1 $\mathrm{mm}$) on the simulated yearly ET${}_{\mathrm{a}}$ (data not shown). Besides simulations of moss evaporation with a cover of Campylopus introflexus, soil physical characteristics of Hypnum cupressiforme were used in the first 2 $\mathrm{cm}$ of the model domain to analyze the effect of different moss species on the water balance. Soil parameters used in the model are explained in more detail in Appendix A.

Since the grass and heather lysimeters fully covered the soil, soil evaporation was neglected for these surfaces. The root profile for the grass and heather lysimeters was 30 cm deep, with the highest concentration of roots in the upper layer decreasing linearly with depth. A water stress reduction function (Feddes et al., 1978) was used to simulate the closure of leaf stomata during water-stressed periods. Vegetation parameters are explained in more detail in Appendix B. Modeled actual evapotranspiration (ET${}_{\text{a,mod}}$) was aggregated to daily values and compared to field measurements of ET${}_{\mathrm{a}}$ during moist (ET${}_{\text{a,mod}}=$ ET${}_{\mathrm{p}})$ and dry conditions (ET${}_{\text{a,mod}}\ne$ ET${}_{\mathrm{p}}$).

Model performance assessment

Model performance of $R_{\text{ns}}$, $R_{\text{nl}}$, $G$ and ET${}_{\text{a,mod}}$ simulations were tested with the Nash–Sutcliffe model efficiency coefficient (NSE): $$\text{NSE}=\mathrm{1}-\frac{{\sum }_{t=\mathrm{1}}^N{\left(x_{\mathrm{o},t},-,x_{\mathrm{m},t}\right)}^{\mathrm{2}}}{{\sum }_{t=\mathrm{1}}^N{\left(x_{\mathrm{o},t},-,{\stackrel{\mathrm{‾}}{x}}_{\mathrm{o}}\right)}^{\mathrm{2}}}.$$ where $N$ is the total number of observations, $x_{\mathrm{m},t}$ is the model-simulated value at time step $t$, $x_{\mathrm{o},t}$ is the observed value at time step $t$, and $\overline{x}$ is the mean of the observations. NSE $=$ 1 corresponds to a perfect match of modeled to observed data. If NSE $<$ 0, the observed mean is a better predictor than the model. To assess the magnitude of error of model simulations, the root mean square error (RMSE), the mean difference (MD) and the mean percentage difference (M%D) were used.

Results and discussion

Parameterization of the Penman–Monteith equation

Net shortwave radiation

The measured incoming and reflected solar radiation were used to compute the albedo of the four surfaces by linear regression (Fig. 6; Table 5). This single value for the albedo slightly overestimates the reflected solar radiation at large incoming solar radiation (Fig. 7) because of a dependency of the albedo on solar elevation angle $\mathit{\beta }$ (Yang et al., 2008; Zhang et al., 2013). Nonetheless the use of a single value for the albedo hardly affects the error in modeled $R_{\text{ns}}$; the mean difference (MD) between measured and modeled $R_{\text{ns}}$ lies between $-$0.23 and 1.63 ${\mathrm{Wm}}^{-\mathrm{2}}$ (Table 1), which is equal to the energy required to evaporate 0.008 to 0.057 $\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$. The NSE for estimating $R_{\text{ns}}$ is close to 1 (Table 1), showing almost a perfect match of modeled to observed data.

Linear regression between incoming and reflected solar radiation.

[Figure omitted. See PDF]

The dense moss mat Campylopus introflexes entirely covers the underlying mineral soil, which results in a low albedo (0.135) due to the dark green surface. The albedo of bare sand (0.261) is comparable to values found in literature for bare dry coarse soils (Qiu et al., 1998; Van Bavel and Hillel, 1976; Linacre, 1969; Liakatas et al., 1986) and the albedo for grass (0.179) is consistent with values reported in other studies during summer time (Hollinger et al., 2010) or for dried grass (Van Wijk and Scholte Ubing, 1963). Heather has a somewhat lower albedo (0.078) than was found in the literature: Miranda et al. (1984) report an albedo of 0.13 (Calluna, LAI ca. 4); Wouters et al. (1980) report an albedo of 0.102 (Calluna). The heather vegetation in our study was in a later successional stage with aging shrubs having a relatively large fraction of twigs and a smaller LAI (3.47) than found by Miranda et al. (1984). Furthermore, the albedo data of heather vegetation were collected primarily past the growing season from September till November. The darker surface after the growing season and the lower LAI explains the small albedo compared to other studies.

Modeled compared to measured net solar radiation (panels a–d, dashed lines are $\mathrm{1}:\mathrm{1}$ lines) and deviations from the $\mathrm{1}:\mathrm{1}$ line (panels e–h; dashed lines indicate 5, 50 and 95 percentiles).

[Figure omitted. See PDF]

Net longwave radiation

The fitted function of Eq. (6) describes the dynamics of the surface temperature relative to air temperature (Fig. 8, Table 5). All surfaces have a similar average nighttime surface temperature ($T_{\text{s,offset}}$) relative to $T_{\mathrm{a}}$, ranging between $-$7.47 and $-$10.21 ${}^{\circ }$C. The solar elevation angle at which the surfaces become warmer than the air temperature (${\mathit{\mu }}_{\mathit{\beta }}$), as well as the speed at which the surface warms up or cools down (${\mathit{\sigma }}_{\mathit{\beta }}$), are comparable between the surfaces. The main difference between the surfaces is observed at high solar elevation angles. Sand and moss show a clear increasing slope during the day, while grass and heather are able to attenuate the increase in surface temperature, possibly due to a larger latent heat flux (Fig. 8). The moss surface shows the largest increase in surface temperature during the day. Although organic layers, e.g., dry peat, have a larger specific heat (1600 $\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$) than dry sand (693 $\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$) (Gavriliev, 2004), the energy required to heat up the moss material is much smaller than for sand, because of the small dry bulk density of ca. 26.8 $\mathrm{g}{\mathrm{L}}^{-\mathrm{1}}$ (derived for Campylopus introflexus from Voortman et al., 2014). Therefore, the surface temperature and the emitted longwave radiation are largest for the moss surface.

Model performance of $R_{\text{ns}}$ simulations.

Surface	$N$	NSE	RMSE	MD	M%D
			($\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$)	($\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$)	(%)
Sand	218	0.998	5.99	$-$0.23	$-$0.10
Moss	1317	0.999	5.46	1.18	0.46
Grass	1203	0.998	7.78	1.63	0.55
Heather	407	0.999	3.00	0.24	0.09

Calibrated net emissivity parameters of the FAO-56 $R_{\text{nl}}$ submodel (Eq. 7).

	$a$	$b$
Sand	0.31	$-$0.00
Moss	0.33	0.02
Grass	0.36	$-$0.06
Heather	0.24	0.02

Model performance of $R_{\text{nl}}$ simulations for hourly time steps.

Surface	$N$	NSE	RMSE	MD	M%D
			($\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$)	($\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$)	(%)
Using Eq. (3)
Sand	5891	0.65	27.37	0.92	1.52
Moss	5997	0.74	28.57	3.73	5.19
Grass	6113	0.71	25.66	1.41	2.36
Heather	2424	0.63	27.63	$-$0.21	$-$0.40
Using FAO-56 Eq. (7)
Sand	5891	0.41	35.39	4.34	7.14
Moss	5997	0.31	46.67	14.84	20.65
Grass	6113	$-$0.07	49.41	$-$18.23	$-$30.38
Heather	2424	0.29	38.24	10.50	19.54

Measured surface temperature relative to air temperature ($T_{\mathrm{s}}$–$T_{\mathrm{a}}$) for clear hours ($f_{\text{cd}}>$ 0.9) as a function of solar elevation angle $\mathit{\beta }$. Relationships (red lines) were fitted to the data using Eq. (6).

[Figure omitted. See PDF]

Our $R_{\text{nl}}$ model (Eqs. 3 and 6) simulates $R_{\text{nl}}$ much better than the calibrated (Table 2) FAO-56 $R_{\text{nl}}$ submodel (Table 3). For the natural grass surface, the NSE even becomes negative using the calibrated FAO-56 approach. Several studies showed that the FAO-56 $R_{\text{nl}}$ submodel underestimates the magnitude of $R_{\text{nl}}$ for reference grass vegetation and poorly describes the diurnal pattern (Matsui, 2010; Blonquist Jr. et al., 2010; Yin et al., 2008; Temesgen et al., 2007). As mentioned, the FAO-56 $R_{\text{nl}}$ submodel was originally developed for reference grass vegetation under well-watered conditions for daily time steps, but is commonly applied at hourly timescales (ASCE-EWRI, 2005; Perera et al., 2015; Gavilán et al., 2008; López-Urrea et al., 2006; Irmak et al., 2005). At daily time steps, $T_{\mathrm{s}}$ is close to $T_{\mathrm{a}}$, since the warmer daytime $T_{\mathrm{s}}$ is compensated by the cooler nighttime $T_{\mathrm{s}}$. For hourly time steps, the assumption that $T_{\mathrm{s}}$ follows $T_{\mathrm{a}}$ is not valid, which explains the poor performance of the FAO-56 $R_{\text{nl}}$ model for hourly time steps. This poor performance cannot be compensated by calibrating the net emissivity parameters, since the diurnal pattern remains unaffected.

In this analysis a typical pattern in $T_{\mathrm{s}}$ relative to $T_{\mathrm{a}}$ is used to estimate $T_{\mathrm{s}}$ (Eq. 6), and subsequently $R_{\text{nl}}$ (Eq. 3). This relationship (Fig. 8) is sensitive to local weather conditions, which implies that the parameters of Eq. (6) (Table 5) are not directly transferable to other locations or climates. The applicability of the presented approach to simulate $R_{\text{nl}}$ should be tested before it is used for other surfaces or climates. It should be noted that the number of parameters that are required to simulate $R_{\text{nl}}$ is relatively large. However, ${\mathit{\mu }}_{\mathit{\beta }}$ as well as ${\mathit{\sigma }}_{\mathit{\beta }}$, are comparable between the surfaces. These parameters might be assumed similar for every surface, reducing the species-specific model parameters to three (one more than the FAO-56 approach). More data of different vegetation types are required to generalize these results and to assess the number of parameters that are required to accurately simulate $R_{\text{nl}}$.

Soil heat flux

The soil heat flux $G$ as a fraction of $R_{\mathrm{n}}$ ($F_{\text{day}}$ and $F_{\text{night}}$) decreases with vegetation cover (Table 5). The nighttime fractions are larger than the daytime fractions, as $R_{\mathrm{n}}$ becomes smaller in magnitude during the night, which simultaneously corresponds to a change in direction of $R_{\mathrm{n}}$ and $G$, from downward (positive) to upward (negative). Relatively small systematic errors are made using daytime and nighttime fractions of $R_{\mathrm{n}}$ to simulate $G$ (MD between 1.92 and 0.69 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$) (Table 4). In remote sensing algorithms $G$ is often simulated as fraction of $R_{\mathrm{n}}$, depending on the LAI or the fractional vegetation cover. In e.g., the SEBS algorithm, the soil heat flux fraction ($F$) is interpolated between 0.35 for bare soil and 0.05 for a full vegetation canopy (Su, 2002). These limits are close to the bare sand (0.270) and heather (0.066) $F_{\text{day}}$ fractions (Table 5). The heather $F_{\text{day}}$ (0.066) was close to the value found by Miranda et al. (1984) of 0.04.

The analysis of the relationship between $R_{\mathrm{n}}$ and $G$ was based on the average of two sets of soil heat flux plates per surface. These sets of measurements showed on average a good agreement: a MD below 1.07 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$, with a RMSE ranging between 5.02 and 9.40 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$.

Model performance of $G$ simulations.

Surface	$N$	NSE	RMSE	MD	M%D
			($\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$)	($\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$)	(%)
Sand	6080	0.820	20.06	1.92	22.16
Moss	5335	0.901	12.02	1.65	24.29
Grass	6046	0.868	8.97	1.60	43.42
Heather	2028	0.641	11.39	0.69	40.27

Parameters of the four different surfaces used for the calculation of ET${}_{\mathrm{p}}$ for hourly time steps.

Parameter	Sand	Moss	Grass	Heather
Albedo (–)	0.261	0.135	0.179	0.078
${\mathit{\mu }}_{\mathit{\beta }}$ (radians)	0.10	0.10	0.13	0.09
${\mathit{\sigma }}_{\mathit{\beta }}$ (radians)	0.09	0.09	0.11	0.08
$T_{\text{s,amp}}$ (${}^{\circ }$C)	11.26	14.21	19.70	15.89
$T_{\text{s,offset}}$ (${}^{\circ }$C)	$-$7.47	$-$8.14	$-$10.21	$-$9.67
$T_{\text{s,slope}}$ (${}^{\circ }\mathrm{C}{\mathrm{radians}}^{-\mathrm{1}}$)	7.83	11.82	0.00	0.00
$F_{\text{day}}$ (–)	0.270	0.211	0.129	0.066
$F_{\text{night}}$ (–)	0.761	0.647	0.527	0.462
$r_{\text{swet}}$ ($\mathrm{s}{\mathrm{m}}^{-\mathrm{1}}$)	–	–	0	0
$r_{\mathrm{s}}$ ($\mathrm{s}{\mathrm{m}}^{-\mathrm{1}}$) before drought	10	10	181	107
$r_{\mathrm{s}}$ ($\mathrm{s}{\mathrm{m}}^{-\mathrm{1}}$) after drought	10	10	351	331

Energy balance

All the terms in the energy balance can be defined using daily lysimeter measurements of LE (latent heat flux) and an estimate of the sensible heat flux ($H$) as a residual term of the energy balance. For daytime measurements (between sunrise and sunset), the LE, $H$, $G$, $R_{\mathrm{s}\uparrow }$ and $R_{\text{nl}}$ can be expressed as fraction of the $R_{\mathrm{s}\downarrow }$. Table 6 summarizes the average fraction of $R_{\mathrm{s}\downarrow }$ attributed to these five different energy fluxes during the measurement campaign. The net longwave radiation is for most surfaces the largest energy flux during daytime (Table 6).

The LE of most surfaces is the second largest flux during daytime, of which fraction increases with vegetation cover. Despite the large difference in albedo between bare sand and moss, the moss surface has only a slightly larger LE fraction than bare sand (Table 6). This is primarily caused by the larger $R_{\text{nl}}$ flux of moss, which compensates the smaller amount of reflected solar radiation.

Potential and actual evapotranspiration

The modeled ET${}_{\mathrm{a}}$ is in agreement with the measured ET${}_{\mathrm{a}}$, with some exceptions at the onset of dry out events (Fig. 9). In general, the reduction of ET${}_{\mathrm{p}}$ to ET${}_{\mathrm{a}}$ is modeled a few days later than it emerges from measurements. The cumulative ET${}_{\text{a,mod}}$ over the measurement period (May–October 2013) deviates 21 $\mathrm{mm}$ (13 %), $-$13 $\mathrm{mm}$ ($-$7 %), 5 $\mathrm{mm}$ (2 %) and $-$3 $\mathrm{mm}$ ($-$2 %) from the measured ET${}_{\mathrm{a}}$ of the sand, moss, grass and heather lysimeters, respectively. The results of modeled vs. measured ET${}_{\mathrm{a}}$ for non-water-stressed (ET${}_{\mathrm{a}}=$ ET${}_{\mathrm{p}}$) and water-stressed conditions (ET${}_{\text{a,mod}}<$ ET${}_{\mathrm{p}}$) are summarized in Table 7.

We did not calibrate our model, e.g., by adjusting soil hydraulic properties, because several processes outlined by Allen et al. (1991) and wall flow (Cameron et al., 1992; Corwin, 2000; Till and McCabe, 1976; Saffigna et al., 1977) affect lysimeter measurements of ET${}_{\mathrm{a}}$ and drainage. We suspect that wall flow caused the slightly earlier reduction of ET${}_{\mathrm{p}}$ to ET${}_{\mathrm{a}}$ at the onset of dry out events than was simulated by the model. Wall flow leads to a quicker exfiltration of rainwater and a subsequent lower moisture content in the lysimeter, and therefore a slightly earlier timing of drought compared to the model. Since wall flow does not occur in the undisturbed vegetation outside the lysimeters, calibrating e.g., soil hydraulic properties using measured surface and drainage fluxes in the objective function could lead to biased characterizations of the soil hydraulic properties and erroneous simulations of soil water flow and ET${}_{\mathrm{a}}$.

Average fractionation of the incoming shortwave radiation ($R_{\mathrm{s}\downarrow }$) between different energy fluxes during daytime.

Surface	LE	$H$	$G$	$R_{\mathrm{s}\uparrow }$	$R_{\text{nl}}$
Sand	0.22	0.13	0.10	0.26	0.28
Moss	0.24	0.17	0.09	0.14	0.36
Grass	0.27	0.21	0.06	0.18	0.29
Heather	0.35	0.20	0.05	0.08	0.32

Measured and modeled daily ET for the four lysimeters. Gray bars indicate time periods where ET${}_{\text{a,mod}}$ is smaller than ET${}_{\mathrm{p}}$, i.e., when evapotranspiration was water-limited.

[Figure omitted. See PDF]

In our simulations, we neglected vapor flow within the soil and moss layer. Due to temperature and potential gradients, vapor fluxes may occur through the soil and moss layer in upward and downward direction by diffusion. Vapor flow may occur by advection as well, e.g., through macropores. Water and vapor flows act together and are hard to distinguish. Modeling and lab experiments show a minor cumulative effect of vapor flow on evaporation for moist and temperate climates. Soil evaporation in a temperate climate for loamy sand in Denmark was only slightly smaller (1.5 %) than a simulation excluding vapor flow (Schelde et al., 1998). Experiments of Price et al. (2009) show that only 1 % of the total water flux was caused by vapor flow in columns of Sphagnum moss. Nevertheless, for a dry and warm Mediterranean climate – different from ours – Boulet et al. (1997) found a dominant vapor flux down to a depth of 25 $\mathrm{cm}$ in a bare soil during 11 days in a dry and warm Mediterranean climate. Because large temperature and potential gradients occur when ET${}_{\mathrm{a}}\ne$ ET${}_{\mathrm{p}}$, vapor flow could especially become dominant in the water-limited phase of evaporation. We compared the model performance between dry (ET${}_{\text{a,mod}}\ne$ ET${}_{\mathrm{p}}$) and wet (ET${}_{\text{a,mod}}=$ ET${}_{\mathrm{p}}$) days in Fig. 10. The model performance in both moisture conditions is comparable (RMSE of sand when dry was 0.40, when wet 0.46; RMSE of moss when dry was 0.30, when wet 0.39), suggesting that our simplified model could describe the dominant processes and the simulation of vapor flow was not required for the temperate climate of our study area.

Modeled ET${}_{\mathrm{p}}$ and ET${}_{\mathrm{a}}$ for different surfaces in a lysimeter (lys.) and for a situation with deep groundwater levels (gw. ind.) for the year 2013.

	ET${}_{\mathrm{p}}$	ET${}_{\mathrm{a}}$ lys.	ET${}_{\mathrm{a}}$ gw. ind.
	(mm)	(mm)	(mm)
Bare sand	400	295	258
Moss (Campylopus int.)	468	312	272
Moss (Hypnum cup.)	468	–	182
Grass	392	350	333
Grass, no dieback	429 ($+$9 %)	382 ($+$9 %)	362 ($+$9 %)
Heather	549	460	391
Heather, no dieback	610 ($+$11 %)	499 ($+$8 %)	420 ($+$7 %)

Measured vs. modeled ET${}_{\mathrm{a}}$ of the lysimeters for all, wet (ET${}_{\text{a,mod}}=$ ET${}_{\mathrm{p}}$) and dry (ET${}_{\text{a,mod}}\ne$ET${}_{\mathrm{p}}$) days. Dotted lines represent the $\mathrm{1}:\mathrm{1}$ lines.

[Figure omitted. See PDF]

One would expect oasis effects to occur in the vicinity of the lysimeters, because freely draining lysimeters must saturate at the bottom of the lysimeter tank before water drains out. This enlarges the water availability inside the lysimeters compared to its groundwater-independent surroundings and occasionally leads to a situation in which the vegetation inside the lysimeters is still transpiring, while the vegetation outside the lysimeters becomes water-stressed and heats up. In such a situation, advection of sensible heat generated in the vicinity of the lysimeters could contribute to the available energy for lysimeter evapotranspiration. However, calculated ET${}_{\mathrm{p}}$ was seldom smaller than measured lysimeter ET${}_{\mathrm{a}}$, indicating that oasis effects were absent. Furthermore, if oasis effects were prominent, systematic underestimation of modeled lysimeter ET${}_{\mathrm{a}}$ would occur, since we ignored the possible contribution of heat advection. Note that it is very unlikely that oasis effects affected the back-calculated surface resistances (Table 5), since these were based on days after rain events for which we may assume ET${}_{\mathrm{a}}$ to be equal to ET${}_{\mathrm{p}}$ for both the lysimeters and their surroundings.

Neglecting feedbacks of drought on the transpiring leaf area and thereby the surface resistance (i.e., using a fixed $r_{\mathrm{s}}$) of heather and dry grassland vegetation leads to an overestimation of cumulative ET${}_{\mathrm{a}}$ of 7–9 % for years with relatively severe drought (Table 7). The delayed drought response of these vegetation types is therefore of importance to water balance studies, especially when, according to the expectations, summers become dryer as a result of a changing climate. Longer recordings of ET${}_{\mathrm{a}}$ in heathland and grassland are required to understand and parameterize the drought response of these vegetation types in coupled plant physiological and hydrometeorological models.

To our knowledge, this paper describes for the first time the evaporation characteristics of a moss surface in a dune ecosystem in a temperate climate. The evaporation rate of the dense moss mat Campylopus introflexus is 5 % larger than the evaporation rate of bare sand. Campylopus introflexus forms dense moss mats and of the moss species investigated by Voortman et al. (2014), it has the largest water-holding capacity. Voortman et al. (2014) hypothesized that moss-covered soils could be more economical with water than bare soils, since the unsaturated hydraulic properties of moss layers reduce the magnitude of evaporation under relatively moist conditions. Our simulations of evaporation from the more open-structured Hypnum cupressiforme moss species (common in coastal dunes), which primarily differs in moisture content near saturation compared to Campylopus introflexus (0.20 instead of 0.61), confirms this hypothesis. The simulated evaporation rate for this species was 29 % lower than the evaporation rate of bare soil. From both our measurements and model simulations, xerophytic (drought-tolerant) mosses appear to be very economical with water; their evaporation rate is comparable with that of bare sand, or lower.

Campylopus introflexus is considered an invasive species in the Northern Hemisphere and was first discovered in Europe in 1941 (Klinck, 2010). Considering the large difference in yearly evaporation between Hypnum cupressiforme and Campylopus introflexus species (90 $\mathrm{mm}$), the invasion of the Campylopus introflexus could have had negative impacts on water resources in specific areas which were previously dominated by more open-structured moss species with poorer water retention characteristics. For sustainable management of groundwater resources in coastal and inland sand dunes, an accurate estimate of the groundwater recharge is required. For consultancy about the availability of water, moss species cannot be categorized in a singular plant functional type, since the modulating effect of the moss cover is species-specific. However, in terms of water retention characteristics, the species investigated by Voortman et al. (2014) are distinguished from each other by the water-holding capacity near saturation (${\mathit{\theta }}_{\mathrm{0}}$, Appendix A), which is easily measured in a laboratory. Moss species could be categorized by this characteristic.

Mosses and lichens are common in early successional stages after colonizing and stabilizing drift sand or as understory vegetation in heathlands or grasslands. Vascular plants might benefit from the presence of certain moss species as more water may be conserved in the root zone. On the other hand, field observations show that moss- and lichen-rich vegetation can persist for many decades (Daniëls et al., 2008). Detailed measurements of understory evaporation in heathlands and grasslands are required to unravel the ecological interactions between mosses and vascular plants.

Conclusions

In this study, the net longwave radiation ($R_{\text{nl}}$) appeared to be one of the largest energy fluxes in dune vegetation. The poor performance of the calibrated FAO-56 approach for simulating $R_{\text{nl}}$ for hourly time steps illustrates that this energy flux has attracted insufficient attention in evapotranspiration research. The novel approach presented in this study to simulate $R_{\text{nl}}$ outperformed the calibrated FAO-56 approach and forms an accurate alternative for estimating $R_{\text{nl}}$.

A relatively simple hydrological model could be used to simulate evapotranspiration of dry dune vegetation with satisfactory results. Improvements in terms of climate robustness would be especially achieved if plant physiological processes were integrated in the hydrometeorological model. Without considering the effects of dry spells on the surface resistance ($r_{\mathrm{s}}$) of grassland and heathland vegetation, ET${}_{\mathrm{a}}$ would be overestimated with 9 and 7 % for years with relatively severe drought (drought events with a reoccurrence of once per 5 years ).

Moss species are very economical with water. The evaporation of moss surfaces is comparable or even lower than bare sand. By promoting moss-dominated ecosystems in coastal and inland dunes, the evapotranspiration could be reduced considerably, to the benefit of the groundwater system. Differences in evaporation between moss species are large and should be considered in water balance studies.

Long-term measurements of ET${}_{\mathrm{a}}$ in heathland and grassland are required to study feedbacks between climate and plant physiological processes in order to integrate the drought response of natural vegetation in coupled plant physiological and hydrometeorological models. To understand the ecological interaction between mosses and vascular plants, detailed measurements of understory evaporation in heathlands and grasslands are required.

Soil hydraulic properties for the simulation of unsaturated flow with Hydrus-1D

Unsaturated flow in Hydrus 1D is described by a modified form of Richards' equation: $$\frac{\partial \mathit{\theta }}{\partial t}=\frac{\partial }{\partial z}\left(K,\frac{\partial h}{\partial z},+,K\right),$$ where $K$ is the unsaturated conductivity ($\mathrm{L}{\mathrm{T}}^{-\mathrm{1}}$), $z$ is the vertical coordinate (L) and $t$ is the time (T). The soil hydraulic properties were assumed to be described by the Mualum van Genuchten functions: $$\begin{array}{}& \mathit{\theta }\left(h\right)=& {\mathit{\theta }}_{\mathrm{r}}+\frac{{\mathit{\theta }}_{\mathrm{0}}-{\mathit{\theta }}_{\mathrm{r}}}{{\left[\mathrm{1},+,\mathrm{|},\mathit{\alpha },h,{\mathrm{|}}^n\right]}^m}& K\left(\mathit{\theta }\right)=& K_{\mathrm{0}}{S}_{\mathrm{e}}^l{\left[\mathrm{1},-,{\left(\mathrm{1},-,S_{\mathrm{e}}^{\mathrm{1}/m}\right)}^m\right]}^{\mathrm{2}}\end{array}$$ with $$S_{\mathrm{e}}\left(h\right)=\frac{\mathit{\theta }\left(h\right)-{\mathit{\theta }}_{\mathrm{r}}}{{\mathit{\theta }}_{\mathrm{0}}-{\mathit{\theta }}_{\mathrm{r}}},$$ where $\mathit{\theta }$ is the volumetric water content (${\mathrm{L}}^{\mathrm{3}}{\mathrm{L}}^{-\mathrm{3}}$), $h$ is the soil water pressure head (L), ${\mathit{\theta }}_{\mathrm{0}}$ is an empirical parameter matching measured and modeled $\mathit{\theta }$ (${\mathrm{L}}^{\mathrm{3}}{\mathrm{L}}^{-\mathrm{3}}$), ${\mathit{\theta }}_{\mathrm{r}}$ is the residual water content (${\mathrm{L}}^{\mathrm{3}}{\mathrm{L}}^{-\mathrm{3}}$) and $\mathit{\alpha }$ (${\mathrm{L}}^{-\mathrm{1}}$) and $n$ (–) are empirical shape parameters of the retention function. $K_{\mathrm{0}}$ is an empirical parameter, matching measured and modeled $K$ ($\mathrm{L}{\mathrm{T}}^{-\mathrm{1}}$), $S_{\mathrm{e}}$ is the effective saturation (–), $l$ is the pore-connectivity parameter (–) and $m$ ($=\mathrm{1}-\mathrm{1}/n$) (–) is an empirical parameter. Drying retention data of two soil samples collected next to each lysimeter at 5 and 15 cm depth were used to fit a retention function with the RETC code (Van Genuchten et al., 1991). Hysteresis in the retention function was accounted for by assuming the retention curve parameter $\mathit{\alpha }$ for the wetting curve (${\mathit{\alpha }}_{\text{wet}}$) to be twice as large as $\mathit{\alpha }$ of the drying retention curve (${\mathit{\alpha }}_{\text{dry}}$) (Šimůnek et al., 1999). The unsaturated hydraulic conductivity parameters $l$ and $K_{\mathrm{0}}$ were estimated using the Rosetta database and pedotransfer functions, providing the fitted drying retention curve parameters as input (Schaap et al., 2001). Average parameter values per lysimeter are summarized in Table A1.

The hydraulic properties of the 15 cm long wick, guiding drainage water below the lysimeter into the tipping bucket, were taken from Knutson and Selker (1994) who analyzed the same brand and type of wick, i.e., Peperell $\mathrm{1}/\mathrm{2}$ inch. The $K_{\mathrm{0}}$ of the wick was adjusted to correct for the smaller cross-sectional area of the wick compared to the cross-sectional area of the lysimeter in the 1-D model simulation (Table A1).

The heterogeneous pore structure of the moss material was described by the functions of Durner (1994): $$\begin{array}{cc}& & S_{\mathrm{e}}=w_{\mathrm{1}}{\left(\mathrm{1},+,[,{\mathit{\alpha }}_{\mathrm{1}},h,{]}^{n_{\mathrm{1}}}\right)}^{-m_{\mathrm{1}}}+w_{\mathrm{2}}{\left(\mathrm{1},+,{\left[{\mathit{\alpha }}_{\mathrm{2}},h\right]}^{n_{\mathrm{2}}}\right)}^{-m_{\mathrm{2}}}& K\left(S_{\mathrm{e}}\right)=\\ & & K_{\mathrm{s}}{\frac{{\left(w_{\mathrm{1}},S_{{\mathrm{e}}_{\mathrm{1}}},+,w_{\mathrm{2}},S_{{\mathrm{e}}_{\mathrm{2}}}\right)}^l\left(w_{\mathrm{1}},{\mathit{\alpha }}_{\mathrm{1}},\left[\mathrm{1},-,{\left(\mathrm{1},-,S_{{\mathrm{e}}_{\mathrm{1}}}^{\mathrm{1}/m_{\mathrm{1}}}\right)}^{m_{\mathrm{1}}}\right],+,w_{\mathrm{2}},{\mathit{\alpha }}_{\mathrm{2}},\left[\mathrm{1},-,{\left(\mathrm{1},-,S_{{\mathrm{e}}_{\mathrm{2}}}^{\mathrm{1}/m_{\mathrm{2}}}\right)}^{m_{\mathrm{2}}}\right]\right)}{{\left(w_{\mathrm{1}},{\mathit{\alpha }}_{\mathrm{1}},+,w_{\mathrm{2}},{\mathit{\alpha }}_{\mathrm{2}}\right)}^{\mathrm{2}}}}^{\mathrm{2}},\end{array}$$ where $w_{\mathrm{1}}$ and $w_{\mathrm{2}}$ are weighting factors for two distinct pore systems of the moss layer, a capillary pore system (subscript 1) and a macropore system active near saturation ($h>$ $-$1 $\mathrm{cm}$, subscript 2), and $K_{\mathrm{s}}$ is the hydraulic conductivity at saturation. Average hydraulic parameters of the capillary pore system and the volumetric portion of the macropore system of the moss species Campylopus introflexus and Hypnum cupressiforme were taken from Voortman et al. (2014) (illustrated with dotted lines in Figs. A1 and A2). The ${\mathit{\alpha }}_{\mathrm{2}}$ parameter was fitted to the functions of Voortman et al. (2014) using $K_{\mathrm{s}}=1000$ $\mathrm{cm}{\mathrm{d}}^{-\mathrm{1}}$ and $n_{\mathrm{2}}=\mathrm{2}$ by minimizing the RMSE by generalized reduced gradient nonlinear optimization. Hydraulic parameter values are listed in Table A2.

Hydraulic parameter values of lysimeter soils.

	${\mathit{\theta }}_{\mathrm{r}}$	${\mathit{\theta }}_{\mathrm{0}}$	${\mathit{\alpha }}_{\text{dry}}$	${\mathit{\alpha }}_{\text{wet}}$	$n$	$K_{\mathrm{0}}$	$L$
	(–)	(–)	(${\mathrm{cm}}^{-\mathrm{1}}$)	(${\mathrm{cm}}^{-\mathrm{1}}$)	(–)	($\mathrm{cm}{\mathrm{h}}^{-\mathrm{1}}$)	(–)
Bare sand	0.01	0.367	0.023	0.046	2.945	1.042	$-$0.401
Moss	0.01	0.397	0.019	–	2.335	0.734	$-$0.173
Grass	0.01	0.401	0.025	0.050	2.071	1.119	$-$0.278
Heather	0.01	0.392	0.018	0.036	2.581	0.679	$-$0.186
Wick	0.00	0.630	0.098	0.196	3.610	2.180	0.500

Hydraulic parameter values of the two moss species.

	${\mathit{\theta }}_{\mathrm{r}}$	${\mathit{\theta }}_{\mathrm{s}}$	${\mathit{\alpha }}_{\mathrm{1}}$	$n$	$K_{\mathrm{s}}$	$l$	$w_{\mathrm{2}}$	${\mathit{\alpha }}_{\mathrm{2}}$	$n_{\mathrm{2}}$
	(–)	(–)	(${\mathrm{cm}}^{-\mathrm{1}}$)	(–)	($\mathrm{cm}{\mathrm{h}}^{-\mathrm{1}}$)	(–)	(–)	(${\mathrm{cm}}^{-\mathrm{1}}$)	(–)
Campylopus int.	0.060	0.936	0.080	2.25	41.67	$-$2.69	0.371	45.89	2.00
Hypnum cup.	0.010	0.971	0.013	2.17	41.67	$-$2.37	0.800	16.61	2.00

Water retention functions of two moss species: Campylopus introflexus and Hypnum cupressiforme. The dotted lines indicate the contribution of the capillary pore system, characterized by Voortman et al. (2014).

[Figure omitted. See PDF]

Hydraulic conductivity functions for two moss species: Campylopus introflexus and Hypnum cupressiforme. The dotted lines indicate the contribution of the capillary pore system, characterized by Voortman et al. (2014).

[Figure omitted. See PDF]

Feddes function used in the Hydrus 1-D model to simulate the closure of leaf stomata during water-stressed periods

The Feddes function (Feddes et al., 1978) describes the relative transpiration rate in relation to the soil water pressure head (Fig. B1) (being 0 if transpiration ceases and 1 if it equals potential rate). Near-positive pressure heads, root water uptake ceases due to oxygen stress (P0). At the dry end of the function, root water uptake ceases (P3). The moment at which transpiration becomes limited due to moisture stress is dependent on the potential transpiration rate. At a high potential transpiration rate (5 $\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$ in the model simulation), leaf stomata start to close earlier (P2H) than under a low potential transpiration rate (P2L, 1 $\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$ in the model simulation). Values for the parameters of Fig. B1 are listed in Table B1.

Parameters of the water stress reduction function used in the Hydrus 1-D model.

P0	P1	P2H	P2L	P3	r2H	r2L
$($$\mathrm{cm}$$)$	($\mathrm{cm}$)	($\mathrm{cm}$)	($\mathrm{cm}$)	($\mathrm{cm}$)	($\mathrm{mm}{\mathrm{h}}^{-\mathrm{1}}$)	($\mathrm{mm}{\mathrm{h}}^{-\mathrm{1}}$)
$-$10	$-$25	$-$300	$-$1000	$-$8000	5	1

The relative transpiration rate as a function of soil water pressure head according to Feddes et al. (1978).

[Figure omitted. See PDF]

Acknowledgements

This study was carried out within the joint research programme of the Dutch Water Utility sector (http://www.kwrwater.nl/BTO), and the project Climate Adaptation for Rural Areas (CARE), which was funded by the Knowledge for Climate Program (http://knowledgeforclimate.climateresearchnetherlands.nl/climateadaptationforruralareas). We thank the staff of Vitens for their permission to perform hydrometeorological measurements in one of their drinking water extraction sites. Edited by: I. Neuweiler