The influence of long-term changes in canopy structure on rainfall interception loss: a case study in Speulderbos, the Netherlands

The evaporation of intercepted water by forests is a significant contributor to both the water and energy budget of the Earth. In many studies, a discrepancy in the water and energy budget is found: the energy that is needed for evaporation is larger than the available energy supplied by net radiation. In this study, we analyse the water and energy budget of a mature Douglas fir stand in the Netherlands, for the two growing seasons of 2015 and 2016. Based on the wet-canopy water balance equation for these two growing seasons, derived interception losses were estimated to be 37 and 39 % of gross rainfall, respectively.

We further scrutinized eddy-covariance energy balance data from these two consecutive growing seasons and found the average evaporation rate during wet-canopy conditions was 0.20 mm h⁻¹. The source of energy for this wet-canopy evaporation was net radiation (35 %), a negative sensible heat flux (45 %), and a negative energy storage change (15 %). This confirms that the energy for wet-canopy evaporation is extracted from the atmosphere as well as the biomass.

Moreover, the measured interception loss at the forest was similar to that measured at the same site years before (I = 38 %), when the forest was younger (29 years old, vs. 55 years old in 2015). At that time, the forest was denser and had a higher canopy storage capacity (2.4 mm then vs. 1.90 mm in 2015), but the aerodynamic conductance was lower (0.065 m s⁻¹ then vs. 0.105 m s⁻¹ in 2015), and therefore past evaporation rates were lower than evaporation rates found in the present study (0.077 mm h⁻¹ vs. 0.20 mm h⁻¹ in 2015). Our findings emphasize the importance of quantifying downward sensible heat flux and heat release from canopy biomass in tall forest in order to improve the quantification of evaporative fluxes in wet canopies.

1 Introduction

Rainfall interception is the portion of precipitation temporarily captured by vegetation before evaporating back into the atmosphere. It is by definition unavailable for soil infiltration or runoff. Evaporation from intercepted rainfall is an important component of the water balance, and in coniferous forests it can represent around 25–45 % of gross rainfall (P_G) (Rutter et al., 1975; Gash et al., 1980; Carlyle-Moses and Gash, 2011). Starting from the early twentieth century (Horton, 1919; benchmark papers in Gash and Shuttleworth, 2007) much research has been done to quantify the magnitude of rainfall interception over different forest ecosystems (cf. Carlyle-Moses and Gash, 2011; Muzylo et al., 2009; Miralles et al., 2010; Llorens and Domingo, 2007). Despite decades of research, the physical processes and atmospheric conditions that allow a large fraction of rainfall to return to the atmosphere are still poorly understood (van Dijk et al., 2015). A key issue is that the water budget often suggests a higher evaporation rate than the available energy permits. In addition, little research has focused on the long-term evolution of rainfall interception loss with forest growth and development, and on the implications for the forest water balance. Improving knowledge of the evolution of rainfall interception with forest growth, in particular the canopy storage and wet-canopy evaporation rate, will provide insight into the role of forests regarding moisture recycling, and will assist when developing both forest management strategies against threats to forest water resources and remote sensing techniques to monitor rainfall interception based on forest structure (i.e. Miralles et al., 2010; Cui et al., 2015; Hassan et al., 2017).

In this study, we revisited a site with a Douglas fir plantation in the centre of the Netherlands, the “Speulderbos”, for which historical measurements of rainfall interception are available. During the 1990s research focused on the impact of air pollution on forest growth (“Aciforn project”) at this site (Evers et al., 1991a). Several hydrological studies at that time found high interception losses (∼ 40 %) from the forest (Tiktak and Bouten, 1994; Klaassen et al., 1998; Bouten et al., 1996), contributing to the debate about forest plantations in the centre of the Netherlands affecting the ground water recharge (Moors, 2012). Detailed measurements and modelling studies at the site included canopy growth and architecture (Evers et al., 1991b), soil water dynamics (Tiktak and Bouten, 1994), modelling evaporation and transpiration (Bosveld and Bouten, 2001), spatial patterns of throughfall (Bouten et al., 1992), and measuring and modelling canopy water storage (Bouten et al., 1991, 1996). It was found that the high interception losses in Speulderbos at that time were due to a high leaf area index, resulting in high canopy storage capacity (2.5 mm) as measured by microwave transmission (Bouten et al., 1991). Evaporation during rainfall contributed minimally to the overall interception loss (Klaassen et al., 1998).

The historical data collected in Speulderbos offer a good opportunity to evaluate the effects of long-term changes in canopy structure on the rainfall interception process and their implications for other phases of the water balance (i.e. soil infiltration, plant water uptake, ground water recharge, and stream discharge). Changes in forest structure can be resultant from different factors such as tree phenology, management practices (Bormann et al., 2015), changes in species composition (Thom et al., 2017), and stand development (Franklin et al., 2002; Freund et al., 2015). Although several studies have analysed the effects of short-term changes (i.e. seasonal) in forest canopy structure on rainfall interception (Dolman, 1987; Price and Carlyle-Moses, 2003; Deguchi et al., 2006; Herbst et al., 2008; Muzylo et al., 2012), far fewer studies have considered the effects of forest growth over longer timescales (i.e. decades, centuries) in forest canopy structure on rainfall interception. In some cases comparisons have been limited to monitoring the water balance components of two (or more) stands of different ages concurrently (i.e. Pypker et al., 2005; Keim et al., 2005).

Due to the long-term maintenance of the research infrastructure at the site, a unique opportunity arose to evaluate the changes in eco-hydrological processes over several decades. We were particularly interested in the water and energy budget of the forest. For this, two variables are of key importance: the water storage capacity of the forest, and the evaporation rate of the wet canopy, which is energy limited. The evaporation rate, in turn, depends on the radiation budget, the aerodynamic roughness, and the heat storage capacity. We expected that thinning of the forest stand in combination with a substantial increase in forest height must have affected the water and energy storage capacity, as well as the forest roughness.

For this study, we combined data collected using an eddy-covariance flux tower with precipitation, stemflow, and throughfall data to obtain an estimate of the interception loss from a mature Douglas fir plantation (ca. 55 years old). The objectives of the present study were to

• i.

  assess two indirect methods for estimating canopy storage capacity,

• ii.

  quantify the sources of energy that drive the latent heat flux involved in the evaporation of intercepted rainfall,

• iii.

  examine the effect of long-term changes in canopy structure on the rainfall interception losses, and

• iv.

  explore the relative importance of climatic and forest structural factors to overall rainfall interception loss using a physically based interception model.

2 Materials and methods

2.1 Research site

The study was conducted within a 2.5 ha evergreen Douglas fir (Pseudotsuga menziesii) stand located in the forested area of “Speulderbos” (52^∘15^′04${}^{\prime \prime }$ N, 05^∘41^′25${}^{\prime \prime }$ E) at an elevation of 50 m a.s.l., near the settlement of Garderen, the Netherlands (Fig. 1). The site is equipped with a 47 m scaffolding tower, which supports measurement of a range of micrometeorological data. The plot is surrounded by several stands of other species such as beech, oak, and hemlock. The climate is classified as temperate–humid. Based on “de Bilt” weather station data, located at 38 km south-west of the plot, the average (± SD) annual precipitation for the period 2000–2015 was 864 (±92) mm. In general, July is the wettest month, with about 12 % of the annual rainfall, and April the driest month, with 4 % of the annual rainfall. The mean annual temperature is 10.6 ^∘C (±0.6), with January being the coldest month (3.7 ± 2 ^∘C) and July the warmest month (18.2 ± 1.6 ^∘C) (KNMI, 2015). The type of soil in the study area is Typic Dystrochrepts on thick heterogeneous sandy loam and loamy sand textured ice-pushed river sediments (Tiktak and Bouten, 1988).

