Forest canopy water fluxes can be estimated using canopy structure metrics derived from airborne light detection and ranging (LiDAR)

Agricultural and Forest Meteorology 203 (2015) 131–141

Contents lists available at ScienceDirect

Agricultural and Forest Meteorology journal homepage: www.elsevier.com/locate/agrformet

Forest canopy water fluxes can be estimated using canopy structure metrics derived from airborne light detection and ranging (LiDAR) Johannes Schumacher a,∗ , Jesper Riis Christiansen a,b a b

University of Copenhagen, Department of Geoscience and Natural Resource Management, Rolighedsvej 23, DK-1958 Frederiksberg C, Denmark University of British Columbia, Department of Forest and Conservation Sciences, 2424 Main Mall, Vancouver V6T 1Z4, Canada

a r t i c l e

i n f o

Article history: Received 14 March 2014 Received in revised form 31 October 2014 Accepted 15 December 2014 Available online 16 January 2015 Keywords: Ecosystem fluxes Laser scanning Remote sensing Precipitation Canopy rain interception Throughfall

a b s t r a c t Forests contribute to improve water quality, affect drinking water resources, and therefore influence water supply on a regional level. The forest canopy structure affects the retention of precipitation (Pr) in the canopy and hence the amount of water transferred to the forest floor termed canopy throughfall (TF). We investigated the possibilities of estimating TF based on bulk Pr and canopy structure estimated from airborne light detection and ranging (LiDAR) data. Bulk Pr and TF fluxes combined with airborne LiDAR data from 11 locations representing the most common forest types (mono-species broadleaf/coniferous and mixed forests) in Denmark were used to develop empirical models to estimate TF on a monthly, seasonal, and annual basis. This new approach offers the opportunity to greatly improve predictions of TF on catchment wide scales. Overall, results show that TF can be estimated by Pr and a canopy density metric derived from LiDAR data. In all three types of TF data sets Pr was the variable explaining the majority of the variance in TF. The proportion of explained variance adhering to the LiDAR variable increased from 1.7% for the monthly data set to 12.2% and 19.5% for seasonal and annual data sets, respectively. Although the analysis was limited to Denmark the model was successful in estimating TF for contrasting tree species (broadleaf vs. conifers) and points to a potential for extending our model approach to other similar regions. Our approach can help to assess how forest cover impacts water resources on a large scale in regions where forests play a major role in water resource management. © 2015 Elsevier B.V. All rights reserved.

1. Introduction Trees act as an important link in the forest water balance as they effectively intercept, store, and transfer water across the canopy–soil interface. Precipitation (Pr) is intercepted by tree crowns and is evaporated back to the atmosphere (= interception loss), or is transferred to the ground by falling through canopy gaps or by leaf drop-off (in the following both are referred to as throughfall, TF), or by stemflow. Canopy interception has been found to range from 10 to 50% of gross rainfall (Roth et al., 2007) and the amount of TF is a function of the interception properties of forest canopies (Rutter et al., 1971). Not surprisingly, several studies of forest water balance have concluded, that interception of Pr is a key factor determining water fluxes in forest ecosystems (e.g. Christiansen et al., 2006, 2010, and references therein). Especially, TF water flux is important as it impacts infiltration and availability

∗ Corresponding author. E-mail addresses: [email protected], [email protected] (J. Schumacher). http://dx.doi.org/10.1016/j.agrformet.2014.12.007 0168-1923/© 2015 Elsevier B.V. All rights reserved.

of water in the root zone, groundwater regeneration, and the deposition of nutrients from the atmosphere and canopies to the soil (Kristensen et al., 2004). However, in studies of water balances at larger scale, it is still desired to integrate forest vegetation composition and retain a detailed physical description of the trees as this has a pronounced impact on the forest water fluxes (Carlyle-Moses and Gash, 2011). Thus, it has been proposed that remote sensing data could be used to estimate forest canopy water fluxes (Roth et al., 2007). Using satellite-based remote sensing data it was possible to estimate changes in canopy-interception under various land-uses in Indonesia (Nieschulze et al., 2009). Extending this idea airborne light detection and ranging (LiDAR) offers an excellent opportunity in this regard as it can cover areas relevant to catchment wide studies while retaining details on canopy structure at the stand level. Recently, Mitchell et al. (2012) used airborne LiDAR data to improve estimates of catchment scale evapotranspiration from eucalyptus forests in Australia. A key feature of forest canopy structure is their density and the spatial distribution and configuration of foliage that impacts the amount of precipitation that is captured, stored and transferred

132

J. Schumacher, J.R. Christiansen / Agricultural and Forest Meteorology 203 (2015) 131–141

to the ground surface. In this context LiDAR data used to describe canopy structure can be translated in to a vertical 3D density structure of the forest stand as the laser beam penetrates into the tree crowns and is reflected at different elevations dependent on the density of the canopy cover. Thus, using LiDAR metrics calculated from the returns from forests excluding the ground surface can be used as a measure of the canopy density because the LiDAR integrates the three dimensional configuration of the canopy. Following ˜ et al. (2004), Morsdorf et al. this idea, Lefsky et al. (2002), Riano (2006), Solberg et al. (2009) used LiDAR based variables to estimate effective leaf area index (LAIe, including leaves, branches, and stems), or plant area index (PAI). PAI is as well related to canopy saturation capacity (Rutter et al., 1971; Gash, 1979), and both, PAI and the canopy saturation capacity, are central in order to understand the process of rainfall interception by forest canopies (Muzylo et al., 2009). Despite interception and TF in forests being studied intensively in the past 40 years e.g. Rutter et al., 1971; Muzylo et al., 2012, direct measurements or estimation of canopy interception parameters is still based on time consuming measurements or calibration of generic models (Rutter and Morton, 1977; Muzylo et al., 2009). Roth et al. (2007) suggested that LiDAR offered an opportunity to estimate key parameters relevant to interception models, such as canopy storage capacity and LAI in forests. Different tree types display different canopy structures which relate to their canopy storage capacities and TF fluxes (Holder, 2013). Following these arguments we hypothesized that canopy density metrics derived from LiDAR scans can be used to estimate canopy water fluxes for different forest types and tree species. Our aims were (1) to test if LiDAR based variables can represent canopy interception abilities and therefore can be used for estimating TF at different temporal scales, and (2) to assess if LiDAR metrics can describe forest stand development and associated water fluxes. We further extended our analysis to produce a map of forest canopy TF based on LiDAR metrics and meteorological variables for a selected forest area (Hov Skov) in south-east Denmark.

