Assessing restoration success by predicting time to recovery—But by which metric?

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2019 The Authors. Journal of Applied Ecology published by John Wiley & Sons Ltd on behalf of British Ecological Society 1Department of Environmental Sciences, Western Norway University of Applied Sciences, Sogndal, Norway 2University of Oslo, Geo-Ecology Research Group, Section of Research and Collections, Natural History Museum, Oslo, Norway 3Radboud University, Institute for Water and Wetland Research, Nijmegen, The Netherlands 4Norwegian Institute for Nature Research, Bergen, Norway


| INTRODUC TI ON
Prediction play a minor role in most branches of ecology (Houlahan, McKinney, Anderson, & McGill, 2017), partly because of the complexity of ecosystems (Byers, 2018). However, to respond to the rapid changes in climate and ecosystems due to human activities (Ceballos et al., 2015), ecology needs to become more predictive (Evans, Norris, & Benton, 2012). This also applies to restoration ecology (Brudvig, 2017), since predictions of, for example time to recovery, are directly relevant for evaluation of progress in ecological restoration (Urban, 2006). Until recently, restoration ecology has lacked proper methods to predict time to recovery for informed evaluation of restoration success. The novel ordination regression based approach (ORBA; Rydgren, Halvorsen, et al., 2019b) facilitates prediction of time to recovery through modelling either linear or asymptotic trajectories of compositional change over time. In particular, the asymptotic model appears promising since it accounts for declining successional rates, a common characteristics of successions (Anderson, 2007;Rydgren, Økland, & Hestmark, 2004). Reliable predictions of time to species compositional recovery are relevant for assessing success in restoration because species composition constitutes one of the most important attributes of ecosystems, summarizing the outcomes of all important ecological processes (cf. Clewell & Aronson, 2013).
Although species composition may provide indications of a much wider range of ecosystem properties than aggregated, more generic properties of vegetation such as total cover and species richness, the latter are much more often used to assess restoration success (Ruiz-Jaén & Aide, 2005;Waldén & Lindborg, 2016). Which metrics to use is still vigorously debated (Abella, Schetter, & Walters, 2018;Brudvig et al., 2017;Reid, 2015). In general, we expect species composition to recover more slowly than total cover and species richness. Predictions of restoration success is therefore expected to depend on the metrics considered (cf. Crouzeilles et al., 2016). Still, the view that simpler measures (cover, richness) may be appropriate as proxies for more complex properties such as species composition is not uncommon (Brancalion & Holl, 2016;. A closer examination of the idea that simpler measures can approximate species composition for assessing restoration success by predicting time to recovery, is therefore needed. Soil conditions may play an important role in successful restoration of species composition (Piqueray et al., 2011;Rydgren, Halvorsen, Auestad, & Hamre, 2013). Because abiotic and biotic ecosystem components mutually influence each other, both have to be taken into account in restoration of ecosystems. In harsh environments, such as alpine environments where the majority of hydropower spoil heaps are found, soil processes are slow (Kidd, Streever, & Jorgenson, 2006). Accordingly, alpine spoil heaps have low levels of soil organic matter even after eight decades (Rydgren et al., 2013), and slower recovery of soil properties than of univariate biotic ecosystem properties such as species richness and total cover is therefore expected.
In this study, we use ORBA with both linear and asymptotic models to analyse data from three censuses of species composition and environmental conditions of alpine spoil heaps and their undisturbed surroundings, over a period of 24 years, 7-41 years after the initial disturbance, to predict time to recovery of species composition.
ORBA results are compared with parallel analyses of time-to-recovery predictions for selected properties of the physical environment, total cover and species richness. The studied spoil heaps consist of surplus rock material resulting from tunnel construction for hydropower plants. In alpine areas, such spoil heaps pose serious restoration challenges due to the harsh environment (Rydgren et al., 2013). Previously, full species compositional recovery of the four alpine spoil heaps addressed in the present study was predicted to be achieved in less than 50 years (Rydgren, Halvorsen, Odland, & Skjerdal, 2011). However, the analyses of Rydgren et al. (2011) were based upon time-to-recovery predictions derived from a linear model of compositional change over time between two time-points.
The resulting predictions are likely to be over-optimistic as constant successional rates are unrealistic (Rydgren et al., 2011). The present study is based on three censuses. Since we expect successional rates to decrease with time, we hypothesize that (a) species compositional recovery will be slower than indicated by (Rydgren et al., 2011) (hypothesis 1) and we further expect (b) species compositional recovery to be slower than the recovery of functionally less specific properties such as total cover or species richness (hypothesis 2). Finally, we expect (c) that recovery of the physical environment will be slower than recovery of total cover and species richness (hypothesis 3).

| Field sites
We studied four low alpine spoil heaps in western Norway, situated from 1,000 to 1,360 m a.s.l., within an extent of 30 km (Figure 1, see Figure S1, all of which were also included in a previous study (Rydgren et al., 2011). The climate is relatively wet and cold with encourage further improvement of methods like the ordination-regression based approach that use species compositional data to predict time to recovery. annual precipitation in the range 1,500-1,900 mm and mean July temperatures of 6.6-8.8ºC for the period 1971period -2000period (NVE, 2018. The spoil heaps were constructed between 1974 and 1984, and vary in size from 2.7 to 4.1 ha. They consist of blasted rocks, except Kleådalen which consists of finer-grained substrate (Rydgren et al., 2011). The spoil heaps are made up by bedrock material similar to that of their surroundings, mainly consisting of gneisses, granites, phyllite and other metamorphosed rocks (Skjerdal & Odland, 1995).
Shortly after construction, compound fertilizer was added to all spoil heaps, followed by seeding with commercial seed mixtures annually for at least 3 years (Skjerdal & Odland, 1995). Additional fertilization and seeding also took place later years (see Appendix S1 for details). All spoil heaps have been sporadically grazed by sheep since construction.

| Sampling design and data collection
Data on the physical environment (soil organic matter and pH), total cover (vascular plants, and bryophytes and lichens), species richness (vascular plants, bryophytes and lichens), and species composition were collected at all sites in the early 1990s (Skjerdal, 1993;Skjerdal & Odland, 1995), in 2008 (Rydgren et al., 2011. In the early 1990s (1991in Svartavatn, Fossane and Kleådalen and 1994 in Øydalen), 10-20 non-permanent sample plots (0.5 × 0.5 m each) were placed on each spoil heap by stratified random sampling (64 plots in total), using a baseline approach (Skjerdal, 1993). In 2008, we selectively placed 8 blocks per spoil heap (5 blocks on the heap and 3 in their undisturbed surroundings), each 5 × 10 m. Within each block, three permanent sample plots (0.5 × 0.5 m) were placed at random, with the extra condition that plots had to be separated by at least 1 m (Rydgren et al., 2011). All permanent plots were re-analysed in 2015 except one block at the Kleådalen spoil heap, which was lost due to deposition of rock material. The resulting dataset comprised 253 sample plots with 239 recorded taxa (107 vascular plants, 86 bryophytes and 46 lichens; see Appendix S2 for nomenclature).
The species composition of the plots was recorded in July-September at the three sampling occasions. We divided each plot into 16 equally sized subplots and recorded the abundance of vascular plants, bryophytes and lichens as the frequency in 16 subplots. We Soil samples were collected from the upper 5 cm soil layer at the three sampling occasions, and analysed for soil organic matter in % and pH (see Appendix S3 for details).

| Data analyses
All statistical analyses were performed in R versions 3.2.2 or 3.5.3 (R Development Core Team, 2019).
We used ORBA to predict time to species compositional recovery (Rydgren, Halvorsen, et al., 2019b). This entailed first extracting the gradient structure of a species compositional dataset by parallel use of detrended correspondence analysis (DCA; Hill & Gauch, 1980) and global non-metric multidimensional scaling (GNMDS; Minchin, 1987) Table S1). Second, we calculated successional distance (d jt,0 ) along these gradients, that is the distance between a restoration plot and its reference, using the mean position of plots from the surroundings at a given time point as our dynamic reference. Since 2008 was the first year plots from the surroundings were censused, we used the 2008 plots as references also for plots from the 1990s, implicitly assuming that the reference vegetation was as stable between these years as it was between 2008 and 2015 (Table S1).
Finally, we modelled successional distance as a function of time since disturbance with generalized linear mixed models (GLMMs; Bates, Mächler, Bolker, & Walker, 2015) with identity-link and Gaussian error distribution, by fitting two types of responses, a linear relationship (M L ), and a log-linear, asymptotic relationship (M A ).
From the models we derived three different time-to-recovery predictions for each spoil heap: TR L0 (linear model; restoration target: predicted successional distance = 0), TR As+1 (asymptotic model; restoration target: predicted successional distance = +1 S.D. off the centroid of reference plot scores) and TR Af+0.01 (asymptotic model; restoration target: predicted successional distance = 0.01).
We modelled the physical environmental variables and the biotic variables as functions of time since disturbance (year of spoilheap construction). Soil organic matter, vascular plant cover and bryophyte and lichen cover were expressed as percentages, that is as strictly bounded but non-binomial data. We therefore logit-transformed these variables (Warton & Hui, 2011) before modelling them with identity-link and Gaussian errors. The species richness variables were modelled with log-link and Poisson errors, whereas pH was modelled with identity link and Gaussian errors. Finally, we used the models to predict the number of years before reference levels (mean values in the undisturbed surroundings) are reached.
All models were parameterized using the r package LME4, version 1.1.21 (Bates et al., 2015), accounting for repeated measurements and the spatially nested sampling structure by specifying plot nested within site with block as random effects in the models. The fact that the spoil-heap plots were not set up in a block design in the censuses of the 1990s was handled by allocating the 1990s plots randomly to blocks to fit the structure of the model. For the species richness variables the mixed-effects models did not converge, even when the random effect was reduced to only consist of sites. The preliminary estimates from these models showed that the random effects hardly contained any variation. Therefore, we used generalized linear models without random effects for the species richness variables. Time to recovery estimates from these models therefore have to be interpreted more conservatively.

| Species compositional change during restoration
The first axis of the DCA and GNMDS ordinations represented the restoration successional gradient (Table S1) In the early phase, the species composition of all spoil heaps was dominated by pioneer species like Bryum spp., which decreased rapidly over time both in frequency and abundance (Table S2) Time-to-recovery predictions for the spoil heaps, obtained by ORBA, were between 43 and 60 years for linear models (TR L0 ) and in the ranges 41-105 and 115-212 years for the asymptotic models TR As+1 and TR Af+0.01 respectively ( Figure 3; Table 1). Generally, longer time-to-recovery predictions were obtained by use of successional distances obtained from DCA than from GNMDS ordination, but the differences between ordination methods were generally small (two of the predictions based upon TR Af+0.01 excepted). The shortest and longest predictions were most often obtained for the Kleådalen and Svartavatn spoil heaps respectively (Table 1).

| Dynamics of total cover
Patterns of change in vascular plant cover over time differed considerably among the four spoil heaps (Figure 4). Significant increase over time was observed only at Kleådalen and Svartavatn where reference levels were predicted to be reached in 43 and 55 years after disturbance respectively ( Table 2).
The cover of bryophytes and lichens increased rapidly after disturbance at all spoil heaps and had reached reference levels at all sites in 2015, 31-37 years after disturbance (Figure 4, Table 2).
Particularly rapid recovery of bryophyte and lichen cover was observed at Kleådalen, where the reference level was reached already at the first census in 1991, 10 years after disturbance.

| Dynamics of species richness
The number of vascular plants (species richness) increased significantly over time at all spoil heaps except Øydalen ( Figure 5). Predictions from the significant models indicated that the reference levels would be reached in 36-58 years after disturbance in these spoil heaps ( Table 2).
The number of bryophyte species increased more rapidly than vascular plant species numbers and reached reference levels at the second census, 13-28 years after disturbance ( Figure 5). The number of lichen species followed the same pattern as bryophytes; reference levels were reached at the second or third census ( Figure 5).

| Dynamics of the physical environment
Soil organic matter increased significantly over time at all four spoil heaps (Figure 6), but was still far below reference levels at 100 years was predicted to be required after disturbance before reference levels are reached (Table 2). Shortly after disturbance, pH was much higher at the spoil heaps than in their surroundings but, with the exception of the Svartavatn spoil heap, decreased significantly over time ( Figure 6). For these three spoil heaps pH was predicted to reach reference levels 67-92 years after disturbance.

| Time to recovery estimation using ORBA
Our results show that successional rates in the four studied spoil heaps TA B L E 1 Predicted time to species compositional recovery, measured in years, as obtained from three different ordination-regression based approach (ORBA) models: TR L0 (linear model; restoration target: predicted successional distance = 0), TR As+1 (asymptotic model; restoration target: predicted successional distance = +1 SD off the centroid of reference plot scores) and TR Af+0.01 asymptotic model; restoration target: predicted successional distance = 0.01) 2010; Myster & Pickett, 1994). Furthermore, we obtain a prominent successional gradient along the closely similar first axes of the GNMDS and DCA ordinations which justifies their use for predicting time to recovery (Rydgren, Halvorsen, et al., 2019b) and ensures that the basic assumption of the ordination-based approach for predicting time to recovery (ORBA), that a proxy for the successional gradient is available, is satisfied (Rydgren, Halvorsen, et al., 2019b).  time into account by applying asymptotic models of successional distance, obtaining predictions for time to full compositional recovery from 115 to 212 years, indicates that the major reason for the difference is that successional rates decrease over time (Foster & Tilman, 2000;Myster & Pickett, 1994;Rydgren, Halvorsen, et al., 2019b;Rydgren, Halvorsen, Töpper, & Njøs, 2014). Although it is not possible from our data to conclude which of the models, linear or asymptotic, that gives the most accurate predictions, two strong arguments point in favour of the asymptotic models. First, only the asymptotic models account for decreasing successional rates over time, that is that successions stand out as logarithmic processes.

Site
Second, older spoil heaps in the region are still far from full recovery and more than 100 years seem to be needed for full recovery due to the common construction practice of using a coarse top substrate (cf. Rydgren et al., 2013). The definition of recovery also influences predictions of time to recovery. To reach near-zero successional distance from the reference state (TR Af+0.01 in this study) will of course take longer time than to satisfy a more relaxed criterion (TR As+1 in this study). In fact, full recovery in the sense defined here may not be possible, feasible or even relevant in all cases (McDonald, Gann, Jonson, & Dixon, 2016). Therefore, adopting a relaxed criterion like +1 standard deviation off the centroid of reference plot scores may sometimes be appropriate.
Our results show that ORBA predictions are influenced by the robustness of ordination methods (DCA and GNMDS), which in turn depends on properties of the data (Eilertsen, Økland, Økland, & Pedersen, 1990;Minchin, 1987;Økland, 1990;Rydgren, 1993) such as the distribution of species´ frequencies. In particular, species that occur in few plots are known to influence the ordination result, and especially so in species-poor or otherwise deviant sample plots (Økland, 1990). We therefore encourage studies on the robustness of ORBA to differences in dataset properties. So far, ORBA has only been tested on a boreal forest dataset (Rydgren, Halvorsen, et al., 2019b). Also in that case, asymptotic models were superior to linear models, and the asymptotic GNMDS-based models fitted the data better than the asymptotic DCA-based models (Rydgren, Halvorsen, et al., 2019b). More studies using ORBA are, however, needed before eventual final conclusions about choice of ordination method can be drawn. A likely outcome of further tests is, however, that no method is consistently better than the other (as for ordination of species data; Økland, 1990;van Son & Halvorsen, 2014) and that two methods should always be used in parallel to secure robustness of the results.

| Time to recovery using univariate metrics
Predictions of time to recovery based upon more generic, univariate metrics are substantially lower than those obtained from ORBA predictions of compositional recovery. Furthermore, the uncertainty of estimates based upon these metrics is larger, confirming hypothesis 2. This result may be due to (a) insufficient understanding of the ecological processes, (mis)leading us to opt for too simplistic models, or (b) that the relevant processes do not follow a simple monotonic pattern.
While both species richness and total cover of bryophytes and lichens increase relatively fast after disturbance, reaching reference levels within 13-58 years after disturbance, recovery of soil organic matter and pH is predicted to take up to 111 years. The relatively fast recovery of species richness, particularly for bryophytes, is due to the appearance of pioneer species as well as species typical of established vegetation relatively soon after disturbance. Later, as seen at the last census in 2015, bryophyte species richness drops when pioneer species go locally extinct. This shows that recovery of species richness is decoupled from species compositional recovery and, accordingly, that richness recovery does not imply recovery of an ecosystem with functions and a species composition in dynamic equilibrium with its environment.
The two cover variables follow a similar pattern; in two of the four spoil heaps the reference level is predicted to be reached relatively fast, within 55 years. The reason is that the cover of pioneer and established species follow the same pattern as species richness with shifts in the species groups that contribute to cover over time. In addition, the highly persistent seeded grass species contribute considerably to the total cover over the entire time period (Rydgren et al., 2016).
Our results indicate that univariate metrics for the physical environment (here: soil organic matter and pH) are more relevant for modelling the recovery process than properties of the biotic environment. There are two reasons for this: (a) that environmental conditions known to be important for the species composition are important in their own right when the recovery status is assessed; and (b) that the physical properties seem to follow a monotonic trajectory (e.g. continuously increasing soil organic matter content and decreasing soil pH). Based upon metrics for the physical environment, recovery is predicted to take place many years later than predicted from the univariate, biotic variables, and closer to the predicted recovery of the species composition, confirming hypothesis 3.
Typically, dominant species in the surroundings such as the Vaccinium spp. prefer soils with high content of organic matter and high moisture retention capacity for successful recruitment (Eriksson & Fröborg, 1996). For these species, conditions for establishment and growth at spoil heaps are still rather unsuitable 40 years after spoil heap construction. Soil development therefore appears as a major bottleneck for recovery (Rydgren et al., 2013(Rydgren et al., , 2011 of the species composition on alpine spoil heaps, emphasizing why knowledge about soil development is crucial for understanding recovery processes (Forbes, Ebersole, & Strandberg, 2001).

| Time to recovery-which metric to use?
The four categories of metrics used to predict time to recovery in this study give widely different results. The relatively fast recovery of species richness and total cover contrasts the slow recovery of properties of the physical environment and the species composition. Although species richness is the most commonly used metric to assess restoration success (Waldén & Lindborg, 2016), it is not regarded as a core metric for evaluating restoration outcomes (SER, 2004). Our study clearly demonstrates that species richness and total cover may reach pre-disturbance levels and suggest successfully accomplished restoration, whereas the succession of plant species and re-establishment of environmental conditions are still in relatively early stages. In contrast, the plant species composition is an indicator of the basic properties of the ecosystem and should therefore be regarded as 'the principal obligation of restorationists' (Clewell & Aronson, 2013). We advocate using species compositional data and soil conditions to predict time-to-recovery, and advise against using species richness and total cover as unique metrics of successful restoration (cf. Piqueray et al., 2011). However, we acknowledge that species richness and other generic variables may be valuable metrics shedding light on the restoration process when restoring the species composition is not feasible (Halpern, Antos, Kothari, & Olson, 2019). For example when ecosystems are considerably degraded, appropriate reference systems do not exist and full recovery will neither be attainable or desirable (McDonald et al., 2016). Another promising development for estimating time-to-recovery may be to use functional-trait based metrics, which may serve as indicators of ecosystem function (cf. Funk et al., 2017) and, therefore, also provide relevant indicators of recovery status.
The ordination-regression based approach (ORBA) used for predicting time to species compositional recovery in this study provides promising results, but further exploration, most notably of later phases of restoration successions, is needed to generalize about the precision of the predictions. Anyway, our results clearly support the view that, to advance restoration ecology, we should move beyond simple biodiversity measures (Brudvig, 2017) and develop methods for predicting time to recovery that take the full plant species composition into account (cf. Urban, 2006).

AUTH O R S ' CO NTR I B UTI O N S
K.R., I.A., L.N.H. and J.S. collected the data. K.R. and J.P.T. analysed the data; K.R. led the writing of the manuscript. All authors contributed critically to the draft and gave final approval for publication.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data available via the Dryad Digital Repository https ://doi. org/10.5061/dryad.hdr7s qvcr (Rydgren, Auestad, et al., 2019a).