Defining and quantifying the resilience of responses to disturbance: a conceptual and modelling approach from soil science

There are several conceptual definitions of resilience pertaining to environmental systems and, even if resilience is clearly defined in a particular context, it is challenging to quantify. We identify four characteristics of the response of a system function to disturbance that relate to “resilience”: (1) degree of return of the function to a reference level; (2) time taken to reach a new quasi-stable state; (3) rate (i.e. gradient) at which the function reaches the new state; (4) cumulative magnitude of the function (i.e. area under the curve) before a new state is reached. We develop metrics to quantify these characteristics based on an analogy with a mechanical spring and damper system. Using the example of the response of a soil function (respiration) to disturbance, we demonstrate that these metrics effectively discriminate key features of the dynamic response. Although any one of these characteristics could define resilience, each may lead to different insights and conclusions. The salient properties of a resilient response must thus be identified for different contexts. Because the temporal resolution of data affects the accurate determination of these metrics, we recommend that at least twelve measurements are made over the temporal range for which the response is expected.


Model development
Following this analogy of a mechanical spring and damper system and applying Newton's second law of motion, the dynamics of the response to disturbance are described by the ordinary differential equation: where y (m) is the displacement, (kg) is the mass, is the damping constant (which when multiplied by the velocity gives the damping force), is the spring constant (which when multiplied by displacement gives the force from the spring), (m) is the final displacement (which is used in order to represent the expected small changes in soil functioning after a perturbation) and H is a unit Heaviside step function. To eliminate one parameter this can be written as: where = √ ⁄ (s -1 ) is the natural frequency (which is the frequency at which the system would oscillate if no damping force were applied) and = 2√ ⁄ is the damping factor.
Note that for values of damping factor lower than one the mass would oscillate before returning to equilibrium, but we do not consider this response here. The initial condition = implies that an impulse disturbance is imposed on the system at time 0, such that the gradient of the function changes instantaneously at this point. The disturbance to the system is thus represented by a combination of this impulse and the step change represented in equation S2, both of which are considered to occur simultaneously.
Metrics of the resilience characteristics are quantified from this fitted model. As such, the degreee of return is defined as: In some contexts it may be desirable to calculate | |, if it is only a change in value from the reference level that is important rather than the absolute value.
The damping factor ( ) provides an indicator of the temporal response (on a dimensionless scale) after perturbation. = 1 (critical damping) corresponds to the quickest return to new equilibrium and increasingly larger values correspond to longer times for return to equilibrium. When applied to data, the peak response may occur at different times (e.g. and the damping factor should be scaled by the frequency parameter, . We therefore quantify the return time characteristic as: Whilst the expression ⁄ is a characteristic time associated with Equation 1, the multiplication by a factor 5 is used so that this time corresponds to approximately 95% return of the signal from to the largest perturbation (peak height) to the final return point (degree of return) for a dimensionless, critically damped system with = 0. Comparisons in the characteristic time identified for different soils could be made without using this factor, but we use this 'return time' to develop metrics for the other resilience characteristics.
An alternative to scaling the damping factor relative to the frequency parameter would be to scale relative to the time to peak which can be derived by finding the point at which the gradient of the curve is equal to zero. However, this would increase the sensitivity of results to the accurate estimation of the time to peak, which is likely to require higher resolution of data around the expected time of this peak.
The return rate is quantified by calculating the average magnitude of the gradient during the response period: Which can be solved numerically or analytically by identifying a turning point (if any) at which the gradient equals zero.
As the degree of return is already quantified by , here is calculated as the area under the response curve assuming that it returns to it's initial state (i.e. = 0). If the degree of return was also an important component of the 'efficiency' the two resilience characteristics could be considered to give a more complete description of the response. Note that we calculate the area within the time range from zero to twice the return time as beyond this value the reponse tends to zero, thus: In addition to reducing the sensitivity to the time range, this formulation removes the difficulty of interpretation that could occur if the response crosses the reference state, as in this situation areas below the x-axis will have a negative value and above the axis will be positive and may thus cancel each other out.
If the system is assumed to return eventually to its initial state (i.e. x=0 in S1), a simpler set of solutions is found with an overdamped, critically damped and underdamped responses given respectively as Parameter Rr was large in almost all of the measured data and thus justifies the more complex formulation (S3).

Model fitting and estimation of resilience characteristics
The model (Eq. 1) was fitted to data in order to estimate the four parameters ( , , and ). This fitting was conducted in two stages. First, the parameters were optimised within a constrained range of parameters using the MATLAB function lsqcurvefit which uses iterative least squares. This optimisation was started from multiple points (50-100) within the constrained region, so that a global rather than a local optimimum was found. The parameter was constrained such that 1/ was larger than half the measurement interval, or half of the time to the first measurement. This was done so that the fitted parameter values did not simulate an immediate peak in the simulated response. Whilst low resolution data could mean that a peak of this kind does in fact occur more than half of the time before the initial measurement is made, this cannot be identified by the model without additional data. The optimum parameters from the first, constrained routine were then used to initialise a further, unconstrained fitting algorithm (nlfit) in order to estimate the parameters and the covariance matrix of these parameters. If the first stage of the optimisation process identified that the response was likely to be critically damped (i.e. = 1). This parameter was fixed and other three parameters were optimised in the second stage. This improved the convergence of the model.
The parameter estimates and the covariance matrices were then used in a Monte Carlo simulation of 10000 points in order to estimate the probability density function of the four resilience characteristics calculated using equations A1.4-A1.7. This error propagation simulation was performed in @Risk software. In the Monte Carlo simulation, the parameter values for and were each sampled from a normal distribution. Values of the parameters and , which have constraints at 0 and 1 respectively, were sampled from lognormal distributions. The correlation between the four parameters was accounted for during sampling. The median, interquartile range and the range (excluding outliers) of each of these simulated populations were then calculated and represented estimates of the different characteristics.

Synthetic data generation
Synthetic data was first generated using equation A1.1 by adding two types of noise to the signal.
The first component of the noise was assumed to occur in the four model parameters ( , , and ), simulating a time series of replicates from which repeated measures were taken. A second term ( ) represented white noise in the data and was sampled from a normal distribution with = 0, = 2. Each synthetic time series ( ℎ ) was calculated using the four parameters ( , , and ) and the error terms at a series of time points, , using: ℎ = ( + , + , + , + , ) + ( ) Data was simulated for three replicates (i.e. three values of the parameter errors).