Active reforestation in the area, previously sand dunes, started at the end of the nineteenth century. The current stand was planted with 2-year old seedlings in 1962. For the study period canopy height was about 34 m, whereas stem density and mean diameter at breast height (DBH) were 571 trees ha⁻¹ and 34.8 (±8.9) cm, respectively. The leaf area index of the plot (LAI, using a LI-COR LAI 2000 Plant Canopy Analyser) was 4.5 (±0.38) (Fig. 1b). No other tree species were recorded in the plot and understory was largely absent (Fig. 1c).

Figure 1 Study area in Speulderbos: (a) location map in the Netherlands; (b) top view of the canopy; (c) funnel-type collector used to quantify throughfall in the study site.

2.2 Field measurements

2.2.1 Rainfall

Gross rainfall (P_G, mm) was measured in a nearby well-exposed clearing (ca. 250 m from the centre of the plot) using two tipping bucket rain gauges (Rain Collector II, Davis Instruments, USA) with a resolution of 0.2 mm per tip. The orifice of the rain gauge was positioned at 1.5 m above the ground to avoid ground-splash effects. The automatically recorded data were stored by a HOBO event logger at 1 min intervals (Onset Computer Corporation, USA).

Gross rainfall was also collected at the top of the 47 m scaffolding tower operated by the University of Twente (ITC-UT) (at ca. 200 m distance from the clearing), using a tipping bucket rain gauge (Onset HOBO-RG3, resolution 0.2 mm). The data collected at the top of the tower were only used to fill gaps in the data from the clearing (from 23 July to 12 August 2015, as well as 24 May to 9 June 2016) using a linear regression equation linking 10 min rainfall totals of the two locations (R² = 0.93, n = 500).

Table 1 Main micro-meteorological instruments installed on the Speulderbos flux tower.

2.2.2 Throughfall

Throughfall (Tf, mm) was measured by an automated gutter system and validated by an arrangement of manual (roving sampling) funnel-type collectors. The automated gutter system consisted of four stainless steel gutters (200 cm × 30 cm each), randomly located in the plot and connected by pairs to two tipping buckets (V2A UP Umweltanalytische Produkte GmbH). As no apparent alignment of the trees was observed in the planted stand, no specific orientation of the gutters was considered. The gutters were mounted on a wooden frame, about 60 cm from the forest floor and at an inclination of 10 % to facilitate drainage to the tipping buckets. Combining two gutters and correcting for the inclination provided a total catch area of 1.2 m² yielding 0.084 mm per tip. The tipping buckets were connected to a data logger (CR23X, Campbell Scientific Ltd.) and tip pulses were recorded at a 1 min resolution. The gutters and the tubes were cleaned every 7 to 15 days to avoid clogging due to falling litter. In addition, Tf was measured using funnel-type collectors. The manual array of collectors was operated from 17 February to 2 November 2015. A stratified random sampling approach was used to ensure an even spread of sampling locations. We defined a plot size of 32 m × 64 m, which was divided into 32 square sub-plots of 8 m × 8 m each; each sub-plot was marked in its centre. Collectors (32 in total) were placed at some distance from each marked point, by generating random values for an azimuth angle and distance from the grid point. The azimuth angles ranged from 0 to 360^∘ and the distances from 0 to 4 m. In case the randomly selected position coincided with the position of a stem, the azimuth angle was maintained while the distance was adjusted until the collector was located next to the tree base and the adjusted distance was recorded. The funnel-type collectors consisted of a 2 L collector and a funnel (165 cm² orifice area). The orifices of the gauges were positioned 50 cm above the forest floor to avoid splash-in from the ground. The funnel-type collectors were read (and relocated) ∼ bi-weekly (i.e. roving sampling; Ritter and Regalado, 2014). Measured Tf volumes were converted to equivalent depth (in mm) by dividing the volume of water in each gauge by the orifice area.

Using the manual array, Tf was measured 15 times (cf. Table A1). For the first 10 measurement periods, we applied a roving sampling method by randomly relocating the position of the funnel-type collectors after each Tf measurement (resulting in 320 different positions of the funnel-type collectors). However, the coefficient of variation (CV) in the 10 initial measurement periods was low (∼ 15 %) and therefore we did not use the roving technique for the remaining 5 measurement periods (i.e. the gauges were not relocated after measurements were taken).

2.2.3 Stemflow

Following a stratified sampling approach, stemflow (Sf, mm) was measured on four trees with differing DBHs, representative of the whole stand. The four diameter size classes were < 30, 31–40, 41–50, and > 50 cm. Each set-up consisted of a halved plastic tube wrapped around the tree stem in a spiral fashion, starting at a height of ca. 80 cm; the lower end of the tube was connected to a closed tipping bucket (Onset HOBO^® S-RGA Rain Gauge, resolution 0.254 mm). Silicon sealant was applied between the stem and the plastic tube to seal the gaps (and hence avoid stemflow loss). Stemflow proved to be only a minor component of the wet-canopy water balance. As the sampled trees covered the whole range of diameter classes within the plot, total stemflow in the plot was calculated by multiplying the stemflow volumes by the number of trees for each diameter class (cf. Levia and Germer, 2015; Eq. 2). Stemflow measurements were carried out over 113 days from 27 July to 11 November 2016. During this period, a total of 240 mm of rain was received at the plot.

2.2.4 Net radiation and soil heat flux

Net radiation (R_n) was measured by a four-component net radiometer (model CNR1, Kipp and Zonen) mounted at 35 m above ground (Table 1), and averages were stored at 10 min intervals, except during three periods, totalling 120 days, encountering instrument failure (from 24 April to 23 June 2016, from 17 July to 4 August 2016, and from 15 September to 14 October 2016).

Soil heat flux (G) was estimated by combining measurements from two heat flux plates (HFP01, Hukseflux) both installed at 8 cm depth, and temperature profile measurements at five different depths (1, 3, 8, 20, and 50 cm). Estimations of G at ground level were carried out following the harmonics method (Verhoef et al., 1996) with the derivatives of the Fourier series calculated analytically according to van der Tol (2012).

2.2.5 Turbulent heat fluxes

Sensible (H) and latent (λE) heat fluxes were estimated by the eddy-covariance technique with a sonic anemometer (CSAT3, Campbell Sci. Inc.) and an open path gas analyser (LI7500, LI-COR Biosciences) mounted at 47 m (Table 1).