2. Materials and methods 2.1. Sites and meteorological data We used TF and bulk Pr data measured monthly between 2005 and 2007 at 11 sites distributed across Denmark (Fig. 1) that included data from 21 plots of approximately 0.15–0.25 ha. These sites represented monospecies stands of one conifer and five broadleaf species and one mixed plot dominated by old broadleaf trees within a relatively narrow age range (Table 1). We included different tree species commonly found in Danish forests in order to represent different canopy and leaf structures. The geographic representation of sites captured a broad range of growing conditions of both high and low productive sites (Gundersen et al., 2009). At all the above mentioned sites the stands had reached canopy cover ≥ 90% before measurements began. The restriction to closed canopy sites is due to the fact that the sites were originally established as monitoring plots for water and nutrient budgets in mature forests in Denmark. Additionally, we added a slightly older TF dataset (2000–2003) from two chronosequences at the two locations Vestskoven-2 and Gejlvang (Fig. 1 and Table 2) to study the relation between stand development, LiDAR metrics and TF. We used data from monocultures of Norway spruce and oak planted in five successive stages since 1969–1970. The stand ages at the Vestskoven-2 chronosequence at the time of LiDAR scan were for oak 13, 18, 27, and 29 years, and for Norway spruce 9, 16, 18, 33, and 37 years. At Gejlvang

the stand ages for Norway spruce at the time of LiDAR scan were 13, 26, 31, and 47 years. Within each stand of one age class three plots were established and measured. Three more stands for which no TF data was available were included to assess the LiDAR variable in stands of different age classes (36 and 206 years old oak stands at Vestskoven-2, and 9 years old Norway spruce stand at Gejlvang). Bulk Pr and TF were monthly sampled in funnels with an internal diameter of 12.8 cm and an area of 128.7 cm2 at 1 m above ground. Five funnels were installed in each plot, and one to two collectors were placed in a nearby open field or clearing of size > 2 ha and in 100–400 m distance for bulk Pr measurement. Each TF funnel was equipped with a 45 ␮m filter placed beneath a sieve to avoid contamination of samples by animals and plant fragments. Filters were replaced at each sampling occasion. Plastic tubes connected bulk Pr and TF funnels to plastic bottles buried in the soil to keep the samples dark and cool to avoid evaporation. At sampling, the volume of collected water was registered for each funnel and converted to L m−2 or mm. The monthly sums of each funnel were averaged for each plot and month, resulting in 268 Pr and 374 TF observation in total which were used for modelling. Monthly bulk Pr and TF for all sites and all tree species used for modelling ranged from 4.82 to 208.6 mm and from 1.41 to 170.97 mm, respectively (Table 1). The maximum capacity of the water collecting buckets was 169.23 mm. On twelve sampling occasions the buckets were found full (two at Als, one at Gludsted, one at Hald Ege, two at Lindet, three at Suserup, and three at Vestskoven-1). Four of these values were replaced with modelled bulk precipitation (Prm ) provided by the Danish Meteorological Institute (DMI), resulting in maximum values of 200.3 mm at Als, 174.8 mm Lindet, 208.6 mm Suserup, and 172.9 mm Vestskoven-1. At the remaining eight occasions, the modelled DMI values were smaller than the maximum measured values which therefore were kept. Thus, two sites in Table 1 (Gludsted and Hald Ege) show the same maximum Pr value. Some measurements were missing and were therefore interpolated based on other measurements and modelled Prm provided by DMI. This resulted in maximum Pr of 199.76 mm for spruce, and maximum TF of 170.97 mm for oak at Thyregod, and maximum Pr of 178.79 mm for oak at Vallø (Table 1). Maximum TF at Vestskoven-2 for oak (230 mm), and maximum Pr at Gejlvang (244 mm) (Table 2) are most likely measurement errors. Errors and uncertainties are part of modelling and the reason for variation, and therefore we kept these values and did not change anything in the data. The monthly Pr and TF for each funnel were summed and averaged to two seasons (winter (November–April) and summer (May–October)) and on an annual basis. The DMI provided monthly climate data (air temperature, wind speeds, global radiation, and potential evapotranspiration, modelled for a 20 km by 20 km grid, and bulk precipitation modelled for a 10 km by 10 km grid (Prm )) based on data from climate stations distributed all over the country (Scharling, 2012). Two cells of Prm data were used for modelling TF for a selected forest area (Hov Skov) in southern Zealand. The mean of the monthly Prm of 2007 was used, since that year corresponded to the time the LiDAR data was collected for this site. 2.2. LiDAR data The LiDAR data was collected using a density of 0.5 pulses per m2 during leaf-off season in spring/fall 2006 and in spring 2007 for the entire country to calculate a digital terrain model for Denmark. The data was collected using a fixed-wing aircraft and preprocessed by a private contractor (COWI, 2007a,b). First and last returns were recorded with a density of 0.5 pulses per m2 , and each return was classified into the classes (1) last echo on vegetation, (2) last echo on terrain, or 3) first echo on terrain or vegetation. Flight parameters and details about the scanning survey are shown in Table 3.

