Quantifying the transient shock response of dynamic agroecosystem variables for improved socio-environmental resilience

In classic resilience thinking, there is an implicit focus on controlling functional variation to maintain system stability. Modern approaches to resilience thinking deal with complex, adaptive system dynamics and true uncertainty; these contemporary frameworks involve the process of learning to live with change and make use of the consequences of transformation and development. In a socio-environmental context, the identification of metrics by which resilience can be effectively and reliably measured is fundamental to understanding the unique vulnerabilities that characterize coupled human and natural systems. We developed an innovative procedure for stakeholder-friendly quantification of socio-environmental resilience metrics. These metrics were calculated and analyzed through the application of discrete disturbance simulations, which were produced using a dynamically coupled, biophysical-socioeconomic modeling framework. Following the development of a unique shock-response assessment regime, five metrics (time to baseline-level recovery, rate of return to baseline, degree of return to baseline, overall post-disturbance perturbation, and corrective impact of disturbance) describing distinct aspects of systemic resilience were quantified for three agroecosystem variables (farm income, watertable depth, and crop revenue) over a period of 30 years (1989–2019) in the Rechna Doab basin of northeastern Pakistan. Using this procedure, we determined that farm income is the least resilient variable of the three tested. Farm income was easily diverted from the “normal” functional paradigm for the Rechna Doab socio-environmental system, regardless of shock type, intensity, or duration combination. Crop revenue was the least stable variable (i.e., outputs fluctuated significantly between very high and very low values). Water-table depth was consistently the most robust and resistant to change, even under physical shock conditions. The procedure developed here should improve the ease with which stakeholders are able to conduct quantitative resilience analyses.


Defining resilience
The term "resilience" has been used in a narrow sense to refer to the rate at which a perturbed system is restored to equilibrium; in a slightly broader context, it has been interpreted as the postdisturbance rebound time or the degree of functional recovery to a baseline of performance. Recently, resilience has emerged as a cognitive framework for understanding how dynamic systems self-regulate and evolve over time. Since the first published definition of ecological resilience by Holling (1973), researchers in the social and natural sciences have gained a vastly improved understanding of "resilience thinking," which has informed the enhancement of disaster mitigation strategies, resilient infrastructure, personal and communal coping mechanisms, and adaptive capacities. Here, we define resilience as the combined ability of a system component variable to resist and efficiently recover from an array of social-environmental shocks (i.e., disturbances in variable behavior that force the system to operate outside of its normal functional paradigm).

