Understanding how ecosystems respond to changes in environmental conditions and how resilient they are to perturbations is an urgent necessity in the context of global change1. Theoretical studies have proposed that the lengthening of the time needed to recover from a small perturbation can be used as an indicator for deteriorating resilience and can possibly be an early warning for the imminent collapse of an ecosystem2,3. This phenomenon (Fig. 1) of an increasing return time with increasing stress is referred to as “critical slowing down”. Although applying this concept seems straightforward, it is largely developed in the realm of theoretical models and translating the concept to practical applications is challenging. For example, stability analysis to get the eigenvalue is often performed in theoretical studies2,3. Here, an equilibrium state is perturbed with a negligible small disturbance. In real-world ecosystems this procedure is difficult to replicate. If too big disturbances are used, the ecosystem can be pushed beyond the tipping point as the results of the experimental disturbance itself and no reliable measure of recovery is obtained. Thus, disturbance type and magnitude need to be carefully chosen. Yet, even more challenging is how to measure and monitor recovery in real-world ecosystems after the disturbance is applied.
Most theoretical and empirical studies follow a nondestructive monitoring approach as proposed by ref.  for estimating the recovery rate, λ. Here, an exponential model is fitted against the time series of the biomass development that is obtained of the recovery after the disturbance (Fig. 2). In the empirical studies that tested resilience indicators so far (e.g. ref. [4, 5]), the state of the system could be relatively simply monitored without disturbing the biomass. For instance, in ref.  the researchers used the light attenuated as a proxy for biomass. By measuring the light attenuation in the mesocosm frequently and at regular intervals the changes in biomass could be traced. Likewise, in ref.  the concentration of cells could be monitored without disturbing the cell numbers. A second method, often used in theoretical studies employing numerical simulations, is to measure the time it takes before the system is recovered to the pre-disturbed state (Tr in Fig. 2, e.g. used in ref. ). The biomass has to be recovered within a certain accuracy around the pre-disturbed value (e.g. within 0.01 or 0.05%) before the disturbance-recovery experiment is terminated and the recovery time is recorded. However, these approaches are often impractical for the in situ assessment of resilience.
Measuring recovery rates reliably to estimate the resilience in situ in real-world ecosystems is often challenging due to: 1) difficulties in getting frequent and regular field observations of biomass; 2) changing environmental conditions such as seasonal dynamics (Fig. 3); and 3) the high level of heterogeneity and stochasticity in many real-world ecosystem (Fig. 3). These above points pose some serious challenges to the approach for estimating the recovery rates or recovery times in perturbed systems as proposed in the theoretical studies3,6 or as performed in controlled laboratory environments4,5. Monitoring biomass (i.e. vegetation) development frequent and at regular intervals can be particularly challenging due to the inaccessibility of the environment and the accuracy of the available methods to do so. For example, the ability to access the field sites of intertidal ecosystems shifts from day to day, and from one week to the other due to the tides, making recurrent regular field visits logistically difficult. Moreover, non-destructive biomass measures, such as canopy height or coverage can be impractical and laborious to acquire and can be inaccurate at the small scale of the disturbances due to the high level of variability. Automated observations, e.g. with fixed camera’s, still need extensive calibration and might suffer from the harsh hydrodynamic conditions and biofouling. For the second methodology (i.e. monitoring of recovery time) all the above objections remain and are supplemented with the fact that the duration of the experiments is determined by the recovery time itself, which is impractical. Therefore, an alternative and easy to execute approach was developed.
To circumvent the above issues, we measure recovery rates using a destructive sampling approach. After a fixed time interval since a disturbance is applied, the recovery in the disturbed plots is compared with a proxy for the equilibrium biomass, Veq, derived from undisturbed control plots at the same stress level (Fig.3, and see next paragraph). This approach is comparable with measuring net primary productivity over a certain time interval Δt. Like in ref. , we assume that during the recovery period the development approximates an inverse exponential:
V(t) = Veq – De-λΔt (1)
Here, V(t) is the biomass at time Δt, Veq is the equilibrium biomass, and D the disturbance as the amount of biomass removed. To find the recovery rate after a fixed time interval Δt we can normalize the function, if we define the relative recovery f as V(t)/Veq. In that case equation  can be written as:
f = 1 - de-λΔt (2)
in which d is the relative disturbance magnitude, defined as D/Veq. This leads to the rearrangement that allows to calculate the recovery rate:
λ = -log((1 - f)/d)/Δt (3)
High levels of heterogeneity and stochasticity in natural ecosystems can pose challenges for the estimation of the equilibrium biomass, Veq, to which the recovery measurements are compared (Fig. 3). Preferably, the recovery of a disturbed plot is immediately compared to a control plot (i.e. paired comparisons). This works well in case the relative recovery f is much lower than 1 for all replicates and when the response of the disturbance and control plot in a pair is closely linked. However, due to the high variability of biomass in the field in both the control as well as the disturbance plots, it can occur that f>1 for some pairs of disturbance and control plots. This leads to incorrect estimation of the recovery rate λ. Dismissing these measurements as outliers is not desirable as it reduces the power of the statistical analysis. However, if the disturbance and control plots are sufficiently independent (e.g. corroborated by the high variability in biomass between plots and the low neighborhood correlation) the issue can be resolved by estimating the mean (or maximum) biomass of the controls as proxy for the equilibrium biomass (see point 4 of Procedure and Troubleshooting for further elaboration).
The protocol, and the equipment and materials used as described below, is developed specifically for disturbance-recovery experiments in tidal marshes where inundation by seawater is the main stressor. However, this methodology can be applied to any vegetated ecosystem in which the main stressor is clearly defined, and it should be easy to modify the protocol for other (non-vegetated) ecosystems.