The response of canopy height diversity to natural disturbances in two temperate forest landscapes

Structural diversity strongly influences habitat quality and the functioning of forest ecosystems. An important driver of the variation in forest structures are disturbances. As disturbances are increasing in many forest ecosystems around the globe, it is important to understand how structural diversity responds to (changing) disturbances. Our aim was to quantify the relationship between forest disturbances and structural diversity with a focus on diversity in canopy height. We assessed diversity in canopy height for two strictly protected Central European forest landscapes using lidar data. We used a multi-scale framework to quantify within-patch (α), between-patch (β), and overall (γ) diversity. We then analysed the variation in canopy height diversity over an extensive gradient of disturbance rates. Diversity in canopy height was strongly driven by disturbance rate, with highest overall diversity between 0.5 and 1.5% of the forest area disturbed per year. The unimodal responses of overall diversity to disturbance emerged from contrasting within- and between-patch responses, i.e., a decrease in within-patch diversity and an increase in between-patch diversity with increasing disturbance. This relationship was consistent across study landscapes, spatial scales, and diversity indicators. The recent wave of natural disturbances in Central Europe has likely fostered the structural diversity of forest landscapes. However, a further increase in disturbance could result in the crossing of a tipping point (at ~ 1.5% of forest area disturbed per year), leading to substantial structural homogenization.