Table 1 Locations used for model development with tree species, age, stem no. (N) and basal area (G) per ha, leaf area index (LAI), max. tree height (h) estimated from LiDAR, measurement periods, no. of canopy throughfall (TF) measurements (= mean of 5 funnels) (x) per plot, coefficient of variation (CV) and confidence intervals (C.I.), and range of bulk precipitation (Pr) and TF. Tree species5

N per ha

G per ha

LAI6

Max h LiDAR

Date (LiDAR)

Dates (TF, Pr)

x per plot

CV [%]

95% C.I.

Als1 Frederiksborg1 Gludsted1 Hald Ege1 Lindet1

Beech Beech Spruce Oak Spruce

110 41 44 147 38

na 721 na na na

na 28.3 na na na

3.4 5.6 3.8 3.0 6.7

31.7 22.7 18.5 24.0 24.0

01.05.’06 07.05.’06 22.03.’07 27.04.’06 02.04.’07

01/’05–12/’07 01/’06–12/’07 01/’05–12/’06 01/’05–12/’07 01/’05–12/’06

35 22 21 35 23

23.9 18.5 27.3 23.4 61.4

20.7, 27.1 16.9, 20.1 23.9, 30.7 16.5, 30.3 49.0, 73.8

Lime Spruce Ash Beech Maple Oak Mix4

34 34 34 34 34 34 ∼5–200

940 1158 989 1130 992 1047 na

31.0 51.3 23.9 24.5 29.8 16.8 na

6.1 5.4 4.2 4.2 3.6 3.4 5.7

20.0 22.1 19.2 13.5 17.1 12.7 31.6

15.6 28.4 18.8 60.2 16.7 58.1 20.8

12.4, 18.7 23.6, 33.1 14.5, 23.1 47.4, 72.9 13.8, 19.7 46.1, 70.1 16.9, 24.6

Spruce Oak Spruce

26 ∼100–200 42

na na 1303

na na 39.0

6.2 3.2 6.5

15.4 20.4 19.2

7.1 11.8 30.4

4.6, 9.6 8.2, 15.4 27.5, 33.3

Oak Maple Spruce Lime Beech Oak

34 34 34 34 34 35

653 933 883 902 729 na

17.6 19.3 45.6 28.7 24.6 na

3.1 4.8 5.9 6.1 5.3 3.0

16.8 17.1 21.1 17.0 17.6 17.4

13.3 16.1 28.4 17.7 16.1 17.5

12.2, 14.5 13.9, 18.3 23.2, 34.5 10.3, 25.2 14.1, 18.0 15.7, 19.3

Age [years]

Range of values Pr [mm]

Mattrup2

Suserup1 Thyregod3 Ulborg1

Vallø2

Vestskoven-11 1 2 3 4 5 6

26.03.’07

01/’06–09/’06

13.03.’07

01/’05–12/’07

12.03.’07

08/’05–07/’07

10.02.’07

01/’06–12/’07

03.01.’07

01/’06–08/’06

26.11.’06

02/’05–12/’07

9

30 25 25

9

34

Data from Gundersen et al. (2009). Data from Christiansen et al. (2010). Unpublished. Dominated by old-growth beech trees interspersed with ash. Beech (Fagus sylvatica L.), Norway spruce (Picea abies (L.) Karst.), Oak (Quercus robur L.), Lime (Tilia cordata), Ash (Fraxinus excelsior), Maple (Acer pseudoplatanus). LAI measured as described in Christiansen et al. (2010).

10.77–200.30 14.52–165.75 20.18–169.23 21.15–169.23 9.54–174.80

30.07–134.19

6.62–208.60 21.15–199.76 34.30–157.08

4.82–178.79

7.10–172.9

TF [mm] 5.04–150.27 11.00–138.91 12.77–133.37 12.70–137.54 3.59–95.37 13.54–92.55 6.99–69.09 16.84–96.39 9.30–91.04 16.98–93.79 9.32–98.75 5.09–167.34 7.48–133.87 11.75–170.97 9.39–155.35 2.99–146.03 4.10–157.75 1.41–110.40 2.97–139.32 2.57–120.76 8.90–161.03

J. Schumacher, J.R. Christiansen / Agricultural and Forest Meteorology 203 (2015) 131–141

Location

133

134

J. Schumacher, J.R. Christiansen / Agricultural and Forest Meteorology 203 (2015) 131–141

Fig. 1. Map of Denmark showing the locations used in this study. The sites represent long- and short-term forest monitoring plots with concurrent precipitation, canopy throughfall (TF), and LiDAR data. At Hov Skov the model was applied to produce a TF map.

Table 2 Locations of the chronosequence plots with tree species, age, max. tree height (h) estimated from LiDAR, measurement periods, no. of collected canopy throughfall (TF) measurements (= mean of 5 funnels) (x) per plot (three replicates per plot), and range of bulk precipitation (Pr) and TF, to asses the effect of stand age. Location

Tree species

Oak Vestskoven21 Spruce

Gejlvang1

1

Spruce

Data from Rosenqvist et al. (2010).

Age [years]

Max h LiDAR

13 18 27 29 9 16 18 33 37

6.9 11.0 13.5 13.3 7.1 10.4 12.6 16.2 20.4

13 26 31 47

5.5 14.5 17.3 20.0

Date (LiDAR)

Dates (TF, Pr)

x per plot

Range of values Pr [mm]

TF [mm]

26.11.2006

10/’00–12/’02

81 81 81 81 63 81 80 81 81

8–141

5–122 4–133 3–136 4–230 3–131 2–143 2–140 1–158 1–128

12.03.2007