Measuring resilience
Many research teams, special interest groups, official government entities, and think tanks have adopted similar yet slightly divergent heuristics for understanding resilience in socialecological systems (Gunderson 2002, Walker et al. 2006, Angeler and Allen 2016, Asadzadeh et al. 2017, Allen et al. 2018, Salomon et al. 2019, Cains and Henshel 2020. Several of these heuristics are used by different groups to frame their specialized definition(s) of systemic resilience; for example, Folke et al. (2010) outline resilience as social-ecological persistence, adaptability, and transformability, whereas Gallopín (2006) defines the linkages between vulnerability, resilience, and adaptive capacity. Although many of these definitions have overlapping elements, there still seems to be a lack of consensus regarding which aspects or behaviors of a system best exemplify resilient patterns. Several researchers have begun measuring resilience with specific respect to dynamic agroecosystems. Agroecosystem resilience has been assessed by applying ecological resilience-based (e.g., Peterson et al. 2018) or behavior-based (e.g., Cabell and Oelofse 2012) indicator frameworks or, as in our case, by using the stability, resistance (robustness), and recovery of system processes as a basic framework for resilience monitoring (e.g., Hodgson et al. 2015, Oliver et al. 2015, Lamothe et al. 2019, Bardgett and Caruso 2020).

Quantitative measurement methods
Here, we develop an approach for quantifying the resilience of three socio-environmental variables (farm income, water-table depth, and crop revenue; Appendix 1) in the Rechna Doab watershed in Pakistan. The approach was informed by the work of previous research teams (e.g., Hodgson et al. 2015, Nimmo et al. 2015 who developed independent but related methods for successfully quantifying the effects of anthropogenic pressures on ecological systems. Hodgson et al. (2015) and Nimmo et al. (2015) suggest mapping behavioral-https://www.ecologyandsociety.org/vol26/iss2/art17/ response metrics onto a bivariate state space with joint consideration for the "resistance" and "recovery" characteristics of a system. Ingrisch and Bahn (2018) propose a similar method for disturbance response measurement by using the normalized impact of disturbance and the normalized recovery rate to define the bivariate space. Similar to the preceding studies, we use quantifiable metrics for developing a comprehensive understanding of the resilience of a given system or variable. The five metrics identified as salient characteristics of a resilient functional response to disturbance (i.e., degree of return, rate of return, perturbation, time to return, and corrective impact) are based on the metrics used by Cimellaro et al. (2010), Todman et al. (2016), and Ingrisch and Bahn (2018) for analyzing the resilience of systems (using a physical model of soils, and computational and statistical models of industries, and ecosystems, respectively) and system components to different exogenous and endogenous shocks. These characteristics could be applied to any time series data set that measures a change in variable functionality after a disturbance. In fact, when it comes to identifying robust, replicable methods for quantifying resilient behavior, one of the most reliable methods involves the monitored application of shock scenario simulations (e.g., Hodgson et al. 2015, Nimmo et al. 2015, Bitterman and Bennett 2016, Todman et al. 2016, Meyer et al. 2018, Schibalski et al. 2018.

Shock scenario application
Shocks are disturbance events that have the capacity to reduce the baseline (i.e., normal functional state) of any or all components within a dynamic system. The process of measuring and analyzing resilience using a shock-response regime is not conceptually complex but requires a systematic approach and sufficient knowledge of the variables involved. When a system or entity is "shocked," its response can be quantified based on the behavior of its constituent outputs over a known time series. In other words, when the average behavior of system components is known, a pronounced deviation in that behavior (as a result of shock application) belies inherent system vulnerabilities (Carpenter et al. 2009, Angeler et al. 2010, Anderies et al. 2013, Choularton et al. 2015, Todman et al. 2016. According to Sagara (2018), there are two primary benefits of incorporating shock-based measurements into the monitoring and evaluation process of a comprehensive resilience assessment. First, shock scenario analysis improves the conceptual understanding of complex relationships among disturbances, critical capacities, and socio-environmental well-being. Second, shocks and stressors pose significant operational threats to development gains; as such, acknowledging and understanding the capacity for efficient hazard responses is a vital step in the assessment of overall resilience for any complex system (Sagara 2018). The shock scenarios that we apply include various durations and intensities of (1) market inflation and (2) canal water supply reduction (Appendix 2).

Modeling resilience
Quantitative resilience assessment methods often involve the use of statistical or computational modeling techniques (Cimellaro et al. 2010, Cumming 2011, Tyler and Moench 2012, Hodgson et al. 2015, Nimmo et al. 2015, Polhill et al. 2015, Bitterman and Bennett 2016, Todman et al. 2016, Meyer et al. 2018, Schibalski et al. 2018). These methods allow for an explicit description of system processes, enabling the user to obtain concrete, replicable data related to the specific vulnerabilities and adaptive capacities of individual variables within a system. Several authors have explored the concept of quantifiable resilience characteristics through the application of system dynamics (SD) models (e.g., Simonovic and Peck 2013, Candy et al. 2015, Gotangco et al. 2016, Herrerra 2017, Herrera and Kopainsky 2020 and physically based models (e.g., Fowler et al. 2003, Cox et al. 2011, Miller-Hooks et al. 2012. However, it can be argued that the dynamic nature of complex socioenvironmental systems is most reliably represented using a coupled physical-SD modeling approach because coupled models are able to incorporate the concrete nature of physical data modeling with the connectivity and feedback flow of SD models. Coupled modeling (i.e., the use of an integrated system of two or more models in which the communication and interchange of information between constituent models is facilitated computationally) with respect to resiliency analysis is still in its developmental infancy; however, there are several authors who have led the way in terms of coupled model applications for hazard vulnerability assessments (e.g., Schibalski et al. 2018). Through the use of a coupled modeling approach within the resilience and stability landscape domains, Bitterman and Bennett (2016) were able to measure select aspects of agroecosystem resilience using a pre-and post-disturbance comparative functionality procedure. We employ a similar, baseline-reference methodology for analyzing resilience, with five important distinctions: First, our methodology was developed in a participatory context, i.e., the resiliency assessment procedure has been devised with the ultimate goal of encouraging uninhibited, non-expert, stakeholder use; therefore, the methods described herein are intentionally user friendly. Second, the methods employed by Bitterman and Bennett (2016) focus directly on system-level resilience with respect to stability landscapes, whereas we attempted to quantify system component variable resilience concretely, with discrete values relating to the variables' transient shock response, as opposed to average basins of behavior. Third, the integrated model that we employed was developed by coupling a stakeholder-built SD model with a biophysical model using the dynamic coupling software Tinamït, which allows the models to exchange information at runtime (Malard et al. 2017). This innovative form of model coupling allows for the improved exploration of complex relationships among various system elements, as well as the resulting behavioral dynamics of the system, while retaining stakeholder values and inputs (Inam et al. 2017a). Fourth, Bitterman and Bennett (2016) performed repeated scenario simulations based on a set of contemporary farm data, whereas we extracted 30 years of historical dynamics and trends to elucidate how system variables have interacted over time and how they are likely to respond to disturbances in the future. Finally, whereas Bitterman and Bennett (2016) were primarily interested in understanding how cross-scale processes within and between social and ecological domains contribute to overall system resilience, we sought to analyze the resilience of specific system variables in a comparative context (i.e., the extent to which certain variables exhibit resilience compared to other variables under identical shock conditions). The coupled model that we employed has been used to show that certain variables are https://www.ecologyandsociety.org/vol26/iss2/art17/ more critical to overall system stability than others (e.g., canal water supply, government subsidies), and that shocks applied to these keystone variables have stronger effects on the system as a whole than those applied to variables with fewer adjacent connections or dynamic feedbacks. This approach is beneficial in that it allows the model user to pinpoint specific, variable-level failures in a system during a shock or disturbance event, thereby facilitating the development and application of tailored damagemitigation and adaptive capacity measures.
Our primary objective involves the application of a dynamically coupled modeling framework for the development of a stakeholder-friendly, replicable methodology to quantify resilience metrics in a dynamic agroecosystem experiencing a range of socio-environmental shocks. This objective includes the use of these metrics for: (1) comparative variable resilience analyses, and (2) identification of potential regime shifts, transformations, and previously unidentified system vulnerabilities. The application of a dynamically coupled P-GBSDM (physicalgroup built system dynamics model) to quantify socioenvironmental resilience is particularly unique, as is the choice to analyze specific variables within a complex system instead of solely assessing the system itself (e.g., Bitterman and Bennett 2016, Todman et al. 2016). The coupled model used here contributes feedbacks and incorporates complex variable linkages that other models are not able to produce reliably, and the incorporation of the five metrics allows for the notable malleability and adaptability of the resilience assessment procedure for other socio-environmental systems or variables (Appendix 3).

Study site
Rechna Doab is a sub-watershed located in the Indus Plain of central-northeastern Pakistan. The study area lies in a region defined by the latitudinal range 30° 32' N to 31° 08' N, and the longitudinal range 72° 14' E to 71° 49' E. The area of interest covers approximately 732.50 km² and was divided into 215 discrete polygons ( Fig. 1), each with its own unique topology, agricultural divisions, and soil composition. The Rechna Doab ("two waters") basin lies just above the confluence of the Ravi and Chenab Rivers and sits within the Haveli Canal command area (Fig. 2). Approximately 30% of the potentially cultivatable land in the Rechna Doab watershed is presently unexploited because of high soil salinity. The culture and economy of this region are highly dependent on agriculture, and many inhabitants' livelihoods are affected directly by socio-environmental change as a result of climatic or socioeconomic disturbances (Inam et al. 2017b).

Resilience metrics
To determine the degree of resilience exhibited by each variable in each unique shock scenario, five metrics, each describing a unique feature of a resilient response to disturbance, were applied to the normalized data for each intensity and duration combination. The five metrics chosen were based on metrics used by several authors (e.g., Cimellaro et al. 2010, Todman et al. 2016 to quantify resilience based on functional response curves. Reproduced from Inam et al. (2017b) with permission. https://www.ecologyandsociety.org/vol26/iss2/art17/ Several metrics can be used to assess the response of a particular variable in a resiliency context when analyzing a functional response curve, i.e., the function that computationally represents the outputs of a system under different shock scenarios. We applied the following five metrics and their associated equations to the normalized data sets obtained after the shock simulation procedure (Fig. 3): (1) time to baseline-level recovery (Rt), (2) rate of return to baseline (Rr), (3) degree of return to baseline (Rd), (4) overall post-disturbance perturbation (Rp), and (5) corrective impact of disturbance (Rci). We chose these five metrics because of their combined potential for accurately describing the resilience of the study variables based on a transient shock response. Considered individually, these metrics only give a partial explanation of the overall resilience of a shocked entity; however, when examined concurrently, these five metrics describe three important aspects of a completely resilient shock-response. First, these metrics demonstrate the capacity of a variable to withstand and resist stress. Second, they indicate the efficiency with which a variable can recover from disturbance. Third, they account for the very real possibility that a variable will not return to a predisturbance functional equilibrium, allowing for the recognition of potential regime shifts and transformations in the variables of interest.

Fig. 3.
Theoretical shock-response curve and the data boundaries determining each of the five resiliency (R) metrics. Rt = return time to baseline-level recovery, Rr = postdisturbance rate of return to baseline, Rd = degree of return to baseline, Rp = overall post-disturbance perturbation, and Rci = corrective impact of disturbance on system functionality, accounting for an overshooting of the response after the first return to baseline.

Return time, or time to baseline recovery (Rt):
The time to baseline recovery is the amount of time it takes a system or variable of interest to return to a pre-disturbance state of functionality after a shock event. The return time is similar to Holling's (1996) original definition of engineering resilience; this metric quantifies the length of the transient response period of the observed function. If a variable exhibits a relatively high level of resilience in response to the disturbance scenario, then we would expect to see a relatively low time of return to the baseline functionality level (Fig. 4).

Rate of return (Rr):
A resilient system or variable will return to a stable level of functioning more quickly (at a faster pace or steeper gradient) than one that is not resilient. Rate of return is a combined measure of the return time and the magnitude of the impact of the function during the transient response ( Fig. 4; Cimellaro et al. 2010, Hodgson et al. 2015, Todman et al. 2016. Fig. 4. Five dimensions of resilience. In the equations, t 0 is the initial time measurement at the beginning of the simulation period (i.e., t 0 = 0); t r is the time measurement for the point at which the functionality curve returns to baseline postdisturbance; S i is the functional state of the system at maximum shock impact; S r is the functional state of the system after 30 years of simulations, i.e., the functional value at time t = 60 seasons, and may be equivalent to S 0 if the final state of the system is equal to that of the baseline case (i.e., S 0 and S r = 1); t i is the time (x-value) at maximum impact; S 0 is the functional state of the system before the first moment of disturbance (i.e., S 0 = 1); and t r2 is the time of second baseline return in the case of overshooting. Adapted from metrics suggested by Cimellaro et al. (2010), Todman et al. (2016), and Ingrisch and Bahn (2018).
Functional degree of return (Rd): Variables that exhibited a resilient response (with respect to the baseline functional state) were represented by a curve (i.e., data set) that returned to a stable level of functionality closer to that of the reference level. To be more explicit, the degree of return is a measure of the extent to which the observed function comes back to a prescribed reference level; this reference level could be the level of a baseline function before the disturbance (as was determined in this study), or the level of a completely controlled, hypothetical system. The degree of return was measured as the difference between the baseline case (i.e., unshocked, pre-disturbance state) functionality value and the final output value after 30-year simulations ( Fig. 4; Todman et al. 2016).

Post-disturbance perturbation (Rp):
The perturbation was measured using the area above the output response curve but below the baseline case data line. Using this metric, a more resilient system produces a smaller area between the functional response curve and the baseline boundary (Cimellaro et al. 2010, Todman et al. 2016). If the shocked variable never returned to the baseline state of functioning, then the vertical boundary for the perturbation metric was drawn at the point where the variable settled into a new equilibrium. For our purposes, a new https://www.ecologyandsociety.org/vol26/iss2/art17/ equilibrium was defined as the occurrence (post-maximum perturbation) of identical variable outputs (to the nearest 100th degree) for at least two consecutive years (i.e., four seasons). If a new equilibrium was never reached, the boundary was drawn at the end of the viewing window, i.e., after 30-year simulations (Fig.  4). This methodology was developed with the presupposition that the transient changes in function as a result of variable disturbance were undesirable; as such, the most resilient response, as described by the perturbation metric, was that which produced a function unperturbed by disturbance, i.e., perturbation = 0. With respect to Rp, if there was any loss of variable functionality, the area above the functional response output curve was negative, indicating a cumulative loss in function when the system was shocked (Fig. 4).
Corrective impact (Rci): This metric accounts for the potential overshooting of a variable upon return to baseline after disturbance. Rci was calculated as any area above the functional response baseline curve in the event of a post-disturbance increase in functional behavior (Fig. 4;Ingrisch 2018, Yeung andRichardson 2018). Upon initial inspection, a high Rci may appear to be a constructive response to a shock event; however, this response could also bely a systemic inefficiency if it means that a limited resource is being used more quickly. It is important to keep these complex feedback dynamics in mind when analyzing socio-environmental data for resilient characteristics.
The degree of shock resistance exhibited by a variable influences the Rr and Rp of a response data set; thus, these two metrics sufficiently embody the key consequences and outcomes of disturbance resistance. Rt most consistently coincides with system or variable resilience when all five metrics are taken into account separately; however, Rr, Rp, and Rci provide a more reliable measure of overall variable resilience (as opposed to calculating Rr alone) because they make use of all of the available data, rather than a single point. They also negate the need to identify an additional fragility metric for measuring system vulnerability (Todman et al. 2016). Rd rounds out the set of five metrics by accounting for the very real possibility of a shocked variable not returning to a state of pre-disturbance functionality.

Simulation and analysis procedure
Links between the shock nodes in the coupled P-GBSDM and each of the three study variables (farm income, water-table depth, and crop revenue) were initially tested manually using Vensim software (Ventana Systems, Harvard, USA) to ensure feedbacks and connecting loops between adjacent variables were sound and reasonable. Using this manual testing method, we first noticed the trend indicating that severe reduction in precipitation coupled with soil salinity and no policy changes resulted in canal supply reduction. The variables present in the Vensim diagram were initially programmed with their own equations and values based on the stakeholder-built causal loop diagrams, as well as up-todate socioeconomic data from the government of Pakistan. To improve the ease with which data could be transferred between the Vensim model and the Tinamit programming interface, shock switches were incorporated into Vensim, allowing for the manual adjustment of shock durations and intensities. The intensity switches were created by establishing new dimensionless constants, which were subsequently written into the equations for the system components they were modifying. The duration switches were created as new constants in the Vensim model with seasonal time units. The initial equations assigned to each shock factor were altered to include if-then-else statements, accounting for the changes experienced by the variables when intensity and duration values were modified. The original, unaltered equations for the shock variables were written back into the modified equations to account for the baseline case state of the variable; these loops also accounted for the times when the shock needed to be strategically "shut off," i.e., intensity = 1 (Fig. 5). Shock and variable connections were tested manually in Vensim for a second time to ensure the accuracy of the modified equations and to test the links and feedbacks. Once the shock switches were confidently assessed in Vensim, the new duration and intensity metrics were coded into the Tinamit programming interface using a simple Python script.
To streamline the simulation and data procuring process, four Python codes were written to interface with the Tinamit package for automatic data modification. The first code involved the intensity and duration shock switches initially established in Vensim. When run using the P-GBSDM, this script creates unique CSV files that contain the shocked data outputs for each simulation of the switch scenarios. The second code was written for the normalization of the shocked data to the baseline case state of the interest variables, which allows the CSV files created using the first code to be read, organized, and subsequently converted to baseline case normalized data frames. All data points were included for the full 30-year time range across all 215 polygons. Time, as a unique metric standardized across all variables and shock types, was not normalized to a baseline case scenario but was instead recorded as a raw data figure in each run.
Once the data were normalized to the baseline case scenario for each variable, a third code was written for the procurement of the five resilience metric outputs for each variable and shock scenario combination (32 files in total). This code assigned the five metric equations to each polygonal data series in the variable set for each shock combination; in other words, each of the 32 baselinenormalized shock files (containing data for all three study variables in each of the 215 polygons over 30 years) was modified by the five resiliency metric equations, and each of the three variable data sets in each file received five new output values pertaining to the resilience metrics for that polygon. The outputs were then assessed for their comparative levels of resilience based on ideal values. The ideal metric values for a perfectly resilient system are Rt = 0, Rr > 0, Rd = 0, Rp = 0, and Rci ≥ 0. For a perfectly shock-resistant system Rci would be zero, but in the case that a variable is not perfectly robust (not resistant to shock damage, which is likely to be the case), then a high Rci is ideal. These metric values are consistent with our definition that a perfectly resilient result will exhibit similar behavior patterns to those of a system that has not been shocked, i.e., a pre-disturbance state. However, exceptionally large values for Rd or Rci are likely indicators that a regime shift has taken place; a regime shift may have positive or negative consequences based on the specific shock conditions and variables involved. In contrast, a notably large value for the Rt coupled with a large value for Rd may indicate that the system has fallen out of functionality altogether and may https://www.ecologyandsociety.org/vol26/iss2/art17/ be irreparably damaged. Once the final simulation had been run, 32 files containing the five resilience metric values for each variable type and shock combination were available for further analysis (Fig. 6).

RESULTS
As predicted, the three variables assessed each exhibited unique reactions to shocks of varying duration and intensity. We provide select examples of these outputs based on the normalized time series responses of each variable (i.e., the curves used to quantify the resilience metrics) to two different shock types (market inflation × 10 for 10 years and 50% canal supply reduction for 5 years) for the lower (Figs. 7 and 8), middle (Figs. 9 and 10), and upper (Figs. 11 and 12) watershed regions. These shock scenarios were chosen because of their realistic probability of real-world occurrence and because their outputs are representative of the overall variable behaviors for the respective shock types.
Farm income exhibited high perturbation values and long return times, even under the least severe inflation scenario (×2 for 1 year). Water-table depth and crop revenue began exhibiting vulnerabilities under an inflation factor of ×5 for 5 years; under these shock conditions (Figs. 7,9,and 11), all three variables exhibited very high Rp values, mostly from upper-watershed polygons (Figs. 7 and 8), indicating that they were all pushed very far out of their normal functional patterns before returning to baseline. At this point, there were also very high instances of 61 Rt (indicating that the variable did not return to the baseline case functional state after 60 seasons) and 0 Rci for water-table depth (indicating that water-table depth performed very poorly overall for this trial). With an inflation factor of ×10, the pattern emerged of crop revenue exhibiting high Rp values and extremely high Rci values, with a very large degree of return; these results indicate high fragility and instability in these variables for this shock type. Water-table depth began to perform the best of the three variables under this shock scenario, whereas the two socioeconomic variables became increasingly vulnerable to the inflation shock type at higher intensities. With an inflation intensity of ×10 for 20 years, the farm income variable lost the capacity to return to baseline after shock application, i.e., farm income was never able to recover from inflation of that intensity and duration. With respect to the ideal resilience values, farm income performed the worst under inflation shock conditions. Crop revenue was the most inconsistent variable, and water-table depth was the most stable variable, although the latter also nominally faltered as duration increased. All five metrics of water-table depth increased as the duration of the shock increased; this trend included the Rci value, indicating that the shock induced erratic variable behavior and that shock duration has a greater effect on water-table depth than intensity. Rci decreased notably for crop revenue as the shock duration increased.
Under canal supply shock conditions (Figs. 8, 10, and 12), high Rp and Rci values frequently occurred at the head of the watershed, indicating that these polygons are both highly adept at recovering from a shock but also highly unstable throughout all canal supply shock runs. All "lowest" metric values actually decreased from run S2,10,01 to run S2,10,10 (i.e., as shock duration increased from 1 to 10 years for a reduction of 10% in canal water supply), and most, but not all, of the metrics on the high end increased from run S2,10,01 to S2,10,10; this result indicates higher rates of fluctuation or instability in functionality Ecology and Society 26(2): 17 https://www.ecologyandsociety.org/vol26/iss2/art17/ for variables as the shock duration increases. The canal supply shock induced large values of Rt and Rp for farm income compared to the other two variables, indicating that farm income experienced the most difficulty recovering to a pre-disturbance state after the shock event. Rci values for crop revenue were also very high compared to the other two variables, and crop revenue showed comparatively high values for rate of return and degree of return across all intensities and durations of this shock type. These high values indicate that crop revenue is generally resilient under these shock conditions, but some polygons are unpredictably unstable and may have experienced functional regime shifts. All low values for all three variables in Rci were 0 for this shock, which was not the case for the inflation shock at any intensity; this result indicates that the variables were better      Table 1. Regional resilience metric outputs. S1 = market inflation shock of ×10 sustained for 10 years, S2 = canal supply shock of −50% sustained for 5 years, Rt = time to baseline-level recovery, Rr = rate of return to baseline, Rd = degree of return to baseline, Rp = overall post-disturbance perturbation, and Rci = corrective impact of disturbance. able to maintain a level of homeostasis under these shock conditions than they were for the inflation shock. All metric values for water-table depth increased from run S2,50,01 to S2,50,20 (i.e., 50% reduction in canal water supply sustained from 1 to 20 years); this trend did not occur for crop revenue and farm income, which both saw a decrease in Rci from run S2,50,01 to S2,50,20, and farm income saw negative changes in all metrics from run S2,50,01 to S2,50,20 (i.e., increased Rp, Rt, and Rd, and decreased Rci and Rr). Rci decreased for all variables, and Rp increased for all variables, from run S2,90,01 to S2,90,20 (i. e., 90% reduction in canal water supply sustained for 1 to 20 years). Thus, an increase in intensity seems to decrease the capacity of the variables to commit sufficient corrective behavior. Crop revenue and farm income saw negative changes for every single metric between runs S2,90,01 and S2,90,20 (i.e., increased Rt, Rp, Rd, and decreased Rr and Rci.) Water-table depth remained remarkably consistent throughout all canal supply runs, indicating that this variable is highly resilient, especially when considering that this is a physical shock (canal water supply reduction).
We produced watershed-level heat maps (Figs. 13-16) and a regional resilience metric table (Table 1) using a market inflation intensity of ×10 for a duration of 10 years (Figs. 13 and 15) or a canal supply reduction of 50% for 5 years (Figs. 14 and 16). The relative scale shows the precise measurements taken in each polygon (Figs. 13 and 14), whereas the standardized scale allows comparison among watersheds (Figs. 15 and 16). A regional pattern emerged based on the comparative metrics of the three study variables under identical shock conditions. In the upper watershed, water-table depth displayed the most resilient behavior under both socioeconomic and physical shock conditions. Farm income was unable to return fully to baseline in either shock scenario, and both farm income and crop revenue displayed massive perturbation values for the market inflation shock type. Farm income fared the worst under canal supply shock conditions, and crop revenue displayed signs of overcompensation under the canal supply shock, indicating a potential functional regime shift; both socioeconomic variables showed a comparable level of resilience (or lack thereof) under the market inflation shock. Water-table depth was extremely stable throughout both shock trials in the upper watershed. These trends held true for the middle-watershed polygons as well, with farm income performing poorly under market inflation shock conditions, and crop revenue exhibiting strong over-corrective patterns under canal supply shock conditions. It is not until the lower watershed is examined that a change in water-table depth patterns can be identified. In the lower watershed (i.e., polygons farthest from the head or source of the watershed) water-table depth begins exhibiting higher values for perturbation, degree of return, and return time, all indicating a general loss in resilience for the water-table depth variable for both shock types. In the lower watershed, farm income continued to perform poorly under both shock conditions, whereas crop revenue actually exhibited the greatest robustness (resistance to shock influence) for the shocks in this region. Interestingly, the data (Table 1) indicate regional differences between each of the regions from North to South, and the heat maps indicate additional differences between East and West (particularly the southwestern corner of the watershed, where crop revenue exhibits particularly high resilience). These results reinforce the importance of analyzing regional trends from multiple perspectives. The clear difference in resilience of the study variables based on watershed regions and individual polygons is a textbook example of spatial resilience whereby trends and outcomes at different scales both affect, and are affected by, local system resilience (Cumming 2011).

DISCUSSION
Through each shock trial and for each of the watershed regions, water-table depth most consistently aligned with the ideal resilience metric values. This result was expected under socioeconomic shock conditions (i.e., market inflation), but it is notable that water-table depth continued to exhibit the greatest resilience and robustness even during canal supply shock scenarios. One reason for this seemingly inherent robustness is that water-table depth is a "slow" variable, meaning that it reacts less dramatically (at least initially) in response to social-ecological drivers, and also has the capacity to influence "fast" variables (e. g., farm income, crop revenue), which are adjacent in the system and experience the same system-level stresses (Walker et al. 2012).
Our study addresses issues related to internal drivers that are Ecology and Society 26(2): 17 https://www.ecologyandsociety.org/vol26/iss2/art17/ Fig. 13. Heat maps showing levels of the five resilience metrics across the Rechna Doab watershed for market inflation shock of ×10 sustained for 10 years, in relative scale. Rd = degree of return to baseline, Rr = rate of return to baseline, Rci = corrective impact of disturbance, Rt = time to baseline-level recovery, Rp = overall post-disturbance perturbation. incorporated on a variable-level scale but, owing to the dynamic nature of the coupled model, act as system-level drivers affecting most variables in the system. This function is evidenced by the time-series plots showing the marked responses of each variable type to each different shock type. However, water-table depth showed the greatest signs of resilience loss in the lower watershed regions. This result makes intuitive sense because the lower regions are farthest from the source stream and least likely to receive water in amounts copious enough to reserve for times of scarcity. The extreme robustness (i.e., shock resistance) exhibited by water-table depth in the upper and middle regions of the watershed can be explained by the capacity of these regions to exercise water reservation practices based on their more advantageous location in the basin closer to the head of the watershed. These results indicate areas for improvement in watershed-level supply allocation, irrigation infrastructure, and water banking policy; the identification of these specific sectors requiring improvements are supported by the findings of Inam et al. (2015Inam et al. ( , 2017b in the Rechna Doab basin. Interestingly, farm income showed the lowest capacity for resilience and shock resistance of any of the three study variables, regardless of shock-type, duration, or intensity. This result indicates that farm income is itself a very fragile variable, influenced by watershed-level disturbance events of both socioeconomic and biophysical origin. Likewise, crop revenue, the other socioeconomic variable, exhibited extreme fluctuation in perturbation values while also maintaining consistently large corrective impact values. This erratic behavior indicates that crop revenue is the variable most likely to experience regime shifts in times of stress. Transition to a different baseline level of functionality can be an indicator of extreme functionality loss and also tremendous adaptive capacity; whether it is the former or the latter depends on the response of closely adjacent variables in the system. For a socioeconomic variable, however, it is likely the former. That is to say, if there were a notable increase in crop revenue as the result of high inflation or a drop in water supply, this revenue increase would be quite unsustainable over a long time period. The spatial resilience exhibited by the study variables (Figs. 7-14, Table 1) has roots in several socioeconomic and ecological processes, most notably, the unequal distribution of water resources from upper to lower watershed polygons. Not only do upper watershed farmers have more reliable access to fresh https://www.ecologyandsociety.org/vol26/iss2/art17/ Fig. 14. Heat maps showing levels of the five resilience metrics across the Rechna Doab watershed for canal supply shock of −50% sustained for 5 years, in relative scale. Rd = degree of return to baseline, Rr = rate of return to baseline, Rci = corrective impact of disturbance, Rt = time to baseline-level recovery, Rp = overall post-disturbance perturbation.
water than their downstream counterparts, but government subsidies have incited farmers of all regions to increase cropping intensities, leading to unsustainable drawdown of groundwater resources, especially in the middle watershed (Inam et al. 2015). The depletion of groundwater resources exacerbates the fragility of variables in the middle and lower watershed regions by reducing the adaptive capacities of these regions in times of systemic stress. The inability to adapt effectively to changing conditions is reflected in the observed spatial trends. Keeping these trends of resource allocation in mind, it is highly likely that the observed resilience metrics would hold true (on average) for the Rechna Doab at both finer and courser scales. It is worth noting that even within the current polygonal structure, there is variation in the resilience of individual farms, and the trend of reduced resilience from upper to lower watershed is not linear; however, the methodology described herein could easily be used at such a scale as to determine the resilience of individual agricultural operations, which is useful information.
The P-GBSDM was built with multiple feedback loops and complex socio-environmental relationships linking the variables in the system. It is therefore unsurprising that the variables that reacted strongly to one shock type also experienced changes when exposed to another, even starkly different, type. This innovative resilience assessment scheme will allow stakeholders and model users to understand better the unique vulnerabilities and adaptive capacities of certain variables in dynamic agroecosystems. The stakeholder knowledge used to develop the GBSDM half of the dynamically coupled model was absolutely critical to the understanding of the study system and its constituent variables and feedbacks. The methodology we developed was designed to test the capacity of the dynamically coupled model to produce realistic scenarios that can be used in real-world resilience assessments. Application of the methodology by stakeholders is a continuation of our research and will be explored in detail in future publications. The methods and results described herein are directly applicable and relevant to classic resilience theory whereby resilience is understood as the capacity of a system to absorb or withstand perturbations, disturbances, and various stressors such that the system remains within the same regime (or stability landscape), essentially maintaining its structure and function (Holling 1973, Gunderson and Holling 2002, Walker et al. 2004. The use of a stakeholder-informed, dynamically coupled model for variablelevel resilience quantification provides a valuable contribution to Ecology and Society 26(2): 17 https://www.ecologyandsociety.org/vol26/iss2/art17/ Fig. 15. Heat maps showing levels of the five resilience metrics across the Rechna Doab watershed for market inflation shock of ×10 sustained for 10 years, in standardized scale. Rd = degree of return to baseline, Rr = rate of return to baseline, Rci = corrective impact of disturbance, Rt = time to baseline-level recovery, Rp = overall post-disturbance perturbation.

CONCLUSION
Using the integrated P-GBSDM, discrete variable-level shock scenarios were simulated to determine the dynamic response patterns of farm income, water-table depth, and total crop revenue in each unique regional polygon of the Rechna Doab basin. Following shock-scenario simulations, the output data from each variable were analyzed using five metrics describing a resilient response to disturbance. The five resiliency metric outputs were subsequently analyzed for each of the three interest variables under identical shock conditions. Each polygon in the watershed was assessed separately, each receiving a comprehensive analysis of the comparative resilience of the study variables according to the five calculated resiliency metrics. A comprehensive assessment relating to regional and watershed-level resilience was conducted based on the outcomes of the metric analysis for each variable, under each shock condition, in each unique polygon. Our results show that the methodology allows the user to examine the intricate differences and discrepancies between variable reactions to stress for both socioeconomic and environmental-physical variables and shock scenarios. Because of the realistic outputs provided by the dynamically coupled model, this approach for variable-level resilience quantification has some beneficial realworld applications in the spheres of disaster mitigation policy, vulnerability and adaptive capacity assessments, and long-term risk analyses.
The practical limitations of our methods could present some challenges for the widespread application of the approach to other systems. These limitations include potential deficiencies in the data required to build an accurate model of a chosen study system, as well as any difficulties related to the inclusion of stakeholderdefined variables and processes in a coupled model, which relies https://www.ecologyandsociety.org/vol26/iss2/art17/ Fig. 16. Heat maps showing levels of the five resilience metrics across the Rechna Doab watershed for canal supply shock of −50% sustained for 5 years, in standardized scale. Rd = degree of return to baseline, Rr = rate of return to baseline, Rci = corrective impact of disturbance, Rt = time to baseline-level recovery, Rp = overall post-disturbance perturbation.
heavily on accurate feedbacks. With these limitations in mind, several elements of our study could provide the basis for further scientific investigation. For example, we explored shock scenarios from a discrete perspective, i.e., simultaneous disturbances were not taken into account. We recommend that future research, applying similar methods, should be conducted using compound disturbance scenarios, in different regions, climates, and with additional focus variables.
Responses to this article can be read online at: https://www.ecologyandsociety.org/issues/responses. php/12354