Introduction
Halting the loss of global biodiversity is one of the current challenges faced by humanity (Steffen et al. 2015). A key contribution of science to addressing this issue is quantifying biodiversity and assessing its relationship to drivers of global change. However, biodiversity is multi-faceted and while some aspects such as species richness have been studied widely (Gaston 2000), other aspects-such as structural diversity (Franklin et al. 2002)-have received less attention to date. Structural diversity, here defined as the variability in the dominant canopy layer, is tightly related to ecosystem functioning. For instance, structural diversity is intrinsically linked to habitat diversity as it provides a variety of niches and more diverse opportunities of exploiting the environmental resources available within a forest (Stein et al. 2014). High structural diversity at the patch scale is also typical for late-successional ecosystems (Franklin et al. 2002), which not only harbour high levels of biodiversity, but are also important contributors to global forest carbon storage (Carey et al. 2001;Zhou et al. 2006;Luyssaert et al. 2008). Structural diversity can also increase the resilience of forest ecosystems (Seidl et al. 2014(Seidl et al. , 2016a, which is increasingly important as both the environment and society change rapidly (Seidl et al. 2016b). As such, understanding the drivers of structural diversity is of crucial importance.
Disturbances are thought to increase structural diversity in forests (Franklin et al. 2002;Turner 2010;Donato et al. 2012), especially if considering diversity at both the patch and landscape scale Kulha et al. 2019Kulha et al. , 2020). Yet, the exact response of structural diversity to disturbances is not well quantified. The intermediate disturbance hypothesis (IDH) has been used to predict the response of biodiversity to disturbances, suggesting that maximum species richness is reached at intermediate disturbance rates (Connell 1978;Wilkinson 1999). While the usefulness of the IDH has been heavily debated in the context of biodiversity (Fox 2013), it might provide a valuable framework for formulating a-priori expectations of how structural diversity can respond to disturbance: Based on the IDH we would expect overall structural diversity to be maximized at intermediate disturbance rates. More specifically, low disturbance rates favour structurally complex (latesuccessional) forests that have high within-patch structural diversity (Franklin et al. 2002;Franklin and Van Pelt 2004). Yet, structural diversity at the landscape scale is low, as large areas will converge towards similar late-successional stages. Conversely, high disturbance rates increase the diversity of structures at the landscape scale (Turner 1989), but at the same time limits the development of structurally complex (late-successional) patches. Consequently, highest structural diversity can be expected to emerge where the trade-off between patch-and landscapescale diversity is maximized, which is the case for intermediate disturbance rates.
A better understanding of the response of structural diversity to disturbance across scales is urgently needed in order to predict how diversity might respond to ongoing changes in natural disturbance regimes . For instance, it remains unclear if an increase in disturbance rates beyond current levels will increase or decrease structural diversity in forests. Better understanding the response of structural diversity to disturbance, however, requires large scale approaches quantifying structural diversity from the patch to landscape scales, which still remains challenging. Here, the emergence of active remote sensing systems, such as Light Detection and Ranging (Lidar), can make an important contribution (Mura et al. 2015(Mura et al. , 2016Valbuena et al. 2016). Lidar allows for the assessment of canopy structure across large spatial extents, delivering insights into the vertical distribution of the canopy layer. As such, Lidar offers great potential for testing the multi-scale disturbancestructural diversity relationship described above.
We here aim at quantifying the response of structural diversity to disturbance. We focus on one important dimension of structural diversity, the diversity in canopy height classes (Kuuluvainen et al. 1996), by answering the following three questions: (1) What is the effect of natural disturbance rates on canopy height diversity? (2) What is the contribution of patch-scale diversity (i.e., variability of canopy height classes within the patch) versus landscape-scale diversity (i.e., variability of canopy height classes between patches) on overall diversity in canopy height classes? (3) At which disturbance rates is diversity in canopy height maximized? We quantify canopy height diversity using Lidar data and a multi-scale framework that partitions canopy height diversity into three components describing within-patch (i.e., patchscale), between-patch (landscape-scale) and overall diversity. The analysis is performed for two strictly protected Central European forest landscapes that bracket the natural disturbance regimes of temperate forests, ranging from from small-to large-scale disturbances .

Study landscapes
Most forests within Central Europe are managed by humans, imposing a strong management signals upon the relationship between structural diversity and natural disturbances. To control for management effects in our study, we focused on two strictly protected forest landscapes where management was ceased more than 40 years ago: The Bavarian forest national park (founded in 1970) and Berchtesgaden national park (founded in 1977). Both landscapes are situated in Germany, bordering Czechia and Austria, respectively ( Fig. 1), and are representative for a large variety of Central European forest types from the submontane to the sub-alpine elevation belt. Dominant tree species in both study landscapes are European beech (Fagus sylvatica L.) and Norway spruce (Picea abies H. Karst), with silver fir (Abies alba (L.) Mill.), sycamore maple (Acer psuedoplantanus L.), and mountain-ash (Sorbus aucuparia L.) as other important canopy species. At Berchtesgaden, subalpine forests are characterized by European larch (Larix decidua Mill.) and Swiss stone pine (Pinus cembra L.), forming the upper treeline at about 1800 m above sea level. The landscapes strongly differ in their topographic template: The Bavarian forest landscape is a typical highland landscape (Mittelgebirge), with only moderate elevation gradients and only a few peaks above 1250 m. The Berchtesgaden landscape, in contrast, is a typical high mountain landscape with strong elevation gradients and several peaks extending to more than 2000 m. These differences result in contrasting land cover patterns in both landscapes ( Fig. 1): While forest cover is continuous and fragmentation is low in the Bavarian forest landscape, the Berchtesgaden landscapes is characterized by highly fragmented forests of mostly small to moderately sized patches that are often isolated.
The major natural disturbance agents in both landscapes are bark beetle (mostly Ips typographus L.) and windthrow, with snow avalanches being an additionally important disturbance agent in Berchtesgaden. While both landscapes have experienced high disturbance activity in recent decades (Senf et al. 2017), the spatial patterns of disturbance were highly divergent. Disturbances in Berchtesgaden are characterized by patchy, small-scale disturbance activity, whereas the Bavarian forest features some of the largest natural disturbance patches recorded in Central Europe in recent history .

Data
We acquired forest cover and annual forest disturbance maps from a previous study mapping both forest cover and forest disturbances from Landsat satellite data (Senf et al. 2017). Maps were created annually at 30 m spatial grain and cover the time period 1984 to 2016. They depict all pixels that have experienced a substantial change in canopy cover larger than 0.5 ha regardless of agent, with overall accuracies ranging from 81 to 87%. Our analysis thus focusses on standreplacing disturbances rather than gap-dynamics.
We used canopy height models developed from Lidar to quantify canopy height diversity across scales. The canopy height model for the Bavarian Forest landscape was generated by national park authorities from Lidar point clouds with an average point density of 35 returns per m 2 and represents the state of the landscape in the year 2012. The canopy height model for Berchtesgaden was generated by the Bavarian State Office for Survey and Geoinformation from Lidar point clouds with a minimum point density of 10 returns per m 2 , based on data acquired in 2009. Canopy height models were generated with a horizontal grain of 1 m and classified into 1 m vertical bins representing discrete canopy height classes of 1 m step width. We finally focused our analysis on the tree layer (i.e., canopy heights of C 2 m), removing pixels dominated by shrubs, herbs, or tree seedlings (Latifi et al. 2016).

Sampling design
Our primary unit of analysis was a 30 9 30 m grid cell, which is henceforth referred to as the patch scale. We note, however, that disturbances patches can be composed of several pixels in reality, and the purpose of the sampling unit here was simply to calculate the average canopy height diversity within a standardized patch of 30 9 30 m. To cover a wide range of disturbance and environmental conditions we choose a sampling scheme that randomly drew 150 sublandscapes of N ¼ n Â n patches within each study landscape (see Fig. 2). The sub-landscape extent varied across n ¼ f10; 15; 20; 25; 30; 35g in order to account for effects of scale explicitly. For each sublandscape we calculated the proportion of forest cover and removed sub-landscapes with less than 40% of the area covered by forests. This threshold corresponds to common forest definitions (Chazdon et al. 2016) and excludes areas predominantly characterized by nonforest land cover (e.g., rocks, lakes, etc.).

Quantifying canopy height diversity
To accommodate the multi-scale nature of structural diversity, we calculated alpha, beta and gamma diversity in Lidar-based height classes for each sublandscape following the approach of Jost (2007) and summarized in Table 1. Alpha diversity quantifies the average canopy height class diversity within all patches of a sub-landscape and is subsequently referred to as within-patch diversity. Beta diversity characterizes the difference in height class compositions among the patches of a sub-landscape (canopy height turnover) and thus describes the between-patch diversity. Gamma diversity is the pooled canopy height diversity within the entire sub-landscape and is henceforth referred to as overall diversity. We calculated within-patch, between-patch and overall diversity for three measures of diversity, that is canopy class richness ( 0 D), Shannon entropy ( 1 D) and the Gini-Simpson index ( 2 D), respectively (Table 1). These indicators represent a sequence of increasing importance of rare cases (here: rare canopy height classes) in the quantification of diversity (Hill 1973). All calculations were done using the vegetarian package in the R software for statistical computing (R Core Team 2019).

Disturbance metrics and covariates
We calculated the average annual disturbance rate (percent of forest area disturbed per year) as measure of disturbance for each sub-landscape. We did so by using the Landsat-based disturbance and forest cover maps described above. While other descriptors of disturbance activity could also be used, we here note that disturbance rate is generally well correlated with important spatial disturbance characteristics, such as patch density and patch aggregation (Supplementary Materials Fig. S6). Besides disturbance rate, we also included each sub-landscape's forest proportion as well as a measure of each sub-landscape's topographic complexity (coefficient of variation in elevation) as covariates in our model. We did so to account for the highly contrasting topographic templates and landscape compositions found across our study landscapes (Table 2). We further included the sub-landscape size (n; see Fig. 2) as random effect in the model, allowing the intercept to vary randomly among sub-landscape of variable sizes.

Statistical analyses
We used Bayesian Generalized Additive Mixed Models (GAMMs) to regress within-patch, betweenpatch and overall canopy height diversity over the disturbance measures and covariates described above. As we expected a non-linear relationship between disturbance rates and canopy height diversity, we preferred GAMMs over linear models. In particular, GAMMs allow for modelling potentially complex, non-linear relationships by applying smoothed functions to the regression problem (Wood 2017). We used thin plate regression splines with 10 basis functions as implemented in the mgcv package (Wood 2003). Statistical inference was done via Bayes theorem using Monte-Carlo Markov Chain methods implemented in the Stan language for probabilistic modelling (Carpenter et al. 2017). We used weaklyinformative, regularizing priors for all parametersincluding the degree of smoothing-which prevents over-fitting while simultaneously allowing for some degree of surprise when there is strong evidence in the data (McElreath 2018).
To test whether a unique relationship between disturbance rates and canopy height diversity existed for each study landscapes, or whether a pooled model could be generalized across both study landscapes, we fitted one model with, and one model without the study landscape as modulating covariate. We then compared the predictive performance of both models using a leave-one-out estimate of the expected log predictive density. The leave-one-out estimated log predictive density is a measure of relative model performance when predicting new, unknown data points (Gelman et al. 2014;Vehtari et al. 2016), and can be interpreted-like similar information criteria-as a relative measure of predictive performance when comparing nested models.
Once the final model was chosen (i.e., with or without landscape as covariate), we drew posterior Fig. 2 Example of a sub-landscape (b) of the Bavarian forest landscape (a), composed of N = n 9 n patches of 30 9 30 m (c). Each patch contains canopy height information at a horizontal grain of 1 m predictions from the model over the range of disturbance rates observed in our study landscapes, controlling for differences in forest proportion and topographic complexity by holding them constant at their overall mean (see Table 2). We thus report the relationship between disturbances and canopy height diversity for standardized sub-landscapes of average forest cover proportion and topographic complexity.

Results
Disturbance rate was a strong driver of canopy height diversity across all scales and diversity indices (Table 3; for posterior model checks see Supplementary Materials S1 and S2). There was only weak evidence for an improvement of predictive capacity when including the study landscape as modulating covariate (see Supplementary Materials S3), and we thus used one pooled model for making predictions for both landscapes in the following. However, there was evidence for a difference in maximum canopy height diversity among landscapes ( Fig. 3 and Supplementary Materials S3).
In general, within-patch canopy height diversity decreased with increasing disturbance rate, with maximum values at low disturbance rates and a moderate decrease in the undisturbed systems (Fig. 3). Between-patch diversity responded unimodally to disturbance rate for 1 D (Shannon) and 2 D (Simpsons) diversity (Fig. 3). Thus, the turnover of canopy structure was highest at intermediate disturbance rates. For richness in canopy classes ( 0 D), however, we found an exponentially increasing pattern, indicating that the number of canopy classes on the landscape continuously increased within the range of observed disturbance rates. Overall canopy height diversity, emerging from the combined patterns of within-and between-patch diversity, was always maximal at intermediate disturbance rates and approached minimum values at very high disturbance activity (Fig. 3). With respect to landscape extent (n), we did not find substantial differences between extents ( Supplementary Materials Fig. S4). Hence, the relationships identified here are consistent over different landscape extents. Within-patch diversity (a-diversity)

Between-patch diversity (b-diversity)
Overall diversity (c-diversity) Across all equations, N indicates the total number of patches within a sub-landscape (see Fig. 2), S indicates the total number of distinct height classes, and p i is the proportional abundance of each height class within each patch

Discussion
Natural disturbances are increasing across global forests as a consequence of climate change , which makes quantifying disturbance impacts on biological diversity a key priority of research . We here show that the response of canopy height diversity to changing disturbance rates is unimodal and consistent across indicators of diversity and landscapes with widely varying topography. Our analyses thus provide empirical evidence for the general notion that disturbances foster structural diversity, in our case canopy height diversity. For disturbance rates between 0.5 and 1.5%, increasing between-patch diversity compensates for decreasing within-patch diversity, resulting in relatively stable overall structural diversity for those intermediate disturbance levels. At disturbance rates of [ 1.5% per year, however, overall canopy height diversity decreases rapidly, with increasing disturbance rates resulting in a structural homogenization of forests at both the patch and landscape scale. Between 1986 and 2016, canopy disturbance rates in unmanaged forests of Central Europe ranged between 0.3 and 1.0% per year (Senf et al. 2017) and were thus well within the ranges of maximum canopy height diversity identified herein. Recent natural disturbances are thus unlikely to threaten the structural diversity of Central Europe's temperate forests. However, canopy disturbance rates in the region have doubled in the last three decades in response to changes in climate and management regimes . A further increase along the current trajectory could thus be expected to result in a structural homogenization from disturbances in the coming decades.
Based on the IDH we expected that canopy height diversity would be maximized at intermediate disturbance rates and our results strongly support those theoretical predictions. Yet, the IDH has been criticized based on theoretical grounds, calling into question the processes thought to underline it (i.e., weakened competition, interrupted competitive exclusion, and time lags in the response of species to a fluctuating environment; Fox 2013). We note that all these arguments pertain to the diversity of species, while we here focus on the diversity of structures in an ecosystem, a dimension of diversity that has only rarely been investigated in the context of the IDH. The IDH might-based on our results-indeed be well suited for predicting the response of forest vertical structure to disturbances. We also show that different spatial scales of analyses strongly influenced the Disturbance rate (% yr -1 ) 1.8 ± 1.1 (0.0-3.6) 0.2 ± 0.4 (0.0-3.6) 1.01 ± 1.20 (0.0-3.6) Forest cover proportion (%) 89.6 ± 12.2 (40.0-100.0) 72.4 ± 19.0 (40.0-100.0) 81.9 ± 18.6 (40.0-100.00) Topographic complexity (%) 2.9 ± 1.6 (0.2-9.0) 9.0 ± 5.8 (0.9-28.5) 6.40 ± 5.64 (0.2-28.5) Shown are mean values ± standard deviation, with minimum and maximum values in brackets. Topographic complexity is here expressed as the coefficient of variation of elevation relationship between disturbances and canopy height diversity, underlining the importance of considering both the patch and landscape scale for assessing disturbance effects on biodiversity Kulha et al. 2019Kulha et al. , 2020. For instance, while the number and abundance of canopy height classes decreases rapidly with increasing disturbance rates at the patch scale (i.e., patch-scale homogenization through stand-replacing disturbances; (Franklin et al. 2002)), the turnover of canopy classes at the landscape scale increases (higher variability of successional stages on the landscape; Turner 2010). We thus show that a unimodal pattern consistent with the prediction of the IDH emerges from diverging patterns at the patch and landscape scale, while different patterns would be observed if just focusing on one of each scale. We found that the response of canopy height diversity to disturbances was remarkably consistent across different indicators of diversity. This is particularly noteworthy as the three indicators investigated here differ strongly in the way they consider the occurrence of rare canopy height classes, with increasing weights of rare classes from 0 D to 2 D. The only contrasting pattern found was for between-patch canopy height diversity indicators of canopy class richness ( 0 D) and diversity ( 1 D, 2 D). The divergent trend of the lower order diversity index ( 0 D) can be explained by the fact that even under very high disturbance rates there might be isolated islands of Fig. 3 Effects of disturbances on canopy height diversity. The different colours indicate the two study landscapes, but we note that model predictions were derived from a single, pooled model generalizing across both landscapes. The y-axis shows canopy height diversity among the three indices (richness, Shannon, Gini-Simpson; see Table 1) and scales (within-patch, betweenpatch and overall; see Table 1). Individual relationships for each sub-landscape extent are shown in Supplement Materials Figure S4 and raw data is shown in Supplement Materials Figure S5 surviving trees, increasing the richness of canopy height classes while the abundance will be reduced substantially.
Despite the advances in understanding the response of structural diversity to disturbance presented herein, further research is needed to more comprehensively understand the role disturbances play in structuring forest ecosystems. We here, for instance, neglect the spatial variability of trees at the sub-patch scale, which often is high in naturally disturbed forests (Swanson et al. 2011;Winter et al. 2015). Furthermore, while we here aimed to assess structural diversity with a more comprehensive set of indicators, several important aspects of forest structure, such as variation in tree diameters, are not accounted for. Also, our indicators of canopy structure focuses exclusively on live trees, excluding important structural legacies of disturbances such as snags and downed deadwood, which considerably contribute to patch-scale biodiversity after disturbance (Swanson et al. 2011;Donato et al. 2012). We furthermore use a very simple descriptor of disturbance activity (i.e., the average disturbance rate), which might not distinguish between frequent small disturbances and rare large disturbances. Likewise, the satellite-based disturbance maps used herein only capture stand-replacing disturbances and we thus miss canopy gaps resulting from small-scale gapdynamics important in many old-growth forest ecosystems (Kulha et al. 2019). Finally, while our study landscapes were selected to bracket recent natural disturbance activity in Central Europe , they both were affected primarily by disturbances from wind and bark beetles. As wildfire is also an important disturbance agent in temperate forest ecosystems (Sommerfeld et al. 2018), future work should test the generality of our results to other disturbance regimes and ecosystems.
Our findings are highly relevant for ecosystem managers, both in protected areas and commercial forests. We here highlight that naturally developing forests in National Parks are generally well adapted to their disturbance regimes, and that natural disturbances increase the diversity of structures in forests across a wide range of disturbance rates. Natural disturbances are thus fostering conservation goals. Management for resource extraction often results in substantially increased disturbance levels compared to natural disturbance regimes (Sommerfeld et al. 2018). In managed forests, structural diversity might thus be decreased compared to their un-managed counterparts. In this regard, our findings suggest that newer forest management approaches such as retention forestry (Gustafsson et al. 2012;Mori and Kitagawa 2014) and landscape-scale management of forests (Schall et al. 2018;Seidl et al. 2018) could help to retain structural diversity, e.g., by retaining trees during harvesting operations and applying landscape management schemes that vary silvicultural prescriptions between stands.

Conclusions
We here empirically show that canopy height diversity is strongly driven by natural disturbances in temperate forest ecosystems, with disturbances generally fostering diversity in canopy height classes. Yet, we also show that the relationship between disturbance rates and canopy height diversity is non-linear and that increasing disturbance rates can substantially reduce structural diversity once pushed beyond levels of historic variability. We further highlight the importance of considering structural diversity at multiple spatial scales and underline the power of increasingly available Lidar-based canopy height models-as well as novel space-borne Lidar data (Dubayah et al. 2020)-for such assessments. The novel assessment framework developed herein can be easily applied to other regions, allowing for a consistent estimation of the effects of disturbances on canopy height diversity across a range of disturbance regimes. Finally, we highlight that a continued increase in canopy disturbance rates as observed in the past three decades could substantially decrease forest structural diversity. Such a decrease might have negative impacts on important ecosystem functions, such as habitat quality and forest carbon storage, and could potentially erode the resilience of Europe's forests. medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Author contributions CS and RS conceived the ideas and designed methodology; JM provided the data; CS analysed the data; CS, RS and AM led the writing of the manuscript. All authors contributed critically to the drafts and gave final approval for publication.
Data and code accessibility Data and code to reproduce the results of this article are available under following link: https:// doi.org/10.5281/zenodo.3904260.