04/’01–03/’03

68 72 69 71

18–244

4–143 2–121 2–149 3–148

J. Schumacher, J.R. Christiansen / Agricultural and Forest Meteorology 203 (2015) 131–141

135

Fig. 2. Example of the vertical distribution of LiDAR returns from three forest stands (Lime (Tilia cordata), Norway spruce (Picea abies (L.)), and Maple (Acer pseudoplatanus)) at Vallø, viewed from the side for a ∼15 m wide strip. Table 3 Summary of the laser scanner and flight data. System Flying altitude Pulse repetition frequency Scanning angle Average point density Footprint size Horizontal accuracy Vertical accuracy

Optech ALTM 3100 1600 m 70 kHz ±24◦ 0.5 pulse/m2 50 cm 80 cm (1) 10 cm (1)

The resulting point cloud was tested for geometric accuracy and a digital terrain model (DTM) was calculated using the TerraScan software. Based on the DTM, height above ground was calculated for each point (Fig. 2). The LiDAR data we used only captured one state of forest stand development, but the TF dataset used for model building was collected over three years for some of the plots. We assumed that during these three years the trees did not grow to an extent as to have any significant impact on canopy structure and hence TF, since all of these stands had established a canopy cover of ≥ 90% at the time of the measurements. No forest management activities or other disturbances as for example windthrow occurred within this period on the plots used for model development. We extracted all LiDAR returns of class 1 and 3 and calculated a large set of height and density metrics from all returns with height >1 m (Table 4) for all plots (0.15–0.25 ha) using the program Cloud Metrics of the Fusion LDV software (McGaughey, 2012). These metrics describe maximum stand height, vertical distribution of foliage, foliage density in height layers, and canopy thickness in terms of crown length and density. Næsset (2002) related similar LiDAR based metrics to forest stand properties as mean and dominant height, stem number, basal area, and volume. Such properties also contribute to the ability of intercepting precipitation and ˜ et al. (2004), Morsdorf et al. (2006), Solberg therefore to TF. Riano et al. (2009) used similar or related LiDAR variables to estimate LAI, which also is related to rainfall interception and TF in forest stands (Muzylo et al., 2009). In our study, we used raw LiDAR variables as a Table 4 List of LiDAR based variables that were assessed for model development. LiDAR derived metrics

Description

Max Std L-moments (1, 3, 4)

Maximum height Standard deviation of heights product moments corresponding to mean, skewness, and kurtosis (Hosking, 1990) of heights height below which a certain percentage of all returns were reflected ((mean − min)/(max − min)) calculated from the point cloud Percentage of all returns above 1 m Percentage of all returns above mean No. returns per relative height layer 1–10 (1–max height)/total no. of returns

Height percentiles Canopy relief ratio Interception ratio (IR) Interception ratio 2 (IR2) Density metrics

˜ et al. (2004) showed proxy for stand canopy density (Fig. 2), as Riano Pearson correlation coefficients of up to 0.85 between raw LiDAR variables (50th, 75th, 95th height percentiles, mean and maximum height, percentage of canopy hits) and LAI. 2.3. Modelling We used the statistics software R (R Core Team, 2013) to fit mixed linear models using maximum likelihood. Combining the field measurements of square root transformed canopy TF, meteorological data, and independent LiDAR metrics, we found linear mixed regression models to estimate TF on a monthly, seasonal, and annual basis. A starting model with the structure: y = Xˇ + Zu + , with u∼N(0, G) and ∼N(0, R), (1) √ where y = TF, X and Z = design matrices for fixed and random effects, respectively, ˇ = fixed effects parameters, u = collection of all random effects coefficients, G = covariance matrix for random effect coefficients, R = covariance matrix for measurement errors, was stepwise reduced by forward and backward selection of variables based on Akaike Information Criterion (AIC) as stopping rule (stepAIC function in R), and was further reduced by backward selection based on p-values (p< 0.05). The response variable TF was square-root transformed to avoid heteroscedasticity of the model residuals. Transformation bias was corrected during backtransformation by applying the simplified estimator presented in Eq. (14) in Gregoire et al. (2008). Leave-One-Plot-Out-Cross-Validation was applied to assess the performance of the model. All observations of one of the 21 plots at the time were left out and the model was fit with all observations from the remaining 20 plots. Using this model, estimations for the plot that was left out were made. With the differences between these estimated values and the measured values the root-meansquare error (RMSE) was calculated. We tested monthly models for winter and summer separately to identify possible effects caused by the LiDAR data being collected during leaf-off condition and the TF data being collected during all seasons. Additionally, to fit models on a seasonal basis (summer (May–October) and winter (November–April)), the observations for each plot and season were averaged. Correspondingly, the observations for each year were averaged, to fit models on an annual basis. We used the additional dataset from Vestskoven-2 and Gejlvang (including the three additional stands lacking TF data) to assess how a canopy density metric derived from LiDAR behaved in stands of different age classes. Further, we tested how well the annual model performed on stands of age classes different from the ones used for model development. It was, however, difficult to compare the estimated TF and the measured TF, since this additional TF data set was collected 3–6 years before the LiDAR data. We assessed how well the models captured the differences between tree types by evaluating the model residuals for conifers and broadleaves.

136

J. Schumacher, J.R. Christiansen / Agricultural and Forest Meteorology 203 (2015) 131–141

Table 5 Parameter estimates, standard errors, p-values, contribution of each variable to the variation explained by fixed effects, and fixed, random, and residual variance components for the monthly models, and the seasonal, and annual model. Pr = bulk precipitation; Tair = mean modelled air temperature; IR = LiDAR variable “interception ratio” (% returns > 1 m). Intercept

Fixed effects

Random effects

Pr

Pr2

Tair

IR