Thirty minute turbulent fluxes were processed with the EddyUH ver. 1.7 software (University of Helsinki, https://www.atm.helsinki.fi/Eddy_Covariance/EddyUHsoftware.php, last access: February 2017). As initial estimates for EddyUH, a displacement height (d) of 0.7 of canopy height (h) (Stull, 2012) and a roughness length (z₀) of 0.06 h (Weligepolage et al., 2012) were used. Furthermore, the following corrections were performed: de-spiking, 2-D coordinate rotation, cross-wind correction to sonic temperature according to Liu et al. (2001), high-frequency spectral corrections according to Moncrieff et al. (1997) and Aubinet et al. (2000), low-frequency spectral corrections according to Rannik and Vesala (1999), correction for humidity effect on sonic heat flux according to Schotanus et al. (1983), and WPL correction according to Webb et al. (1980). We disregarded turbulence data with an overall quality flag above 2 in accordance with the Foken et al. (2005) quality flag system.

2.2.6 Energy storage

Energy storage (Q), composed of energy storage in the canopy air (Q_air) and in the biomass (Q_bio), was estimated based on measurements of a vertical profile of air temperature and humidity. Further details on the sensors and their location are shown in Table 1. Energy storage changes in the air Q_air result from the change in the temperature in the air Q_T plus the component resulting from the specific humidity Q_q (cf. Michiles and Gielow, 2008).

Q_bio was estimated as the sum of energy stored in the trunks (Q_tr), in the branches (Q_br), and in the needles (Q_nd). For Douglas fir in the centre of the Netherlands, allometric equations on biomass and its vertical distribution derived from Bartelink (1996) were used (Table 2). Specific heat of biomass (c_v) was assumed to be equal for all components at 2400 J kg⁻¹ C⁻¹ (Michiles and Gielow, 2008). A moisture content of 44 % for Douglas fir (Nord-Larsen and Nielsen, 2015) was used to estimate the dry matter content.

Table 2 Comparison of stand parameters and biomass dry weight (DW) for the Douglas fir stand in Speulderbos. Above-ground biomass determined by means of stem survey and allometric relationships from Bartelink (1996).

^a Biomass dry weight values were estimated using tree density and DBH from Tiktak and Bouten (1994). ^b LAI measured by using the LI-COR-2000 instrument. ^c Previous studies in Speulderbos report a LAI of 11: that value was estimated by the destructive method (cf. Heij and Schneider, 1991). For comparative reasons, we use the reported value using the LI-COR photometer.

Q_air was estimated by dividing the air column into four sections: Section 1 (0 to 10 m), Section 2 (10 to 20 m), Section 3 (20 to 28 m), and Section 4 (28 to 34 m). Each section was centred on the respective level of temperature and humidity (Table 1) considered representative of the entire section. The following equations suggested by McCaughey (1985) were used to estimate Q_air:

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }{Q}_{\mathrm{air}}=Q_T+Q_q,\text{(2)}& {\displaystyle }& {\displaystyle }{Q}_{\mathrm{air}}=\mathit{\rho }\left(c_{\mathrm{p}}\mathrm{\Delta }\overline{T}+\mathit{\lambda }\mathrm{\Delta }\overline{q}\right)\mathrm{\Delta }z/\mathrm{\Delta }t,\end{array}$$

where ρ (kg m⁻³) is density of air, c_p (J kg⁻¹ K⁻¹) is the specific heat of air, λ (J kg⁻¹) is the latent heat of vaporization, $\mathrm{\Delta }\overline{T}$ (K) and$\mathrm{\Delta }\overline{q}$ (kg kg⁻¹) are the change in mean air temperature and specific humidity over time, respectively, Δz (m) is the height thickness of the considered layer, and Δt (s) is the time interval.

Q_bio was estimated by means of the following equation (Oliphant et al., 2004; McCaughey, 1985):

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }{Q}_{\mathrm{bio}}=m_{\mathrm{bio}}{c}_{\mathrm{v}}\mathrm{\Delta }{T}_{\mathrm{bio}}/\mathrm{\Delta }t,\end{array}$$

where m_bio (kg m⁻²) is the mass of biomass per unit of horizontal area, c_v (J kg⁻¹) is a representative specific heat of the vegetation, and T_bio (K) is a representative biomass temperature. Some studies have used ambient air temperature as a surrogate for T_bio (Oliphant et al., 2004; Thom et al., 1975; Michiles and Gielow, 2008). We used Eq. (3) with T_bio equalling air temperature for intervals without rainfall, and wet bulb temperature (Stull, 2011) for intervals with rainfall (P_G > 0.5 mm) (cf. van Dijk et al., 2015).

2.2.7 Canopy wetness

Three leaf wetness sensors (LWSs) (Model 237, Campbell Sci. Inc.) were installed at 20, 24, and 26 m above ground, in the middle of the living crown. Ringgaard et al. (2014) used 4 LWSs located within the canopy and demonstrated that the parametrization of the analytical model for interception loss by Gash et al. (1995) was improved. An LWS consists of a circuit board that emulates the leaf surface, with interlacing gold-plated fingers on it. The LWS estimates foliage wetness by determining the electrical resistance on the surface of the sensor. Sensors were placed over the needles in the middle of the branches and tilted about 60^∘ to avoid rainwater puddling on the electrodes.

Table 3 Main equations of the analytical Gash (1979) interception model.

2.3 Modelling rainfall interception

2.3.1 The Gash rainfall interception model

The Gash analytical model (Gash, 1979) was used to simulate rainfall interception loss (I, mm). The Gash model assumes that rainfall occurs as a series of discrete events. The Gash model differentiates three phases in a rainfall event: (i) the canopy wetting phase, (ii) the canopy saturation phase, (iii) the canopy drying phase. Table 3 summarizes the equations associated with the respective phases. The Gash model uses four canopy parameters: (i) canopy storage capacity S, which is defined as the amount of water left in a saturated canopy in the absence of evaporation, after rainfall and drainage has ceased, (ii) the free throughfall coefficient p, which is the fraction of incident rainfall that reaches the forest floor without touching the forest canopy, (iii) the coefficient p_t, which is the fraction of rain diverted to the trunks as Sf, (iv) stem storage capacity S_t, which is the amount of water that can be stored on the stems. In addition to the canopy parameters, the Gash analytical model requires two climatic parameters: the mean evaporation rate $\overline{E}$ and the mean rainfall intensity $\overline{R}$.

For modelling purposes, we divided our data set into two parts: data-set 1 included measurements from 19 June to 31 October 2015, and data-set 2 included measurements from 19 June to 31 October 2016. Data-set 1 was used for parameter estimation (model calibration), and data-set 2 was used for validation.

We used two different parametrizations of the Gash model. In the first parametrization (Run 1) the parameters S, p, and $\overline{E}/\overline{R}$ were derived from the water balance (rainfall, throughfall, and stemflow data) by using the mean method (Klaassen, 2001). In the mean method, S, p, and $\overline{E}/\overline{R}$ were derived from linear regressions of measured I vs. measured P_G from multiple events. In the second parametrization (Run 2), the parameters S and p were derived using individual event analysis (IEA) (Link et al., 2004), while the parameter $\overline{E}/\overline{R}$ was calculated with E derived from the energy balance residual and R derived from the tipping bucket measurements. For both parametrizations, values of S_t and p_t as derived with the Gash and Morton (1978) method were used. The methods are discussed in detail in the following sections.

Moreover, the sensitivity of the modelled I to the relevant parameters, namely S and $\overline{E}$ (Gash, 1979; Loustau et al., 1992; Moors, 2012), was evaluated. We used the RMSE as the criterion to test the sensitivity for data-set 1 and data-set 2 separately. By comparing these different runs and the model sensitivity, we were able to evaluate the consistency of the water and energy balance estimates of evaporation as well as the storage capacity of the canopy.

2.3.2 Derivation of canopy parameters

Two methods were used to derive the canopy parameters: the multiple event analysis or mean method (Klaassen et al., 1998) and the individual event analysis (IEA) (Link et al., 2004). As the Gash (1979) model is event based, it is important to discriminate events first. Because the criteria used to discriminate events can have a major effect on the interception modelling, we evaluated two cases of event selection: Case A, considering an event as a period of rain exceeding 0.5 mm preceded by a dry period of at least 3 h (cf. Klaassen et al., 1998), and Case B, considering an event as a period of rain exceeding 0.5 mm, with the preceding dryness validated by three LWSs (all indicating fully dry) within the living crown. In addition, for Case B saturated events were considered only those events with P_G ≥ 5 mm.

