\documentclass[fleqn,10pt]{wlscirep}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}

%% The amssymb package provides various useful mathematical symbols
\usepackage{amssymb}
%% The amsmath package provides various useful equation environments.
\usepackage{amsmath}
\usepackage{gensymb}
%% The amsthm package provides extended theorem environments
%% \usepackage{amsthm}

\usepackage{hyperref}
\usepackage{float}
\usepackage{tabularx}
%\usepackage{cancel}
\usepackage{multirow}
\usepackage{array}

%% track changes
\usepackage[dvipsnames,svgnames,x11names]{xcolor}
\usepackage[marginparwidth=2.5cm]{geometry}
\usepackage[final]{changes} % to show changes: [markup=default]

%%% Alternative definition to have the remarks
%%% in the margins instead of footnotes
\usepackage{todonotes}
\setcommentmarkup{\todo[color={authorcolor!20},size=\scriptsize]{#3: #1}}


%% Rather hacky definition of an "annote"
%% by riding on \added
\newcommand{\note}[2][]{\added[#1,comment={#2}]{}}

\usepackage{natbib}
\usepackage{lineno}
%% The lineno packages adds line numbers. Start line numbering with
%% \begin{linenumbers}, end it with \end{linenumbers}. Or switch it on
\linenumbers

\begin{document}

\title{Recovery of the forest's protective effect after stand-replacing wind disturbances}

\author[1,*]{Christine Moos}
\author[1]{Kaya Dietrich}
\author[1]{Alexandra Erbach}
\author[2]{Christian Ginzler}
\author[1]{Estelle Noyer}
\author[1]{Christoph Schaller}
\author[1]{Luuk Dorren}
\affil[1]{Bern University of Applied Sciences, BFH-HAFL, Zollikofen, CH-3052, Switzerland}
\affil[2]{Swiss Federal Institute for Forest, Snow and Landscape Research WSL, Birmensdorf, CH-8903, Switzerland}

\affil[*]{christine.moos@bfh.ch}


%\keywords{Protective forest, disturbance, recovery, natural hazards}

\begin{abstract}
Forests in mountain areas provide an indispensable ecosystem service by protecting people and infrastructure from natural hazards. However, open questions remain regarding its long-term and reliable provision, which will amplify with climate change. In particular, natural disturbances, such as wind storms, can cause temporary losses of the forest's protective effect \added{against natural hazards, such as rockfall, snow avalanches or landslides}, potentially leading to higher \deleted{natural hazard} risks \added{for settlements and infrastructures}. In this study, we analyzed the recovery of the protective effect of a large number of windthrow areas in the Swiss Alps using high-resolution canopy height models and single tree detection. We observed a considerable recovery of tree height (67 \% on average), forest cover (90 \%) and stem density (46 \%) 31 years after the storm. Basal area showed a lower recovery of 16 \% on average with 25 \% of the studied areas showing a retarded recovery of $\leq$ 6 \%. This highlights the extended timescales required for full recovery of the protective effect against certain hazard types, such as rockfall or landslides. Recovery rates were most importantly influenced by temperature, precipitation, local topography, gap size and lithological substrate, underlining the importance of local environmental conditions when predicting recovery. The study highlights future research gaps, including the role of deadwood as well as species mixture in the recovery of the forest's protective effect after disturbances.
\end{abstract}


\flushbottom
\maketitle

\thispagestyle{empty}


\section*{Introduction}

Throughout the Alps, human settlements have benefited from the widespread forest cover on steep slopes, providing a natural means of protection against natural hazards\citep{motta2000protective}. These protective forests reinforce slopes with their root systems \citep[e.g.,][]{schwarz2010quantifying}, stabilize the snow cover \citep[e.g.,][]{casteller2018assessing}, reduce the energy of falling \deleted{blocks} \added{rocks} \citep[e.g.,][]{dorren2005mechanisms}, and intercept precipitation, improve water infiltration and storage in the soil and can thus regulate runoff \citep[e.g.,][]{bathurst2020forests, shah2020review}. In doing so, forests can reduce the frequency and intensity of gravitational hazards, such as landslides, snow avalanches, rockfall and torrential floods \citep{moos2018ecosystem}. Without this protective service of forests, large engineered structures would be necessary, entailing immense costs for installation and maintenance \citep{getzner2017gravitational, accastello2019asforesee, poratelli2020integrated}. \\
Although the importance of the protection provided by forests in mountains is generally undisputed and has also been proven by numerous scientific studies \citep[see e.g.,][]{sudmeier2021scientific, moos2018ecosystem}, open questions regarding its long-term and reliable provision remain \citep{moos2023mountain}. In particular, natural disturbances, i.e., tree mortality due to forest fires, insect outbreaks or wind storms, can cause temporary losses of the forest's protective effect, potentially leading to higher natural hazard risks \citep[e.g.,][]{vacchiano2016assessing, sebald2019effects}. With ongoing climate change, the frequency and severity of these disturbances are likely to increase \citep{seidl2017forest, patacca2023significant} and may compromise the reliability of protective forests as natural means of protection. To reduce the uncertainties related to the impact of natural disturbances on the forest's protective effect, not only answering the question to what extent the protective effect is reduced is pivotal, but also how fast it recovers after the disturbance event. Despite the considerable interest in the quantification of forest recovery after disturbances \citep[e.g.,][]{cerioni2024recovery, seidl2022post}, we currently lack large-scale data specifically on protective forests to answer this question in a satisfactory and generalizable manner.
Available evidence on the recovery of the protective effect after disturbances largely stems from field-based studies. In Switzerland, for example, great efforts have been made to monitor the regeneration as well as deadwood decay after the two severe storm events\added{:} Vivian in 1990 and Lothar in \added{1999} \deleted{2000} \citep[e.g.,][]{kramer2014site, wohlgemuth2017post, piazza2024post}. Although these studies provide important insights on the forest's recovery after windthrow disturbances and its potential consequences for the protection against natural hazards, they mostly focus on a limited number of case study sites, thus impeding generalizable conclusions. \\
On the other hand, remote sensing data has increasingly been used to investigate forest recovery after disturbances \citep[e.g.,][]{bartels2016trends, white2020quantifying, white2022mapping}. Most of these studies use tree cover and/or height derived from satellite images as measure of forest recovery \citep[e.g.,][]{senf2019post, stritih2024mapping}. However, a closed crown cover has limited significance with regard to the protective effect. To assess the latter, more detailed forest structural properties, such as tree density and tree diameter distribution, are necessary \citep{dupire2016novel}. High-resolution canopy height models (CHM) may thus be a promising alternative to extract more detailed structural variables, while being increasingly available on a large scale and with high temporal and spatial resolution \citep{schaller2023improving}. These CHM's can be derived using data collected by active sensors such as Light Detecting and Ranging (LiDAR) or aerial photogrammetry \citep{liu2020analysis}.\\
In this study, we used nationwide \added{high-resolution} CHM's \added{based on airborne LiDAR data} to assess the recovery of the protective effect after the storm Vivian in 1990 \added{for two time steps}. We quantified recovery \added{13 and 31 years} \deleted{for two time steps} after the storm \deleted{(13 and 31 years)} based on different forest structural variables for more than 1700 \added{mapped} windthrow areas \added{with severe storm damages} in protective forests in 5 cantons in Switzerland. We thereby determined the recovery of the protective effect of the windthrow areas by comparing them to undisturbed \added{adjacent} reference areas, and analyzed it as a function of environmental conditions. The objective of this study was to answer the following two research questions:
(i) How well did the protective effect, in terms of forest cover, tree height, tree density and stand basal area, of the windthrow areas recover during the following 13 and 31 years?
(ii) Which are the main \added{topographic and climatic} drivers influencing the recovery rate of the protective effect?

\section*{Results}
\subsection*{Recovery of protective effect}
A clear increase in forest cover\added{(i.e. ratio of pixels $\geq$ 3 m and pixels $<$ 3m)}, vegetation height, tree density, tree diameter and basal area was found for the windthrow areas 13 years (in 2003) and 31 years after the storm (in 2021), respectively (Figure \ref{fig:distVar}). While the median basal area of the windthrow areas was 0.52 m\textsuperscript{2}*ha\textsuperscript{-1} after 13 years, it increased to 4.4 m\textsuperscript{2}*ha\textsuperscript{-1} in 2021. Median vegetation height in 2021 was 7.6 m (0.4 m in 2003) and forest cover 75 \% (30 \% in 2003). Tree density reached 112 stems/ha in 2021 (30 stems/ha in 2003) and the median diameter at breast height (dbh) 20.7 cm (13.7 cm in 2003). 
Overall, the recovery of vegetation height (determined as the ratio between the windhtrow and the corresponding reference areas) was on average 67 \% 31 years after the storm (26 \% in 2003; Figure \ref{fig:Recov}), whereas 20 \% of the areas reached a height recovery of more than 85 \%. Forest cover recovered by 50 \% after 13 and by 90 \% after 31 years on average. Tree density recovered by 46 \% after 31 years (18 \% in 2003), while basal area reached a mean recovery of 16 \% after 31 years (3.4 \% in 2003). 48 out of the total of 1740 windthrow areas showed a basal area recovery $\leq$ 1 \% after 31 years and about one-fourth of the areas (n=433) reached a recovery of maximum 6.6 \%.
The properties of the reference areas remained relatively constant between 2003 and 2021. Yet an increase in forest cover (median = 62 \% in 2003 and 81 \% in 2021) was observed. On the other hand, median dbh (38.5 cm in 2003, 40.1 cm in 2021), tree density (141 stems/ha in 2003, 149 stems/ha in 2021) and vegetation height (median = 15.0 m in 2003 and 16.0 m in 2021) changed only marginally.

\begin{figure}[H]
\centering
\includegraphics[width=0.6\textwidth]{figures/Variables_ref_storm_ridge_new.png}
\caption{Distribution of median vegetation height (a), stem number (b), median tree diameter (c) and basal area (d) in windthrow and reference areas in 2003 (13 years after storm Vivian) and 2021 (31 years after Vivian).} \label{fig:distVar}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics[width=0.8\textwidth]{figures/Fig2_adapted.png}
%% Use \caption command for figure caption and label.
\caption{Distribution of the recovery (determined as ratio between windthrow and undisturbed reference areas) of forest cover (a), vegetation height (b), stem density (c) and basal area (d) 13 \added{2003} and 31 \added{2021} years after the storm Vivian.}\label{fig:Recov}
\end{figure}


\subsection*{Variables influencing recovery}

We found a clear trend of a decreasing recovery of the protective effect (in terms of basal area recovery) with increasing elevation and decreasing mean seasonal temperature, day degree sum \added{(i.e., mean yearly sum of daily temperature $>5 \degree$)} and minimum winter temperature (Figure \ref{fig:Recov_elev}). Additionally, recovery was significantly higher in larger windthrow areas (p \deleted{value} $\leq$ 0.05). Slope inclination had a significantly negative effect in 2021, but not in 2003. Its influence tended to decrease with increasing altitude. Furthermore, there was a significantly positive correlation with planar curvature with highest recovery on linear slopes. \deleted{Easterness and northerness did not have a significant effect in the univariate analysis. However, a significant negative effect was detected in the multivariate generalized linear model (GLM, Table \ref{tab:glm})} \added{Northerness significantly influenced basal area recovery in the initial generalized linear model (GLM), implying a decreasing recovery on more northern exposed slopes. This effect, however, was not significant any more after including a spatial correlation structure in the model fitted with generalized least squares (GLS) (Table \ref{tab:glm}).} Furthermore, recovery tended to be smaller on limestone \added{and conglomerate} slopes compared to slopes of silicious rocks\deleted{, conglomerate} or blocky scree.
Overall yearly precipitation did not significantly influence basal area recovery, but had a positive effect on stem density recovery in 2021. On the other hand, both basal area and stem density significantly increased with median spring precipitation in 2021. Additionally, minimum summer precipitation was positively correlated with the recovery of stem density, but not basal area. Basal area recovery on slopes above 1000 m significantly increased with medium and maximum drought stress in 2003, while in 2021 only on slopes above 1600 m (Figure \ref{fig:Recov_Dr}). In 2021, stem density recovery decreased with increasing maximum drought stress at all elevations as well as basal area recovery on limestone sites. The maximum drought stress between 2003 and 2021 was significantly higher compared to the previous period between 1990 and 2003.\\
In the final model of basal area recovery, the explanatory variables that remained were elevation, slope, \deleted{northerness,} lithology, mean spring precipitation, drought stress, time since the disturbance as well as gap size (Table \ref{tab:glm}). Furthermore, there was a significant effect of the interaction of the maximum drought stress with the time since the disturbance (p $\leq$ 0.05) \deleted{and an almost significant effect of the interaction of northerness and elevation (p $\leq$ 0.1)}. The final model had a \added{pseudo} R\textsuperscript{2} of 32 \% and a RMSE (root mean squared error) of 9.6 \%.

\begin{figure}[H]
\centering
\includegraphics[width=0.8\textwidth]{figures/BA_elev_DDsum.png}
\caption{Basal area recovery (a)  in 2003 (13 years after Vivian) and 2021 (31 years after Vivian) and median day degree sum (b) for the periods 1990-2003 and 2004-2021 for three elevation zones.}\label{fig:Recov_elev}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics[width=0.8\textwidth]{figures/Recovery_N_BA_DrStr_elev.png}
%% Use \caption command for figure caption and label.
\caption{Recovery of stem density (top) and basal area (bottom) as a function of maximum drought stress in 2003 (13 years after Vivian) and 2021 (31 years after Vivian).}\label{fig:Recov_Dr}
\end{figure}

\begin{table}[H]
	\addtolength{\tabcolsep}{-4pt}
	\fontsize{9.5pt}{10.25pt}\selectfont
	\renewcommand*{\arraystretch}{1.4}
    \caption{Coefficients, standard errors and significance of the intercept and the predictor variables of the final model of basal area recovery \added{fitted with generalized least squares (GLS) including a spatial correlation structure}. Significance level for p-value $<$ 0.05 = *, p $<$ 0.01 = **, p $<$ 0.001 = ***, 0.05 $<$ p $<$ 0.1 = (*); \textit{mean spring prec.}: mean of the precipitation in the months march, april, mai over the observation periods (1990-2003, 2004-2021); \textit{drought stress}: maximum drought stress over the observation periods; \textit{years since}: years since the storm Vivian.}
	\begin{tabular}{@{}llll@{}}
		Predictor  & Coefficient & Standard error & significance \\
        \hline \\
		intercept & -4.3*10\textsuperscript{-1} & 3.6*10\textsuperscript{-1} &  \\
		elevation & -1.7*10\textsuperscript{-3}  & 2.0*10\textsuperscript{-4} & ***\\
	    slope & -2.8*10\textsuperscript{-2}   & 1.0*10\textsuperscript{-2} & ** \\
		mean spring prec. & 1.6*10\textsuperscript{-3}  & 4.4*10\textsuperscript{-4} & *** \\
		conglomerate  & -3.0*10\textsuperscript{-1} & 1.4*10\textsuperscript{-1}  &  * \\ 
        siliceous  & -1.4*10\textsuperscript{-1} & 1.1*10\textsuperscript{-1}  &    \\
        limestone & -3.9*10\textsuperscript{-1}  & 1.1*10\textsuperscript{-1}  &  ***   \\
        sandstone  & -3.1*10\textsuperscript{-1}  & 1.1*10\textsuperscript{-1}  &  **  \\
        log(gap size)  & 2.7*10\textsuperscript{-1} & 2.3110\textsuperscript{-2}  & ***   \\
        plan curvature  & 8.0*10\textsuperscript{-2} & 1.7*10\textsuperscript{-2}  &  ***   \\
        drought stress & 4.9*10\textsuperscript{-1} & 1.1*10\textsuperscript{-1}  & ***   \\
        years since & 2.4*10\textsuperscript{-1} & 6.0*10\textsuperscript{-3}  & ***   \\
        elev:slope & 2.2*10\textsuperscript{-5} & 7.4*10\textsuperscript{-6} & ** \\
        drought stress:years since & -1.8*10\textsuperscript{-2} & 6.2*10\textsuperscript{-3}  & ** \\
        \hline
		\end{tabular}
		\vspace{4\baselineskip}
		\label{tab:glm}
	\end{table}

\section*{Discussion}
The results of this study indicate that protective forests in the Swiss Alps damaged by the wind storm Vivian in 1990 have considerably recovered after 30 years. We observed a significant increase in the recovery of tree height, forest cover, stem density and basal area between 13 and 31 years after the storm. However, the degree of recovery strongly varied between the considered forest structural variables. While tree height recovered on average by about 70 \% and forest cover by about 90 \% after 31 years, the recovery of stem density of trees $\geq$ 3 m (40 \%) and basal area (15 \%) were critically lower. This has different implications for the potential natural hazard processes a specific forest stand protects from: While a largely closed canopy with trees reaching up to two-thirds of the initial tree heights may be sufficient for the prevention of avalanche releases \citep{frehner2005nachhaltigkeit}, basal area is a more meaningful proxy of the recovery of the protective effect against rockfall \citep[e.g.,][]{moos2017quantifying} or shallow landslides \citep[e.g.,][]{cislaghi2021effects, mao2015modelling}. Thus, the results of this study indicate a substantially longer recovery of the protective effect against rockfall or landslides compared to snow avalanches. Similarly, \cite{may2023quantifying} showed that spruce, fir and beech forests regain the maximum protective effect against rockfall after 50 to 200 years based on an analysis of Swiss national Forest Inventory data. Another study reported an optimal protection against small rocks (volume = 0.05 m\textsuperscript{3}) 50 years after the disturbance \citep{vacchiano2015effect}. 
\cite{stritih2024mapping} reported an average probability of 61 \% of protective forests across the Alps for having a sufficient protective forest structure in terms of canopy cover and tree height 30 years after a wind disturbance. 

Our study underlines that forest cover and tree height might not be sufficient as a general proxy for the recovery of the protective effect, as this is often done in large-scale studies \citep[e.g.,][]{stritih2024mapping, vacchiano2016assessing}. Besides a more detailed characterization of forest structure, the hazard disposition (e.g., slope inclination, snow height, rock volumes, slope length) should ideally also be taken into account to determine the recovery of the protective effect for the specific local conditions. \added{Thus, it is important to analyze the impact of a disturbed area on the forest's protective effect in relation to the total slope and the hazard disposition in local scale studies. By way of example, a slow recovery of basal area in a windthrow area in a rockfall protective forest might not be problematic in case the disturbed area is small compared to the total forested slope length.} \\
Although our results indicate overall a good resilience of the areas destroyed by Vivian, we observed a large variability among sites. Almost a quarter of the windthrow areas had a basal area recovery of $\leq$ 6.6 \% and 2 \% of $\leq$ 1 \%. These areas pose a risk of a a reduced protective effect in long-term. However, since a large majority of these areas at risk are rather small ($\leq$ 0.6 ha in 75 \%), the loss of the protective effect might be bearable at slope scale.   

The recovery of the protective effect was most importantly influenced by temperature, with a generally higher recovery with increasing mean temperatures and day degree sum, and thus with decreasing elevation. This can be attributed to a shorter vegetation period at higher elevations \citep{wohlgemuth2008effekte, schonenberger2002post}. It is also in line with \cite{cerioni2024recovery}, who found a negative influence of elevation on the recovery of windthrow areas across Europe mainly for sites above 1000 m a.s.l.\\
The substrate appears to play a role as well with faster recovery on siliceous substrate or blocky debris compared to limestone. This is probably related to soil conditions and could be explained by either differing soil pH values or soil depth. \cite{kramer2014site} found a positive correlation between regeneration density and soil pH in Vivian windthrow areas. This would be contradicting to a lower recovery on limestone sites, which usually exhibit higher pH values compared to siliceous sites, for example. However, limestone slopes typically have shallow soils with a low water storage capacity, which could be a reason for a retarded forest recovery \citep{gessler2007potential, proll2016substrate}. This effect might be enhanced in combination with slope steepness, as limestone sites were significantly steeper compared to windthrow areas on sandstone and blocky debris. A potential limiting effect of drought at limestone sites is further underpinned by the fact that they showed a negative correlation of stem density and basal area recovery with drought stress in 2021. 

Overall, we observed contrasting effects of drought stress for the two observation periods. While we found a positive correlation of recovery with drought stress in 2003,  basal area recovery was not and stem density was negatively correlated with drought stress in 2021. \added{This was underpinned by the significantly negative influence of the interaction between drought stress and the time since the disturbance in the multivariate model.} On the one hand, this could be explained by the enhanced competition for water resources due to the increasing stem density \citep{aldea2024current, steckel2020implications}. On the other hand, this change in the effect of drought could also be related to an increase in drought frequency and intensity between the two observation periods. Maximum drought stress increased significantly between the periods 1990-2003 and 2004-2021, which may have led to drought increasingly becoming a limiting factor for regeneration.\\

The main limitations of this study stem from the relatively coarse temporal and spatial resolution of the large-scale data used. The available CHM's of the two time steps (2003 and 2021) distinctly differed in their point density, which raises questions about their comparability. We tried to overcome this by defining the recovery of the protective effect not as an absolute value (e.g., a specific basal area), but as a ratio to a reference area as an approximation of the pre-disturbance state of the windthrow area. \added{As a consequence, it is possible that the assessed recovery is high in certain areas, but the protective effect may not be sufficient since it was low for the pre-disturbance forest state.}
The predictor variables, such as climate or geological substrate, were also calculated from large-scale data, which did not allow for analyzing the effects of e.g., microclimate or specific soil conditions. Additionally, we were not able to capture other site effects, such as ungulate browsing or competition from ground vegetation, which might play an important role for the recovery after disturbances \citep[e.g][]{kramer2014site}. 
Furthermore, our study did not allow for conclusions regarding the potential temporary protective effect of disturbance legacies, such as lying deadwood. Several studies showed that lying deadwood can help to reduce the protection deficit in the first years after a disturbance event \citep[e.g.,][]{ringenbach2021full, caduff2022large, rammig2007wind}. \cite{baggio2022novel}, for example, developed and analyzed four indices based on a canopy height model to monitor the protective effect against snow avalanches in two wind storm areas in the Alps. They found a minimum level of protection ten years after disturbances, which presumably coincides with the minimum in the protective effect of the decaying deadwood and the still missing protective effect of the emerging regeneration. In a next step, the proposed indices could be applied to a subset of the windthrow areas analyzed in this study, for which more CHM's with a sufficient resolution are available.
Additionally, previous studies evidenced positive effects of lying deadwood on regeneration after disturbances \citep{bottero2013interacting, marangon2024impact, cerioni2024recovery}, albeit open questions regarding the influence of deadwood on forest recovery remain \citep{kramer2014site, thom2023regeneration}.
A major drawback of our study was the missing information on the post-disturbance management since it was not systematically documented on national level. Therefore, the effects of post-disturbance treatment (e.g., salvage logging or planting) could not be assessed.
Finally, the remote sensing data used in this study did not allow to assess and compare the species composition in the windthrow areas. Therefore, we cannot conclude whether the emerging species are well suited for future climate conditions and, thus, enhance the likelihood of a long-term and sustainable provision of the protective services. In general, disturbances may catalyze changes in species distribution and offer opportunities for climate-adapted species \citep{scherrer2022canopy}. \added{To assess the impact of climate change on the recovery of the protective effect, it would be interesting to extend and combine the results of this study with tools like dynamic vegetation models \citep{rammig2006forest, moos2023mountain}.}

\section*{Conclusions}
Our study on the recovery of protective forest stands 31 years after the Vivian wind storm of 1990 in Switzerland revealed an overall high resilience of the forest's protective effect after severe storm disturbances. It showed that the degree of recovery varied clearly among the considered forest structural variables. Especially forest cover and vegetation height demonstrated faster recovery rates compared to basal area and stem density, indicating a disparity in the restoration of the protective effect against different types of gravitational natural hazards. The retarded recovery of basal area highlights the extended time spans required for a full functional recovery of the protective effect against rockfall and landslide hazards. Conversely, the well advanced recovery of forest cover and tree height suggest a rapid restoration of the protective effect against snow avalanches. 
Our results underscore the importance of local environmental conditions, such as temperature, water availability, geological substrate or local topography when predicting forest recovery. \added{This insights on forest recovery after wind disturbances are not only relevant for the forest's protective effect, but also for, e.g., carbon storage, the water and bio-geochemical cycle, timber production or the susceptibility to future disturbances.}

The limitations of this study, including the reliance on large-scale remote sensing data and the lack of detailed post-disturbance management information, highlight areas for future research. Understanding the role of deadwood and the species composition as well as its suitability for a future climate,  remain critical for a comprehensive assessment of the recovery of the forest's protective service. 


\section*{Methods}

\subsection*{Study region}

The study was conducted with windthrow areas damaged by the storm Vivian in the Swiss cantons Valais (VS), St. Gallen (SG), Uri (UR), Glarus (GL) and Grisons (GR) (Figure \ref{fig:site}). They cover a large part of the Swiss Western to Eastern Alps and were chosen based on the following criteria: i) high impact of Vivian with considerable number of mapped windthrow areas; and ii) availability of canopy height models (CHM) of a satisfactory quality (point density of \textasciitilde 2 m\textsuperscript{-2} for two time steps after the storm (2003 and 2021).

\begin{figure}[ht]
\centering
\includegraphics[width=\linewidth]{figures/Site_overview.png}
\caption{Overview of the 1740 windthrow areas considered in this study (a) and an example of a windthrow (green) and reference (red) area (b). Background maps: $\copyright$ swisstopo.}
\label{fig:site}
\end{figure}

\subsection*{Data}

The windthrow areas were extracted from the national mapping product "forest damage areas Vivian 1990" provided by the Swiss Federal Office for Environment \citep{vivian}. It includes all windthrow areas with a high damage severity, defined as a maximum tree coverage of 20 \% of the remaining stocking. They were mapped based on areal images. We considered windthrow areas in the study cantons containing at least 75 \% protective forest based on the national delimitation of forests with a protective function \citep{Losey2013}. \added{Their median size is 0.5 ha (0.01-30 ha) and they are located mainly in the submontane to upper montane elevation belt (75 \% of the areas are between 1200-1600 m a.s.l.).}
They were compared to undamaged "reference areas", which were determined based on a dynamic buffer around the windthrow areas. Between the reference and the windtrow area, a 15-meter-wide zone was omitted to avoid edge effects of the windthrow event in the reference area. Thus, a dynamic buffer was applied to the storm-affected areas, which had already been buffered by 15 meters, using a factor of 2. This means that the resulting buffered area is equal in size to the input area. The resulting reference areas were finally clipped to the forested area based on a forest mask derived from the swissTLM3D ground cover \citep{swisstopo_swisstlm3d} for the cantons SG, UR, and GL, respectively forest masks provided by the cantons GR and VS (Figure \ref{fig:site}). Windthrow areas, where the reference area was smaller than 70 \% or larger than 500 \% of the windthrow area were removed. This resulted in a total of 1740 windthrow and reference areas.

To quantify the recovery of the protective effect, we used individual trees data extracted from airborne laserscanning (ALS) based canopy height models (CHMs) from 2003 (\textasciitilde 13 yrs after storm Vivian) and 2021 (\textasciitilde 31 yrs after storm Vivian). The ALS were available either from the Federal Office of Topography (swisstopo) or the cantons. The average point density of the ALS data ranges between 1.6 and 2.1 m\textsuperscript{-2} in 2003 and 23 and 31.4 m\textsuperscript{-2} in 2021. 
To derive realistic estimations of the diameter at breast height (dbh) of the individiual tree data generated from the CHM's, allometric functions were fitted based on data from the Swiss National Forest Inventory \citep{brandli2020}.

Explanatory variables for the recovery of the protective effect were mainly derived from a digital terrain model with a cell size of 2 m (swissALTI3D; \cite{swisstopo_dtm}). Additionally, we used gridded monthly temperature and precipitation data with a resolution of 100 m derived from the MeteoSwiss weather stations \citep{meteoswiss}. We further extracted geological data from the Swiss lithological-petrographic map \citep{swisstopo_geol}. 

\subsection{Quantification of the recovery of the protective effect}

The recovery of the protective effect was quantified based on different forest structural variables that are strongly linked to the forest's protective effect \citep[e.g.,][]{perzl2021protective, dupire2016novel, may2023quantifying}. A key variable in this context was stand basal area, as it proved to be a very good proxy for the protective effect of a forest against various natural hazard processes \citep[e.g.,][]{moos2017quantifying, berger2013eco}. Additionally, we analyzed vegetation height, stem number, tree diameter at breast height (dbh) and crown cover density (forest cover\added{, determined as the ratio of pixels $\geq$ 3 m and pixels $<$ 3m}). 
The forest structural proxy-variables were extracted from individual tree data (stem number, dbh, basal area) or the CHM (vegetation height, forest cover) for the windthrow and the references areas for two time steps after Vivian (2003 and 2021). The individual trees were detected with a hybrid model combining statistical modeling and Local Maxima-based methods as described in \cite{schaller2023improving}. It is important to note that the single tree data does not include regeneration, but only vegetation above 3 m height. It is further mainly representative of the dominant tree layer, while lower layers are often underestimated in the single tree detection \citep{schaller2023improving}. The dbh of the individual trees was determined based on allometric functions derived from NFI data. We therefore used single tree data from NFI plots in Swiss protective forests and classified them in nine categories according to their ecoregions (following \cite{frehner2005nachhaltigkeit}) and the observed degree of damages (not damaged and damaged).
We firstly selected the best model to predict tree diameter from tree height (Ht) among seven functions adjusted for the entire data set including power (\cite{feldpausch2011power}), logarithm, Michaelis-Menten (\cite{molto2014mmenten}), Chapman-Richards (\cite{pienaar1973crEqn}), exponential (\cite{meyer1940exponential}), Weibull (\cite{zhang2020weibull}) and logistic (\cite{richards1959logistic}) functions using the nls function in R (version 4.4.0, \cite{Rversion}).
The model performance was evaluated using the $R^2$ (i.e. the coefficient of determination of the model), the Akaike Information Criterion (AIC) and the Root Mean Square Error (RMSE).
The Weibull function showed the best performances on the entire dataset ($R^2$=0.86, AIC=22854.77, RMSE=13.04). However, the model did not converge for each of the considered categories. Thus, we used a logistic function weighted by 1/($Ht^2$) ($R^2$=0.86, AIC=23264.89, RMSE=13.06, Eq. \ref{eqn:DH_eq}).
\begin{equation}
\label{eqn:DH_eq}
 \operatorname{dbh} = \frac {\beta_{0}}  {1 + e^\frac { (\beta_{1} + \operatorname{Ht})} {\beta_{2} } }
\end{equation}

The model parameter $\beta_{0}$ is the maximum asymptotic diameter (in cm), $\beta_{1}$ is the midpoint value of the logistic regression, i.e. the inflection point and, $\beta_{2}$ is the scale parameter.

Although only severely damaged windthrow areas were considered in the analysis, it could not be excluded that the data contained occasionally trees that remained from the stand before the storm event. We therefore used yield charts to get a rough estimate of the upper possible limits of dbh and tree height after 13 and 31 years, respectively \citep{Badoux1983}, and removed trees exceeding this limits from the analysis.

Finally, we calculated for each time step (2003, 2021) the ratio of the proxy-variable between the windthrow and the reference area as a measure of recovery (e.g., basal area recovery of windthrow area $i$ $BA_{recovery,i,2003}$ = $BA_{storm,i,2003}/BA_{ref,i, 2003}*100\%$). We assumed maximum recovery to be 100\% in case the forest structural variable should reach a higher value in the storm compared to the reference area. Accordingly, the reference area was considered as representative for the protective effect of the windthrow area before the storm. Thus, "100 \% recovery" means that the protective effect of the forest before the storm has been regained, and not that the forests provides full protection. Accounting for the relative change between windthrow and reference areas allows for omitting biases due to different data quality in the CHM and individual tree data.

\subsection*{Statistical analysis}
We used \deleted{univariate and} multivariate statistical methods to identify factors influencing the recovery of protective forests after the windthrow damage. In particular, we analyzed the effects of topographic and climatic predictor variables as well as geology on the recovery of the proxy-variables for the protective effect (Table \ref{tab:ExplVar}). \added{In a first step,} the influence of \added{each} single predictor was tested \added{separately} using linear regression and non-parametric Kendall correlation tests \added{to identify meaningful variables for the multivariate model}. We \added{then} \deleted{further} built multivariate generalized linear models to analyze the effects of multiple variables on the recovery of the protective effect. Our main focus here was basal area recovery, since it is the most comprehensive proxy-variable of the protective effect. To avoid colinearity between predictors in the regression model, we calculated the Variance Inflation Factor (VIF) of each explanatory variable and excluded variables with a VIF $\geq$ 5. The VIF is defined as $1/(1-R_j^2)$, whereas $R_j^2$ is the coefficient of determination obtained when variable $j$ is regressed on all other predictor variables \citep{daoud2017multicollinearity}. It was calculated using the VIF function in the $regclass$ package in R \citep{regclassR}. \\
The GLM was fitted with a training data set representing 80 \% of the total data. The remaining 20 \% were used as test data set to evaluate its performance. \added{To separate the data in a training and test set, we applied spatial blocking using the function $cv\_spatial$ of the $blockCV$ package in R \citep{valavi2018blockcv}, which has proved as a robust method for estimating the predictive performance of models with spatial dependence structures \citep{roberts2017cross}.} \\
The final \deleted{model} \added{set of predictor variables} was determined based on stepwise backward variable selection aiming at reducing the Akaike Information Criterion (AIC). We also tested for potential interactions between the predictor variables. The model was examined based on customary residual plots \citep{fox2019regression}, which indicated a logarithmic transformation of the target variable (basal area recovery). \added{To detect potential spatial autocorrelation of model residuals, we further calculated Moran's I using k-nearest neighbors and distance-based weights and examined the semivariogram \citep[e.g.,][]{Dormann2007methods}. Since this indicated that the residuals of the GLM were spatially autocorrelated, we re-fitted the model with the $gls$ function of the $nlme$ package in R \citep{nlme} implementing a exponential correlation structure of the geographical position of the windthrow areas}. We finally calculated the $pseudo-R^2$ by \cite{cragg1970demand} \added{based on the likelihood function}\deleted{, given as the proportional reduction in deviance,} as well as the root mean squared error (RMSE) for the prediction of the test data.

\begin{table}[ht]
    \centering
\caption{Explanatory variables potentially influencing the recovery of the forest's protective effect that were considered in the statistical analysis.}
\label{tab:ExplVar}
\addtolength{\tabcolsep}{4pt}
	\fontsize{9.5pt}{10.25pt}\selectfont
	\renewcommand*{\arraystretch}{1.4}
    %\begin{tabular}{|p{0.15\linewidth}|p{0.15\linewidth}|p{0.2\linewidth}|p{0.2\linewidth}|}
    \resizebox{\textwidth}{!}{
    \begin{tabular}{|p{0.15\linewidth}|p{0.15\linewidth}|p{0.2\linewidth}|p{0.2\linewidth}|}
    \hline
    category & variable & explanation & source \\
    \hline
    \textit{geology} & lithology & main lithological group & lithological-petrographic map \citep{swisstopo_geol} \\
    \hline
    \textit{topography} & slope & mean slope & DTM 2m \citep{swisstopo_dtm} \\
         & northerness & cosinus of aspect & DTM 2m \\
         & easterness & sinus of aspect & DTM 2m \\
         & curvature & mean profile and planar  & DTM 2m \\
         & elevation & meter above see level & DTM 2m \\
    \hline
    \textit{climate} & temperature & mean seasonal temperature (1990-2003, 2003-2021) & \citep{meteoswiss} \\
         & daily degree sum & mean yearly sum of daily temperature $>5 \degree$ for 1990-2003 and 2004-2021 & \citep{meteoswiss} \\
         & precipitation & mean seasonal precipitation and min. precipitation during vegetation period (April-October) & \citep{meteoswiss} \\
         & drought stress & mean yearly index based on monthly ratio between evapotranspiration and the evaporative demand from the soil & monthly mean daily precipitation and temperature based on \citep{meteoswiss}; soil water storage capacity derived from \cite{Soilmap} \\
    \hline
         \textit{others} & gap size & size of windthrow area & \cite{vivian} \\
         & years since & time since the disturbance (13 years in 2003 and 31 years in 2021) &
    \end{tabular}
 }     
\end{table}

\section*{Data availability}
The data set of the recovery of different forest structural variables in the windthrow areas as well as the script for analyzing them is available here: https://doi.org/10.5281/zenodo.14794713   . 

\bibliography{sample}


\section*{Author contributions statement}
C.M. developed the study design, analyzed the data and wrote the manuscript draft,  K.D. contributed to the data preparation and analysis, A.E. contributed to the data preparation and analysis, C.G. provided and prepared the vegetation height models, C.S. run the individual tree detection, E.N. derived diameter-height relationships to determine the diameter of detected trees, L.D. contributed to the study concept and the initial draft writing. All authors reviewed the manuscript. 

\section*{Acknowledgments}
This research was supported by the Interreg Alpine Space Programme project MOSAIC co-funded by the European Union (grant number ASP0100014).


\section*{Competing interests}
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.


\end{document}