Monthly model (N = 374), RMSE = 14.95 Estimates 4.3991 Standard error 0.3709 p-Value < 0.0001 Explained variation [%] Variance components [%]

0.0722 0.0039 < 0.0001 96.49 86.50

−0.0001 0.0000 < 0.0001 1.60

−0.0172 0.0070 0.0144 0.23

−0.0249 0.0038 0.0001 1.68

Monthly model winter (N = 183), RMSE = 16.89 Estimates 4.4093 Standard error 0.7474 p-Value < 0.0001 Explained variation [%] Variance components [%]

0.0632 0.0065 < 0.0001 97.26 74.76

−0.0001 0.0000 0.0014 1.82

Monthly model summer (N = 191), RMSE = 14.07 Estimates 3.9379 Standard error 0.3715 p-Value < 0.0001 Explained variation [%] Variance components [%]

0.0756 0.0042 < 0.0001 97.53 92.98

−0.0001 0.0000 < 0.0001 1.07

Seasonal model (N = 81), RMSE = 14.04 Estimates 5.8346 Standard error 0.5850 p-Value < 0.0001 Explained variation [%] Variance components [%]

0.0415 0.0038 < 0.0001 87.76 61.72

Annual model (N = 38), RMSE = 15.65 Estimates 3.5434 Standard error 0.9280 p-Value 0.0009 Explained variation [%] Variance components [%]

0.1067 0.0184 < 0.0001 51.38 66.14

0.00

2.07

11.43

4.52

4.59

16.13

0.02

1.02

5.97

0.00

12.53

25.75

16.87

16.99

−0.0194 0.0083 0.0447 0.91

−0.0251 0.0040 < 0.0001 1.40

−0.0005 0.0001 0.0002 29.10

3.1. Modelling 3.1.1. Monthly model The best model for estimating monthly TF included Pr, Pr2 , mean air temperature (Tair ), and the LiDAR variable percentage of all returns above 1 m, termed interception ratio (IR) (Eq. (2), Table 5).

TF month =  + ˛ · Pr + ˇ · Pr 2 +  · Tair

+ ı · IR + a · Plot + b · Location + 

Residual Location

−0.0242 0.0062 0.0034 12.24

3. Results



Plot

(2)

The RMSE was 14.95 in the original scale (mm). The fixed effects explained 86.50% of the total variation, the random effect plot explained 0%, the random effect location explained 2.07%, and the residual variance was 11.43% (Table 5). The random effect plot did not contribute significantly to the overall variation, but was kept in the model as part of the model design. The standard errors were much smaller than the parameter estimates, and the p-values were significant (Table 5). Pr covered the main part (96.49%) of the variation explained by the fixed effects, while the LiDAR variable IR only covered 1.68%. The RMSE in the original scale for the monthly winter and summer models was 16.89 and 14.07, respectively. The fixed effects for these models explained 74.76 and 92.98% of the total variation, respectively. The random effect plot explained 4.52 and 0.02%, and location explained 4.59 and 1.02% in the winter and summer

−0.0217 0.0058 0.0012 19.52

model, respectively, and the residual variance was 16.13 and 5.97%, respectively (Table 5). The residuals showed no heteroscedasticidity and the histograms of the residuals showed the typical bell shape of a normal distribution for all of the three monthly models (Fig. 3a–c). The back-transformed estimations corresponded to the observed values. removed: up to TF amounts of 125 mm. Above this value all three models tended to underestimate TF. Two estimations around 100 mm TF can be observed, where the monthly and the winter model overestimated observed TF by almost 200% (Fig. 3d and e). The summer model tends to overestimate TF (Fig. 3f). The estimations of the winter and summer models are similar to the ones of the combined monthly model (Fig. 3d–f). 3.1.2. Seasonal model The best model for estimating seasonal TF included only Pr and IR as fixed effects (Eq. (3), Table 5).



TF season =  + ˛ · Pr + ı · IR + a · Plot + b · Location + 

(3)

The RSME was 14.04 in the original scale. The fixed effects explained 61.72% of the total variation, the random effect plot explained 0%, the random effect location explained 12.53%, and the residual variation was 25.75% (Table 5). The parameter estimates and standard errors were similar to the monthly model, and the pvalues were significant (Table 5). Pr and IR accounted for 87.76% and 12.24% of the variation explained by the fixed effects, respectively. The residuals showed no heteroscedasticidity and the histogram of the residuals showed only a minor positive skewness (Fig. 4a). The estimated values corresponded to the observations, except

J. Schumacher, J.R. Christiansen / Agricultural and Forest Meteorology 203 (2015) 131–141

137

Fig. 3. Evaluation of the monthly model (a and d), the monthly winter model (b and e), and the monthly summer model (c and f). The left column (a–c) shows the histograms of the model residuals, and the right column (d–f) the fitted values versus observed values in the original scale.

two outliers, where the model overestimated TF by almost 200% (Fig. 4c). 3.1.3. Annual model The best model for estimating annual TF included Pr, Pr2 , and IR. Only one random effect (geographic location) was included. The random effect plot was not included in the model since the data for the annual model did not provide sufficient complexity to include this additional variance component (Eq. (4), Table 5).



TF year =  + ˛ · Pr + ˇ · Pr 2 + ı · IR + b · Location + 

(4)

The RSME was 15.65 in the original scale. The fixed effects explained 66.14% of the total variation, the random effect location explained 16.87%, and the residual variance was 16.99% (Table 5). The parameter estimates, standard errors, and p-values were similar to the models above. Pr, Pr2 , and IR accounted for 51.38%, 29.10%,