The mean method is a multi-event analysis where S and p are estimated by linear regression of interception loss (I = P_G − Tf − Sf) vs. P_G. To estimate p, the regression is in the form of I = aP_G with a = 1 − p− p_t and only events with a total amount of rainfall unable to saturate the canopy (P_G < $P_{\mathrm{G}}^{\prime }$) are used. The parameters S and $\overline{E}/\overline{R}$ are derived from events large enough to saturate the canopy (P_G ≥ $P_{\mathrm{G}}^{\prime }$) with the linear regression in the form of I = b₁P_G + b₂, where b₁ = $\overline{E}/\overline{R}$ and b₂ = (1 − p − $p_{\mathrm{t}})P_{\mathrm{G}}^{\prime }$. An iterative procedure is employed where the initial value of $P_{\mathrm{G}}^{\prime }$ is visually defined from an I vs. P_G graph. After fitting both equations, $P_{\mathrm{G}}^{\prime }$ is re-calculated as $P_{\mathrm{G}}^{\prime }$ = b₂∕(a + b₁) and the process is repeated until $P_{\mathrm{G}}^{\prime }$ converges (cf. Klaassen et al., 1998; Holwerda et al., 2012).

In contrast, the IEA consists of an analysis of the behaviour of water fluxes during individual events. It is based on the equations proposed in the Gash (1979) analytical model. When rainfall starts (P_G < $P_{\mathrm{G}}^{\prime }$) throughfall increases approximately linearly with P_G (Tf = pP_G) until saturation is reached: P_G ≥ $P_{\mathrm{G}}^{\prime }$. When rainfall saturates the canopy storage capacity (P_G = $P_{\mathrm{G}}^{\prime }$) an inflection point in the Tf vs. P_G plot occurs and water starts to drain from the canopy. Canopy parameters p, $P_{\mathrm{G}}^{\prime }$, and S can be derived using an iterative regression procedure over the plot of cumulative Tf vs. cumulative P_G. The procedure was applied to events selected in Case B (i.e. using LWS and saturation threshold P_G ≥ 5 mm) and with data records at a 15 min time resolution. Events that did not have sufficient 15 min records before saturation was reached were discarded. The inflection point ($P_{\mathrm{G}}^{\prime }$) was calculated as the intersection of two linear regressions, for the wetting-up stage and after canopy saturation.

2.3.3 Wet-canopy evaporation

Three methods were used to estimate wet-canopy evaporation rate, based on (i) a water balance approach, (ii) the energy balance residual, and (iii) the Penman–Monteith equation.

Firstly, the wet-canopy evaporation rate ($\overline{E}$) was derived from the value of $\overline{E}/\overline{R}$ (obtained from the mean method). Because the distribution of rainfall intensity was skewed, we used the median rainfall intensity instead of the mean following the recommendations of Schellekens et al. (1999). The thus derived value of $\overline{E}$ will henceforth be referred to as the “water balance based” evaporation rate.

Secondly, wet-canopy evaporation rates were estimated from the energy balance residual. The quality of the energy balance data (eddy-covariance flux of sensible heat, net radiation, ground heat, and storage terms) was verified by calculating the energy balance closure and the energy balance ratio (EBR) for the dry and wet periods for which high-quality data (quality flag ≤ 2) of λE were available. In addition, the performance of the sonic anemometer (CSAT3) during wet conditions was evaluated by plotting the standard deviation of the vertical wind speed (σw) against friction velocity (u^*) (Gash et al., 1999; van der Tol et al., 2003; Holwerda et al., 2012). According to Monin–Obukhov similarity theory $\mathit{\sigma }w/u^{*}$ in neutral conditions is a universal constant; therefore, the ability of the anemometer to measure $\mathit{\sigma }w/u^{*}$ during wet and dry conditions was tested (Gash et al., 1999).

Wet-canopy evaporation rate was derived using the energy balance residual approach where λE (W m⁻²) is estimated as the residual of the energy balance equation as

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }\mathit{\lambda }E=R_{\mathrm{n}}-H-G-Q-G_{\mathrm{P}},\end{array}$$

with H derived from the eddy-covariance technique and G_P, the photosynthetic energy flux, estimated at −2 W m⁻² during the night and at 6 W m⁻² during the day (Thom et al., 1975). Equation (4) provides a more complete data set than λE based on the eddy-covariance data only, due to the fact that the open path gas analyser is prone to providing low-quality data (causing rejections during filtering) during wet periods. The evaporation estimated with Eq. (4) is hereinafter referred to as $E_{\mathrm{EB}-\mathrm{EC}}.$

Finally, the last and most common method to estimate wet-canopy evaporation, the Penman–Monteith (P–M) equation, was used. P–M estimates latent heat flux (λE, W m⁻²) as

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }\mathit{\lambda }E={\displaystyle \frac{\mathrm{\Delta }A+\mathit{\rho }{c}_{\mathrm{p}}\left(e_{\mathrm{s}}-e\right)g_{\mathrm{a}}}{\mathrm{\Delta }+{\mathit{\gamma }}^{\prime }}}\end{array}$$

with

$$\begin{array}{}\text{(6)}& {\displaystyle }{\displaystyle }{\mathit{\gamma }}^{\prime }=\mathit{\gamma }\left(\mathrm{1}+{\displaystyle \frac{g_{\mathrm{a}}}{g_{\mathrm{s}}}}\right),\end{array}$$

where Δ (hPa K⁻¹) is the slope of the saturated water vapour pressure curve, A (W m⁻²) is the available energy, γ^′ and γ (hPa K⁻¹) are the adjusted and original psychrometric constant, respectively, ρ (kg m⁻³) is the density of air, c_p (J kg⁻¹ K⁻¹) is the specific heat of air at constant pressure, e_s (hPa) is the saturation vapour pressure at ambient temperature, e (hPa) is the actual vapour pressure, g_a (m s⁻¹) is the aerodynamic conductance, and g_s (m s⁻¹) is the surface conductance.

During wet conditions, surface conductance is assumed to be infinitely large (i.e. surface resistance set to zero). Aerodynamic conductance for momentum g_a,M (m s⁻¹), following Thom (1975), for neutral conditions was derived from the regression of observed friction velocity (u^*, m s⁻¹) vs. wind speed (u, m s⁻¹) measured by a sonic anemometer (Gash et al., 1999; van der Tol et al., 2003):

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle }{g}_{\mathrm{a},\mathrm{M}}={\left({\displaystyle \frac{u^{*}}{u}}\right)}^{\mathrm{2}}u.\end{array}$$

It is necessary to differentiate between conductance for heat and momentum (Lankreijer et al., 1993). Following the empirical relation proposed by Garratt and Francey (1978), we used ln(z_0M∕z_0H) = 2 for g_a,H (cf. Gash et al., 1999; Moors, 2012; Lankreijer et al., 1993). In addition stability corrections for non-neutral hours were implemented according to Paulson (1970). The evaporation estimated with the P–M equation is hereinafter referred to as $E_{\mathrm{PM}-\mathrm{EC}}.$

3 Results

3.1 Rainfall

Total rainfall (P_G) measured during 11 months between 19 June 2015 and 31 October 2016 (excluding the winter season from November 2015 to March 2016) was 955 mm. Mean monthly precipitation was 82 mm (±41 SD), with a minimum of 43 mm in May 2016 and a maximum of 156 mm in August 2015.

A total of 157 events (64 in 2015 vs. 93 in 2016) were identified from half hour rainfall time series during the study period. The frequency distributions of event size, duration and rainfall intensity showed a positively skewed distribution. The mean (and median) event-based amount, duration and intensity were 6.0 (3.1) mm, 6.5 (5.0) h, and 1.1 (0.74) mm h⁻¹, respectively. The maximal event size, duration and intensity were 66.3 mm, 62 h and 8.5 mm h⁻¹, respectively.

3.2 Throughfall, stemflow, and derived interception loss

Throughfall (Tf), measured for the same period as the rainfall, was 577 mm in total, corresponding to about 60 % of P_G. The overall standard error (SE) of Tf was 48 mm (i.e. 5 % of P_G). No significant differences (t-test; α = 0.05) were found between the mean cumulative Tf estimated using either the manual array or the automated gutters, confirming that the automated system was representative for the plot. The homogeneity in the plot was illustrated by a relatively low coefficient of variation (CV) of ∼ 15 % (cf. Table A1).

The total Sf at plot scale was 2.6 mm (1.1 % of gross rainfall). The contribution to the total Sf from the four stem diametric classes, < 30, 30–40, 40–50, and > 50 cm, was 0.8, 0.2, 0.1, and < 0.1 %, respectively. The SE value of the overall Sf was 2 % based on the four stemflow collectors.

The total interception loss estimated for the whole study period based on the wet-canopy water balance was 372 mm (39 % P_G). The SE for the interception loss estimated as the quadratic mean of the SEs of Tf and Sf was 52 mm (i.e. 5.4 % of P_G).

3.3 Canopy-related parameters

Using the mean method, we analysed data-set 1 (the calibration data set) and found only minor differences in estimated parameter values with respect to the criteria used to select the events (Cases A and B). For Case A, values of p and S were found of 0.32 (±0.04) and 1.15 (±0.25) mm, respectively, while for Case B, values of p and S were found of 0.28 (±0.05) and 1.37 (±0.51) mm, respectively.

When we selected events based on Case A (P_G ≥ 0.5 mm and dryness separation time of 3 h), then the linear regression expression that relates I to P_G was I = 0.65 P_G (R² = 0.68; n = 24) for non-saturated conditions and I = 0.25 P_G + 1.15 (R² = 0.67; n = 40) for saturated conditions in data-set 1 (Fig. 2a). Case B (i.e. excluding events that did not reach the condition of preceding dryness validated with three fully dry LWSs) resulted in linear regressions of I = 0.69 P_G (R² = 0.68; n = 11) for unsaturated conditions (P_G < $P_{\mathrm{G}}^{\prime }$) and I = 0.23 P_G + 1.37 (R² = 0.75; n = 17) for saturated conditions (P_G ≥ 5 mm) (Fig. 2b).

Figure 2 Determination of canopy-related parameters using the mean method and the individual event analysis. (a) Linear regression using data-set 1; events selected in Case A. Circles represent rainfall events with total rainfall less than that necessary for saturation; crosses represent data with enough rainfall to saturate the canopy. (b) Linear regression using data-set 1; events selected in Case B (similar legend to a). (c) Individual event analysis (IEA) on 17 September 2015, the plot of data used to estimate canopy direct throughfall and saturation storage capacity. Dots represent values of cumulative rainfall vs. cumulative Tf. The direct throughfall regression equation was Tf = 0.07 P_G, and the saturation regression equation was Tf = 0.68 P_G − 1.38. Canopy saturation point was calculated as the intersection of the two linear regressions, $P_{\mathrm{G}}^{\prime }$ = 2.29 mm.

The values of the stemflow-related parameters (S_t, p_t) were obtained using the method of Gash and Morton (1978). The trunk storage capacity (S_t) was estimated at 0.14 mm (±0.05) and the proportion of rain diverted to stemflow p_t at 0.029 (±0.005) (R² = 0.77; n = 12).

Using the IEA method, average values of p and S of 0.17 (±0.06) and 1.90 mm (±0.5), respectively, were found for data-set 1. This result was obtained using the stricter event selection (Case B) to warrant canopy pre-dryness. Seven out of the 17 events in data-set 1 had sufficient data points in the wetting phase to perform the respective regression analysis. One example is shown in Fig. 2c.

3.4 Energy balance closure and performance of the sonic anemometer

The slope of the linear regression of H + λE (from eddy covariance) vs. R_n − G₀ − Q for wet and dry half-hour periods was 0.96, while the energy balance ratio (EBR) defined as the sum of H + λE divided by R_n − G₀ − Q was 0.98. The RMSE of the regression was 66.3 W m⁻² for the 30 min interval values of R_n, H, λE, G, and Q (Fig. 3a).

We studied 124 half-hour periods for wet-canopy conditions P_G > 0.5 mm for the full study period from 19 June 2015 to 31 October 2016. In addition to the eddy-covariance data quality flag filtering (see Sect. 2.2.5), we tested the performance of the sonic anemometer by plotting the standard deviation of the vertical wind speed against the friction velocity. According to the Monin–Obukhov similarity theory, the ratio $\mathit{\sigma }w/u^{*}$ should be constant in neutral conditions. We found that the plot presents a strong linear relation (Fig. 3b). The slope was consistent with previous estimations (van der Tol et al., 2003), and the offset was very close to zero.

3.5 Wet-canopy evaporation rates

The wet-canopy evaporation rates derived from the water balance approach were estimated from the parameter $\overline{E}/\overline{R}$ obtained after fitting the linear regressions. $\overline{E}/\overline{R}$ values of 0.25 (±0.02) and 0.23 (±0.03) were found for Case A and Case B, respectively. Because they were similar, we decided to use Case B (stricter event selection) to derive $\overline{E}$. The parameter $\overline{E}/\overline{R}$ multiplied by the median R of 0.82 mm h⁻¹, results in a water-balance based estimated evaporation rate of 0.19 mm h⁻¹.

Figure 3 (a) Half-hour interval of turbulent heat fluxes (H and λE) vs. available energy (R_n − G − Q) for the study site. The solid line represents the 1 : 1 line and the dashed line represents linear regression forced through the origin. (b) Half-hour averages of standard deviation of the vertical wind speed σw (m s⁻¹) vs. friction velocity u^* (m s⁻¹), wet-canopy conditions P_G > 0.5 mm (30 min)⁻¹, and near-neutral stability (−0.02 < (z − d)∕L < 0.02).

Average micrometeorological characteristics of the wet periods are shown in Table 4. It is noteworthy that both sensible heat flux (H) and energy storage (Q) were strongly negative during wet periods. This implies a strong cooling of the surface (Q), accompanied by a negative (downward) sensible heat flux. This is the case for wet periods both during the day (07:00 to 19:00 UTC+1) and the night (19:00 to 07:00 UTC+1). The energy balance residual, which we assume is the latent heat flux (see Eq. 4), greatly exceeds the net radiation. Other observations were a very low vapour pressure deficit at the top of the tower and similar average wind speed throughout the day and night (Table 4).

Between 19 June 2015 and 31 October 2016 (excluding the winter season from November 2015 to March 2016), the mean wet evaporation rate calculated from the energy balance residual, E_EB−EC, was 0.28 mm h⁻¹ during the day (Fig. 4a), 0.07 mm h⁻¹ during the night (Fig. 4b), and 0.20 mm h⁻¹ for day and night periods combined (Fig. 4c). The main sources of evaporation heat (in the equivalent water depth unit) were net radiation (0.07 mm h⁻¹), sensible heat (0.09 mm h⁻¹), and the release of stored energy within the canopy (0.03 mm h⁻¹). Evaporation rates derived from the energy balance residual (for day and night periods) presented a skewed distribution with values ranging from −0.53 to 2.59 (mm h⁻¹) with a median value of 0.13 (mm h⁻¹).