and 19.52% of the variation explained by the fixed effects, respectively (Table 5). The residuals showed no heteroscedasticidity and the histogram of the residuals showed the bell shape indicating normal distribution of the residuals (Fig. 4b). The estimated values corresponded well to the observations (Fig. 4d). We applied the annual model to a selected forest (Hov Skov) in southern Zealand and produced a TF map (Fig. 5) for 50 m×50 m grid cells. Pronounced small scale differences can be observed especially in the north-western and in the south-eastern parts. In general the south-eastern part had higher TF rates. The north-western part consisted of dense forest stands whereas the south-eastern part consisted of more open and heterogeneous stands. TF ranged between 39.24 and 71.16 mm, with a mean of 55.98 mm. The standard deviation was 6.43 mm.

138

J. Schumacher, J.R. Christiansen / Agricultural and Forest Meteorology 203 (2015) 131–141

Fig. 4. Evaluation of the seasonal and the annual model. The left column (a and b) shows the histograms of the model residuals (note that the scales of the y-axis differ), and the right column (c and d) the fitted values versus observed values in the original scale.

3.2. Stand development and LiDAR derived metrics

3.3. Forest types and canopy structure derived from LiDAR

IR increased with stand age in Norway spruce at Gejlvang (Fig. 6(a)). A non-linear asymptotic increase from around 10% returns at the youngest stage (9 years) to 98% at stand age 47 years was observed. For Norway spruce at Vestskoven-2 this trend was not visible, and for oak this trend was weak. Correspondingly, TF seemed to decrease with stand age for Norway spruce at Gejlvang, but hardly for the stands at Vestskoven-2 (Fig. 6(b)). However, this trend was weak and given the error bars a clear statement cannot be made. One reason for these differences at the two locations might be different growing conditions between the Gejlvang stands in the western part of the country with nutrient poor sandy soils compared to the Vestskoven-2 stands in the eastern part of the country with fertile clayey soils. This could lead to increased growth in Vestskoven (note the high interception of LiDAR returns of almost 100% in Norway spruce, Fig. 6(a)). Furthermore, and in contrast to the plots used for model development, the stands we investigated for the analysis in this section have probably been thinned at different growth stages. That would overrule part of the natural development of the canopy structure and could help to explain why there is no trend of canopy interception with time at Vestskoven. We applied the annual model on the forest stands of different age classes (chronosequence) and compared estimated to measured TF (note that there were no TF measurements available for the 36 and the 206 year old oak stands at Vestskoven-2, nor for the 9 year old spruce stand at Gejlvang). The differences between estimated and observed TF were significant on the 0.05 level (pvalue of 0.01) when using the data from each individual plot and were not significant (p-value of 0.13) when averaging over the three replicates of one age class.

Broadleaf stands generally had significantly higher TF rates than conifer stands (p-value < 0.0001 (Welch Two Sample t-test), median of TF/Pr = 0.72 for broadleaves (n = 262) and 0.56 for conifers (n = 112), and 95% confidence interval = 0.7, 0.74 for broadleaves and 0.53, 0.6 for conifers; note that the Pr effect had been removed so that TF normalized with total Pr was used for this comparison). However, the model residuals of broadleaf and conifer plots were not significantly different from each other in the monthly (p-value = 0.35), monthly winter (p-value = 0.94), monthly summer (p-value = 0.17), seasonal (p-value = 0.82), and annual model (p-value = 0.32). 4. Discussion and conclusion We successfully used remote sensing data from LiDAR to derive canopy density metrics and related these to TF water fluxes for deciduous, coniferous, and mixed forest types in Denmark. We know only of two other studies that used remote sensing in relation to TF (Nieschulze et al., 2009; Mitchell et al., 2012). Nieschulze et al. (2009) used light reflectance in satellite images to derive metrics of canopy structure and successfully correlated these to in situ TF for varying intensities of forest cover in Indonesia. Mitchell et al. (2012) used LiDAR derived canopy height profiles to estimate evapotranspiration, transpiration, interception loss, and forest floor evapotranspiration for a mixed eucalyptus forest in Australia. Similarly, we used LiDAR data to infer the canopy structure for temperate forest ecosystems. In this respect our study represents an addition to existing knowledge for tropical and subtropical forest types and a first attempt to link LiDAR and forest water balance for forests within the temperate climate zone. Although, Nieschulze et al. (2009) and Mitchell et al. (2012) showed a significant potential

J. Schumacher, J.R. Christiansen / Agricultural and Forest Meteorology 203 (2015) 131–141

139

Fig. 5. Modelled throughfall for the forest Hov Skov in southern Zealand for 50 m×50 m grid cells. TF ranged between 39.24 and 71.16 mm, with a mean of 55.98 mm. The standard deviation was 6.43 mm.

to obtain regional predictions of interception, evapotranspiration, and TF using remote sensing, their analyses were confined to a relatively small area and encompassed only mixed forests. Our analyses represented most regions of Denmark and thus captured widely different growing conditions. However, our study comprised short investigation periods of 1–3 years and considered stands in more or less the same development stage. Therefore our model is suited to predict TF for one time period, however, it may not be appropriate for predicting TF into the future. Still we were able to relate canopy structure to TF underlining the robustness of this method. In all our models, monthly, seasonal, and annual, bulk Pr was as expected the dominant factor explaining the variability of TF. However, the LiDAR variable in the annual model contributed to almost 20% of the variation explained by the fixed effects, and 83% of the total variation was explained by this model. This shows that TF can be modeled with a high accuracy on an annual basis, when LiDAR and total Pr data are available. The explanatory power of the