Table 4 Average micro-meteorological characteristics for half-hour periods with more than 0.25 mm (30 min)⁻¹ of P_G for day (07:00—19:00 UTC+1) and night conditions (19:00–07:00 UTC+1).

Figure 4 Distributions of wet-canopy evaporation rates during daytime (07:00–19:00 UTC+1), nighttime (19:00–07:00 UTC+1), and combined day and night. Two different methods applied: (a–c) energy balance residual (E_EB−EC) and (d–f) Penman–Monteith (E_PM−EC).

In order to estimate average wet-canopy evaporation with the Penman–Monteith equation, we first estimated aerodynamic conductance to momentum g_a,M (m s⁻¹). We estimated the aerodynamic conductance to momentum for the predominant south-westerly wind direction and selected the fluxes coming from the wind direction between 180 and 360^∘. In this direction, effects of the tower construction on the wind were minimal. The $g_{\mathrm{a},\mathrm{M},\mathrm{EC}}$ was estimated by means of the regression between wind speed and friction velocity as $g_{\mathrm{a},\mathrm{M},\mathrm{EC}}$ = 0.0318 u (Fig. 5). When we applied the stability correction for non-neutral hours (Paulson, 1970) and used ln(z_0M∕z_0H) = 2 (Lankreijer et al., 1993; Moors, 2012), we obtained an aerodynamic conductance to water vapour as $g_{\mathrm{a},\mathrm{H},\mathrm{EC}}$ = 0.0303 u.

For the whole study period, using the estimated $g_{\mathrm{a},\mathrm{H},\mathrm{EC}}$, and considering Q and G to be part of the available energy (A = R_n + Q + G), the mean and median wet evaporation rates estimated by the Penman–Monteith equation (E_PM−EC) were 0.13 and 0.10 mm h⁻¹, respectively (Fig. 4f, Table 5). The period 19 June–31 October 2015 presented similar mean and median evaporation rates to the values estimated for the period 1 April–31 October 2016 (Table 5).

In general, the mean E estimated with the P–M equation (using $g_{\mathrm{a},\mathrm{H},\mathrm{EC}}$ = 0.0303 u) was 35 % lower than the mean E derived from the energy balance residual. Using the water balance measurements with the mean method resulted in an estimated evaporation rate of 0.19 mm h⁻¹, which is similar to the mean values of E_EB−EC, although 40 % higher than the estimated values of the median E_EB−EC used in the Gash model.

Figure 5 Linear regression of friction velocity u^* against wind speed u for near-neutral hours (−0.02 < (z − d)∕L < 0.02) and from a south-westerly wind direction.

3.6 Modelling rainfall interception

We used two different parametrizations of the Gash model. In the first parametrization (Run 1) the parameters S, p, and $\overline{E}/\overline{R}$ were derived by using the mean method. In the second parametrization (Run 2), the parameters S and p were derived from the IEA; the parameter $\overline{E}/\overline{R}$ was calculated with $\overline{E}$ as the median E derived from the energy balance residual and $\overline{R}$ as the median rainfall intensity derived from the tipping bucket measurements (Table 6).

Table 5 Summary statistics for the wet evaporation rates estimated for the study period by different methods: energy balance (${\overline{E}}_{\mathrm{EB}-\mathrm{EC}}$) and Penman–Monteith equation (${\overline{E}}_{\mathrm{PM}-\mathrm{EC}}$).

^* “All” refers to data from both periods together: 19 June to 31 October 2015 and 1 April to 31 October 2016.

Table 6 Comparison of the performance of modelled interception loss using different parametrization. Data-set 1 refers to the period from 19 June to 31 October 2015, and data-set 2 to the period from 1 April to 31 October 2016. Run 1, all parameters derived from the mean method. Run 2, canopy parameters (S, p) derived from IEA and E from the energy balance residual method.

Table 7 Components of interception loss in mm (and as percentage of total) for data-set 2 (19 June to 31 October 2016) based on the validated Gash analytical original model.

Run 1 underestimated the interception loss by 3 % with an RMSE of 0.93 mm for the calibration data set (Table 6). The model performance based on the Nash–Sutcliffe (NSE) model efficiency was very good (0.90) (Table 6). Run 2, which used the median value of E_EB−EC (0.12 mm h⁻¹) and the median R (0.82 mm h⁻¹) to obtain parameter $\overline{E}/\overline{R}$ (0.15), underestimated I by about 8 % (166.6 mm modelled, vs. 180.4 mm measured I). The RMSE was 1.36 mm and the performance was lower than that of Run 1 (NSE = 0.79). The predicted total interception loss for the validation data set using both Run 1V (V for the validation data set) and Run 2V configurations were in good agreement with I derived from the Tf and Sf measurements, with relative errors of about 1 %. In both cases, the RMSE was about 0.79 mm and the model performance was good, with an NSE of 0.79.

The interception components that contributed most to the overall evaporation interception loss differ between the two parametrizations of model validation, Run 1V and Run 2V. In the first case, the two largest contributors were evaporation from the saturated canopy during rainfall (37 %) and evaporation loss during the drying phase (34 %). In Run 2, the same two components were the greatest contributors, but in opposite order: evaporation loss during the drying phase was the main contributor (44 %), and evaporation from the saturated canopy during rainfall was the second contributor (24 %). The third largest contribution, in both cases (Run 1 and Run 2), was evaporation from small rainfall events (18 and 24 %, respectively), followed by the contribution from the wetting phase (7 and 5 %, respectively). Evaporation from trunks (during saturated and unsaturated conditions) contributed less than 5 % for both data sets (Table 7).

Figure 6 Sensitivity analysis of the parametrized original Gash model. Run 1, all parameters derived from the mean method. Run 2, canopy parameters (S, p) derived from IEA and E from the energy balance residual method. Contour lines representing the RMSE for different combinations of the parameters' canopy storage capacity (S) and the ratio $\overline{E}/\overline{R}$. (a) Sensitivity analysis using calibration data-set 1 (19 June to 31 October 2015. (b) Sensitivity analysis using validation data-set 2 (19 June to 31 October 2016). The red circles represent the corresponding parameters used in the model Run 1 and Run 2.

The sensitivity analysis of the Gash model shows that parameter equifinality (Beven, 1993) occurs between S and $\overline{E}/\overline{R}$ (van Dijk et al., 2015), which implies in this case that an underestimation of S is likely to be compensated by overestimation of $\overline{E}/\overline{R}$ (Fig. 6). This effect can be seen when modelling the validation data set (Fig. 6b): for Run 1 a low value of S (1.37 mm) may be compensated by a high value of $\overline{E}/\overline{R}$ (0.23), leading to a similar RMSE to Run 2, which used a higher value of S (1.90 mm) and a lower value of $\overline{E}/\overline{R}$ (0.15). A similar effect is detected when modelling data-set 1 (Fig. 6a): both parametrizations (Run 1 and Run 2) produce a relative error lower than 10 % (Table 6).

4 Discussion

4.1 Canopy storage capacity

In the same study area, Klaassen et al. (1998) evaluated the most common indirect methods to derive S from multi-event throughfall measurements. They found that the mean method tended to underestimate S compared to direct microwave transmission measuring. However, indirect methods are still largely used due to their low cost and simplicity. As an alternative to the multi-event methods, Link et al. (2004) proposed an analysis of individual events. They found that the assumption of a constant $\overline{E}/\overline{R}$ during multiple events was unsustainable, especially during the wetting phase, and might contribute to the underestimation of S in the mean method.

Our findings confirm that S estimated with the mean method is lower (−30 % lower) than S estimated with the IEA. To avoid underestimation caused by incorporating events not preceded by canopy dryness in the regression analysis, we made use of wetness sensors. However, our findings indicate that, despite using the wetness sensors to eliminate certain events, the results remained similar.

The storage capacity at the study site had been reduced compared to in earlier studies, very likely due to the decrease in tree density and LAI over the years. When the stand was 29 years old, the tree density was 992 trees ha⁻¹ and LAI was 8 (Table 2), while Klaassen et al. (1998) used the principle of microwave attenuation to determine an S of 2.4 mm. At the time of the present study, the tree density as well as the LAI had decreased by about 40 % (Table 2). However, the average total storage capacity (S + S_t) for the study period was reduced much less, only by about 20 %, if we considered S (2.0 mm) derived with the IEA method.

Table 8 Summary of canopy properties and interception parameters for Douglas fir forests.

^a Mixed Douglas fir and Western hemlock; ^b Klaasen et al. (1998) reported an LAI measured by the destructive method, but LAI estimated with Li-COR-2000 was 8 m² m⁻² (cf. Heij and Schneider, 1991). NA = not available.

Our estimation of S is comparable with that of other old Douglas fir stands (Table 8) and supports the hypothesis that LAI might not be the main predictor of S for Douglas fir forests under temperate climatic conditions. Rutter et al. (1975) reported a total storage capacity of 2.1 (S equal to 1.2 mm and S_t of 0.9 mm) for Bramshill Forest (UK), a 42-year old Douglas fir stand with similar density but larger LAI (660 trees ha⁻¹, LAI = 12). In contrast, Link et al. (2004) applied IEA in a 500-year old mixed Douglas fir and Western hemlock forest (560 trees ha⁻¹, LAI = 8.6) located in Washington (USA), and found larger values of S ranging from 2.71 to 4.17 mm. Pypker et al. (2005) in southern–central Washington found significant differences in S for a young (25-year old) Douglas fir forest and an old-growth (> 450-year old) mixed Douglas fir and Western hemlock forest, with S-young being 1.4 mm and S-old-growth being 3.3 mm. Both forests had a similar LAI (LAI-young = 10.2; LAI-old-growth = 9.6); however, despite a notable difference in the stem density (young: 2200 tree ha⁻¹, old-growth: 441 tree ha⁻¹), the larger S was found in the old-growth forest. The high S was presumably caused by the changes in species composition (i.e. presence of understory) and the presence of epiphytes, conditions that were not observed at our study site.

Some authors have found that S can be linearly related to the fraction of vegetation cover, which implies an exponential relation between S and LAI (Moors, 2012). However, for closed canopies (LAI > 5) in temperate climate utilizing the fraction of vegetation cover might not be an option (Moors, 2012). The relation between S and the number of trees and their basal area has also been noted to be valid for several sites (Turner and Lambert, 1987; Teklehaimanot et al., 1991). We speculate that for pine species the tree density and basal area imply a direct relation with the amount of bark present in the forest. This amount will increase as the canopy gets older and taller. Recent research in cedar trees in Japan has found high values of S, and has shown the bark on the stems providing a major contribution (Iida et al., 2017). We only found a value of S_t of 0.14 mm, but part of the trunks' storage capacity may be hidden in the estimated S.

Our results suggest that a long-term decrease in S in a Douglas fir stand does not necessarily imply a decrease in I because the natural process of tree growth as well as forest management practices, such as thinning, all affect other variables that influence the rainfall interception process. Direct extrapolation of I by means of LAI without considering other driving forces (i.e. aerodynamic conductance or change in energy storage) can lead to erroneous approximations of interception loss.

4.2 Wet-canopy evaporation rate

Previous investigations have shown that it is possible to estimate sensible heat flux from the sonic anemometer during rainy conditions. Latent heat fluxes derived from the energy balance residual have been shown to be a good approach to derive evaporation rates during rainfall. However, discrepancies with the water budget approach and the Penman–Monteith equation have been reported in several studies (Ringgaard et al., 2014; Schellekens et al., 1999; Holwerda et al., 2012).

The value of $\overline{E}$ = 0.20 mm h⁻¹ from the present study was derived from the energy balance residual and is in agreement with the canopy structural changes that occurred in Speulderbos due to natural growth and to management practices causing reduced density. At the Speulderbos study site, when the stand was 29 years old (when canopy height was 18 m and LAI was 8 m² m⁻²), Klaassen et al. (1998) used a combination of eddy correlation and psychrometer profile measurements, and reported a wet-canopy evaporation rate of 0.077 mm h⁻¹ (55 W m⁻²). This value was lower than evaporation rates derived at other coniferous forests of the same height and stand configuration and in similar climatic conditions. For example, in the Hafren forest (UK) Stewart (1977) used the Bowen ratio method and found a value of 0.19 mm h⁻¹ from daytime measurements, while Gash et al. (1980) used the Penman–Monteith equation and found a similar evaporation rate of 0.13 mm h⁻¹. Link et al. (2004) reported an average evaporation rate of 0.14 mm h⁻¹ using the P–M equation in a 500-year old mixed Douglas fir forest (60 m height) (Table 8). In a mixed plantation in western Denmark, Ringgaard et al. (2014) found a wet-canopy evaporation rate of 0.21 mm h⁻¹ during the summer season.

In a long-term comparison done by Pypker et al. (2005), similar values of wet evaporation rates were found (young forest: 0.25 mm h⁻¹, old-growth forest: 0.21 mm h⁻¹). However, the studied old-growth forest was a mixture of Douglas fir and Western hemlock with understory present, and the young stand had a very high tree density (Table 8). Our study site offered the advantage of an unchanged species composition, which allowed us to focus on the long-term effects of the changes in tree density and LAI. Over the past 25 years, tree density and LAI at our study site mainly declined as the result of thinning practices, and to a lesser degree due to natural tree mortality. Moreover, forest height increased by about 16 m. The combination of these changes resulted in an increase in aerodynamic conductance from 0.065 to 0.105 m s⁻¹. This change has a direct impact on the exchange of fluxes between the canopy and the atmosphere (Holwerda et al., 2012; Schellekens et al., 1999; Moors, 2012).

$\overline{E}$ (mean) estimated by means of the Penman–Monteith equation was about 35 % lower than with the energy balance residual approach. Several explanations have been proposed in the literature to explain such discrepancies in similar studies. In a comprehensive analysis, van Dijk et al. (2015) pointed out the following possible reasons: (i) underestimation of biomass and heat ground release; (ii) underestimation of aerodynamic conductance; (iii) unaccounted energy advection; (iv) errors in air humidity measurements; (v) mechanical water transport.

In our analysis, we have incorporated estimations of Q and G in our estimate of the available energy A in the P–M equation, and therefore disregard the underestimation of biomass and heat ground release as a main cause of the underestimation of E. However, we have to consider the uncertainty in the estimated Q_bio. This uncertainty, was calculated as the quadratic sum of the relative errors δm_bio, δc_v, and δΔT_bio. The δm_bio is related to the uncertainty of allometric equations (18 % for n = 23, Chave et al., 2004) in combination with the uncertainty in the assumed moisture content (ranging from 44 to 55 %). The combined uncertainty for δm_bio would be about 27 %. Regarding c_v, the range of values used in studies with similar species is from 2400 to 2928 (J kg⁻¹ K⁻¹) (Oliphant et al., 2004), which means a δc_v of 22 %. ΔT_bio, here assumed to be equal to ΔT_air has the largest uncertainty. Based on data presented by Meesters and Vugts (1996) (their Fig. 6) the difference in temperature amplitude between T_air and T_bio would yield to an uncertainty of about 40 %. The error propagation of the product of the three variables in Eq. (3) yields an uncertainty for Q_bio of 53 %.