LiDAR variable decreased, when the temporal resolution increased (12.24% and 1.68% for the seasonal and the monthly model, respectively). This might be due to the fact that the short term effects of intensity and duration of the rain events and the wind regime have a pronounced effect on TF on a higher temporal resolution (Keim et al., 2004; Rutter et al., 1971), but their relative importance to total Pr diminishes over time. In our dataset we generally do not capture these dynamic meteorological phenomena, and therefore total precipitation dominates the monthly model. Tair was included in the monthly model, but not in the seasonal or annual models. The negative effect of Tair on TF (Table 5) is expected as higher Tair means higher evaporative loss of intercepted rainfall and hence less TF. However, of all the included variables in our statistical model Tair explained the least of the variance of the model. Our analyses thus show that the contribution of Tair to explain TF on a monthly basis is low compared to the other factors we included. The lack of any pronounced effect of Tair might

140

J. Schumacher, J.R. Christiansen / Agricultural and Forest Meteorology 203 (2015) 131–141

Fig. 6. Evaluation of different age classes related to estimation of throughfall. (a) Mean values of the LiDAR variable (IR) with upper and lower limit for different age classes (note the broken x-axis); (b) Three years (2000–2003) mean throughfall with upper and lower limit for different age classes.

be explained by the fact that evaporative loss from the canopy is dynamic and depends on short term variation of meteorological variables such as wind speed, Tair , humidity and incoming radiation. This is also reflected in the state-of-the-art canopy interception models (Muzylo et al., 2009) that use daily or even hourly meteorological input data to estimate canopy interception evaporation. The monthly winter model performed worse than the monthly summer model (the fixed effects explained 74.76 and 92.98% of the total variation, respectively). This was not expected as the winter model reflects the vegetation state as captured by the LiDAR data. However, the summer model overestimated TF after correcting for back-transformation bias on average by 7.2 mm (mean overestimation of the winter and the combined models was 3.2 and 5.1 mm, respectively). Reason for this could be that the models were fitted using LiDAR data captured during winter (leaf-off conditions). Consequently, in summer (leaf-on conditions) the model overestimates TF. Looking at the parameter estimates the differences between the monthly models are not large (Table 5). The estimated values of the monthly model were not significantly different from the estimated values of the monthly winter and summer models combined (p-value = 0.464). However, comparing the winter estimations of the monthly model with the estimations of the monthly winter model showed significant differences (p-value < 0.001), and comparing the summer estimations of the monthly model with the estimations of the monthly summer model showed significant differences as well (p-value < 0.001). The differences between the estimations of the different models are small (Fig. 3), and the deviation averaged 3.60, 3.94, and 3.27 mm for the combined, winter, and summer comparisons, respectively. We attribute these small differences to model uncertainties rather than a natural phenomenon. Another reason for the worse performance of the monthly winter model could be that ice and snow complicated the measurements which may have resulted in measurement errors and making the monthly winter model less accurate. However, this does not necessarily explain the better performance of the monthly summer model. In this context the question arises if the 3D configuration of the branches, which is what we must capture with the LiDAR in the leaf-off season for deciduous trees, determines the overall orientation and configuration of leaves in the leaf-on season? For future studies it would be interesting to investigate if the branching exerts a dominant control over how the canopy is able to orientate its leaves and thus captures the rain in the leaf-on season. All stands used for modelling, except the mixed species plot at Suserup, the beech plot at Als, and the oak plot at Hald Ege,

had roughly the same age, and therefore it is difficult to make assumptions for other stages of forest stand development, especially before canopy closure which would be relevant in relation to land use change. We assessed the development of TF with stand age and the development of the LiDAR variable with stand age on a different set of plots (the chronosequence) and observed matching trends but no clear linear relationship. Here the LiDAR data was collected in 2006/2007, but the TF and Pr data for these plots were collected in 2000–2003 which probably caused the underestimation of the TF at the Vestskoven-2 (may be due to forest growth) and the overestimation at the Gejlvang site (may be due to opening of forest stands). Future studies should focus on the relationship between LiDAR variables and TF over a longer period of stand development, especially before canopy closure. Broadleaf stands generally had higher TF rates than conifer stands, and the model residuals of broadleaf and conifer stands were not significantly different from each other in all the models. An additional variable for tree type was not selected during model development, and when looking at Fig. 3(d–f) no differences between the two tree types are apparent. This indicates that on this temporal and spatial resolution the models perform well without a pre-stratification into broadleaf and conifer. This also supports our hypothesis that LiDAR based variables capture crown properties and differences in crown structures between tree types. Therefore this approach seems suitable for assessing canopy water flux on larger scales as, for example, water sheds. Our results point to that using simple LiDAR metrics, a realistic estimate of canopy TF can be assessed, given that the total precipitation and, in case of the monthly model mean air temperatures, are known. These findings open up for a new combination of remote sensing analyses and water resource research as the effects of forest types and tree species composition on water balances in a catchment scale can be assessed. Although the application of the model as depicted in Fig. 5 is only a demonstration of how detailed we can map small scale differences in TF, it presents an example of how spatial variability of stand characteristics can potentially affect throughfall and canopy interception within forests. The relevance of this approach must be viewed in conjunction with the fact that land use and climate change will challenge future water resources, and rational management requires improved assessments of regional water balances for forested regions. However, ground surveys of forest water balance are costly, timeconsuming, and cover mostly limited plot sizes, and developing

J. Schumacher, J.R. Christiansen / Agricultural and Forest Meteorology 203 (2015) 131–141