In the estimation of the aerodynamic conductance, we used g_a,H, calculated using friction velocity derived from the eddy-covariance system, corrected for stability (Paulson, 1970). van Dijk et al. (2015) pointed out that g_a,H and g_a,M require the validity of the Monin–Obukhov similarity theory (MOST). MOST assumes that the turbulent flow (u^*) is the only important velocity scale, and that the height above zero-plane of displacement (z − d) and Obukhov length (L) are the only important length scales (van Dijk et al., 2015). In the case of tall canopies, as at the Speulderbos site, the observations are taken in the roughness sub-layer. In this layer, the turbulence is also influenced by length scales related to the surface (i.e. mixing layer instability at the top of the canopy). Corrections for this effect in MOST during rainfall have not been developed and could lead to an underestimation of aerodynamic conductance. Energy advection has also been hypothesized to be a source of unaccounted additional energy (Shuttleworth and Calder, 1979). Although it appeared to occur mainly at maritime sites (Schellekens et al., 1999), Stewart (1977) advocated that the energy does not necessarily need to come from the ocean, and could be supplied by drier and warmer nearby areas. Although vertical energy advection would not invalidate the P–M equation, horizontally advected energy occurring below the level of energy balance measurements would not be accounted for, and could influence the underestimation of E. Our results suggest that vertical sensible heat flux measured as negative H is the main source of energy that sustains E during rainfall (Table 4); however, this situation was not predicted by the P–M equation. This could be attributed to errors in air humidity measurements. Wallace and McJannet (2008) estimated that 2 % overestimation in RH already leads to an underestimation of E_PM of 36 %. This situation may occur in our case as the accuracy of our RH sensor (CS215, Campbell Sci. Inc.) during wet conditions (> 90 % RH) is low (±4 %). Likewise, enhanced evaporation of rain droplets splashed from the tree canopy has been mentioned as a possible mechanism allowing high interception losses (Murakami, 2006), but given the low rainfall intensities prevailing in the study area, this is not likely to be important.

4.3 Rainfall interception

The interception loss derived for the two consecutive growing seasons of 2015 and 2016 was 37 and 39 % of gross rainfall (P_G), respectively. These values are in agreement with other similar studies (Table 8). For similar climatic conditions, Rutter et al. (1971) investigated a Douglas fir stand of similar age and stem density in Bramshill Forest (UK) and found an I of 39 % of P_G. Soubie et al. (2016) found a value for I of 30 % in a Douglas fir stand in Belgium, which is lower than the value found in the current study and may be attributed to the difference in stem density (145 tree ha⁻¹) since LAI was similar (LAI = 4.2). In a comparison between mixed young and old Douglas fir stands in the north-western Pacific (Washington, US), Pypker et al. (2005) reported similar values of I, namely 21 and 24 %, respectively. They attributed the slightly larger I in the old stand to the higher S which was also linked to the change in species composition and to the presence of epiphytes.

The other interesting finding in the study by Pypker et al. (2015) is that $\overline{E}/\overline{R}$ was similar for both stands (i.e. young and old). Because the older stand was taller, the aerodynamic conductance may have been larger, with larger expected evaporative rates as a consequence (Teklehaimanot et al., 1991). In contrast, and considering that $\overline{R}$ was not variable (the two stands were close to each other), the evaporative fluxes from the wet canopy were similar. Pypker et al. (2005) observed that the variable species composition at their study site (i.e. Douglas fir mixed with Western hemlocks) increased the gap size and influence on g_a. They explained that in this particular situation the use of the average Douglas fir height was likely inappropriate for calculating d and z₀ and hence $\overline{E}/\overline{R}$.

In the case of Speulderbos at a younger stage, the larger stem density and higher LAI meant a larger S (Klaassen et al., 1998), while, at an older stage (2015–2016), the $\overline{E}/\overline{R}$ was larger mainly as an effect of a larger g_a, due to the larger canopy height in combination with lower tree density (Teklehaimanot et al., 1991).

The original version of the Gash (1979) analytical model successfully predicted I for the calibration and validation data sets (Table 6). Several studies have demonstrated the validity of the Gash model in temperate coniferous forests (Muzylo et al., 2009).

The performance of the Gash model was as good as that of some more sophisticated models applied in earlier studies at the same site (i.e. Bouten et al., 1996). The difference between observed and predicted values of I was lower than in previous applications of the Gash model in similar climatic conditions (Lankreijer et al., 1999). The model overestimated the total I for all parametrizations (Table 6). During the calibration, the relative error of total I was lower than 8 %, and for the validation phase, it was lower than 2 %.

We used the RMSE as a criterion to evaluate the sensitivity of the model to the two main parameters in the Gash model, S and $\overline{E}/\overline{R}$, considering that the other independently derived parameters have less influence on the model and can be kept fixed. We observed that for calibration (Fig. 6a) and validation (Fig. 6b) certain combinations of S and $\overline{E}/\overline{R}$ yield an almost equally low RMSE. Such combinations along the major axis of the inner ellipse that encompasses the optimal solution represent the above-mentioned functional equivalence between S and $\overline{E}/\overline{R}$. For instance, in Fig. 6b, parametrization for Run 1 presents a higher $\overline{E}/\overline{R}$ (and lower S) than the one set for Run 2, with a similar RMSE of below 0.8 mm. This confirms that it is necessary to reduce the uncertainty in one of the parameters (S or $\overline{E}$) by independent measurements, before optimizing the second one (i.e. Gash et al., 1980; Ringgaard et al., 2014). The independent estimations of S and $\overline{E}$ by means of the IEA and the energy balance residual give us the confidence to select the parametrization of Run 2 as the most realistic in the case of Speulderbos.

5 Conclusion

Rainfall interception loss (I) was quantified and modelled for a mature Douglas fir stand (ca. 55 years old) in the centre of the Netherlands for two growing seasons in 2015 and 2016. Over the study period, the interception loss was 38 % of gross precipitation. This value was similar to the value reported for the same stand when the forest was 29 years old, indicating the changes in forest structure may not always result in changes in interception loss. Tree density as well as LAI were reduced by about 40 % in comparison with the former study by Klaassen et al. (1998). However, the change in canopy storage capacity (S) was much less (a reduction of about 20 %). Canopy storage capacity remained relatively stable largely due to the increase in the total biomass, and more specifically in stem bark surface. The reduction in stem density and the growth of canopy height resulted in a larger surface roughness and in consequence enhanced the evaporation rate during rainfall.

The main sources of energy supply that sustain evaporation of intercepted rainfall were net radiation (35 %), sensible heat flux (45 %), and change in energy stored in air and canopy biomass (15 %). Downward sensible heat fluxes estimated by means of the eddy-covariance technique were larger than those predicted by the P–M equation, possibly due to inaccuracies in the measurement of the relative humidity in the air.

The Gash model was able to simulate I reasonably well, with relative errors of less than 10 %. A sensitivity analysis of interception losses simulated with the Gash model shows that the presently higher $\overline{E}/\overline{R}$ can indeed compensate for the lower S, confirming the parameter equifinality effect.

Our results confirm that even after a reduction in tree densities old growth stands can maintain similarly high rates of interception. This finding will be useful to improve long-term model predictions that involve structural changes or planned management practices in forested ecosystems.

Appendix A

Table A1 Statistical description of collection periods of throughfall and average amounts for a sample size n = 32.