remote sensing techniques can be a more cost-effective way to address these issues. Acknowledgements This study was funded by a scholarship of the University of Copenhagen. We would like to extend our gratitude to Thomas Nord-Larsen, Lars Vesterdal, and Per Gundersen at the Department of Geoscience and Natural Resource Management for valuable advice and access to data. Also, we would like to thank Karin Hansen at the Swedish Environmental Research Institute for providing throughfall and precipitation data for the chronosequences at Vestskoven-2 and Gejlvang. References Carlyle-Moses, D.E., Gash, J.H., 2011. Rainfall interception loss by forest canopies. In: Levia, D.F., Carlyle-Moses, D., Tanaka, T. (Eds.), Forest Hydrology and Biogeochemistry, Volume 216 of Ecological Studies. Springer, Netherlands, pp. 407–423 (Chapter 20). Christiansen, J., Elberling, B., Jansson, P.-E., 2006. Modelling water balance and nitrate leaching in temperate Norway spruce and beech forests located on the same soil type with the CoupModel. Forest Ecol. Manage. 237, 545–556. Christiansen, J.R., Vesterdal, L., Callesen, I., Elberling, B., Schmidt, I.K., Gundersen, P., 2010. Role of six European tree species and land-use legacy for nitrogen and water budgets in forests. Global Change Biol. 16, 2224–2240. COWI, 2007a. Cowi når nye højder. Technical report. COWI, Lyngby, Denmark, 8 pp. COWI, 2007b. Kvalitetsbeskrivelse af laserscanning data. Technical report. COWI, Lyngby, Denmark, 7 pp. Gash, J.H.C., 1979. An analytical model of rainfall interception by forests. Q. J. R. Meteorol. Soc. 105, 43–55. Gregoire, T.G., Lin, Q.F., Boudreau, J., Nelson, R., 2008. Regression estimation following the square-root transformation of the response? Forest Sci. 54 (6), 597–606. Gundersen, P., Sevel, L., Christiansen, J.R., Vesterdal, L., Hansen, K., Bastrup-Birk, A., 2009. Do indicators of nitrogen retention and leaching differ between coniferous and broadleaved forests in Denmark? Forest Ecol. Manage. 258, 1137–1146. Holder, C.D., 2013. Effects of leaf hydrophobicity and water droplet retention on canopy storage capacity. Ecohydrology 6, 483–490. Hosking, J.R.M., 1990. L-moments: analysis and estimation of distributions using linear combinations of order statistics. J. R. Stat. Soc. Ser. B: Methodol. 52, 105–124. Keim, R.F., Skaugset, A.E., Link, T.E., Iroumé, A., 2004. A stochastic model of throughfall for extreme events. Hydrol. Earth Syst. Sci. 8, 23–34.

141

Kristensen, H.L., Gundersen, P., Callesen, I., Reinds, G.J., 2004. Throughfall nitrogen deposition has different impacts on soil solution nitrate concentration in European coniferous and deciduous forests. Ecosystems 7, 180–192. Lefsky, M.A., Cohen, W.B., Parker, G.G., Harding, D.J., 2002. Lidar remote sensing for ecosystem studies. BioScience 52, 19–30. McGaughey, R.J., 2012. FUSION/LDV: Software for LIDAR Data Analysis and Visualization, FUSION Version 3.21. Technical report. United States Department of Agriculture, Forest Service, Pacific Northwest Research Station. Mitchell, P.J., Lane, P.N.J., Benyon, R.G., 2012. Capturing within catchment variation in evapotranspiration from montane forests using LiDAR canopy profiles with measured and modelled fluxes of water. Ecohydrology 5, 708–720. Morsdorf, F., Kötz, B., Meier, E., Itten, K., Allgöwer, B., 2006. Estimation of LAI and fractional cover from small footprint airborne laser scanning data based on gap fraction. Remote Sens. Environ. 104, 50–61. Muzylo, A., Llorens, P., Valente, F., Keizer, J.J., Domingo, F., Gash, J.H.C., 2009. A review of rainfall interception modelling. J. Hydrol. 370, 191–206. Muzylo, A., Valente, F., Domingo, F., Llorens, P., 2012. Modelling rainfall partitioning with sparse Gash and Rutter models in a downy oak stand in leafed and leafless periods. Hydrol. Process. 26, 3161–3173. Næsset, E., 2002. Predicting forest stand characteristics with airborne scanning laser using a practical two-stage procedure and field data. Remote Sens. Environ. 80, 88–99. Nieschulze, J., Erasmi, S., Dietz, J., Hölscher, D., 2009. Satellite-based prediction of rainfall interception by tropical forest stands of a human-dominated landscape in Central Sulawesi, Indonesia. J. Hydrol. 364, 227–235. R Core Team, 2013. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. ˜ D., Valladares, F., Condés, S., Chuvieco, E., 2004. Estimation of leaf area index Riano, and covered ground from airborne laser scanner (Lidar) in two contrasting forests. Agric. Forest Meteorol. 124, 269–275. Rosenqvist, L., Hansen, K., Vesterdal, L., Van Der Salm, C., 2010. Water balance in afforestation chronosequences of common oak and Norway spruce on former arable land in Denmark and southern Sweden. Agric. Forest Meteorol. 150, 196–207. Roth, B.E., Slatton, K.C., Cohen, M.J., 2007. On the potential for high-resolution Lidar to improve rainfall interception estimates in forest ecosystems. Front. Ecol. Environ. 5, 421–428. Rutter, A.J., Kershaw, K.A., Robins, P.C., Morton, A.J., 1971. A predictive model of rainfall interception in forests, 1. Derivation of the model from observations in a plantation of Corsican pine. Agric Meteo 9, 367–384. Rutter, A.J., Morton, A.J., 1977. A predictive model of rainfall interception in forests. III. Sensitivity of the model to stand parameters and meteorological variables. J. Appl. Ecol. 14, 567–588. Scharling, M., 2012. Climate Grid Denmark. Dataset for Use in Research and Education, Technical report. Danish Meteorological Institute. Solberg, S., Brunner, A., Hanssen, K.H., Lange, H., Næsset, E., Rautiainen, M., Stenberg, P., 2009. Mapping LAI in a Norway spruce forest using airborne laser scanning. Remote Sens. Environ. 113, 2317–2327.