Water clarity and temperature effects on walleye safe harvest: an empirical test of the safe operating space concept

Successful management of natural resources requires local action that adapts to larger-scale environmental changes in order to maintain populations within the safe operating space (SOS) of acceptable conditions. Here, we identify the boundaries of the SOS for a managed freshwater fishery in the first empirical test of the SOS concept applied to management of harvested resources. Walleye (Sander vitreus) are popular sport fish with declining populations in many North American lakes, and understanding the causes of and responding to these changes is a high priority for fisheries management. We evaluated the role of changing water clarity and temperature in the decline of a high-profile walleye population in Mille Lacs, Minnesota, USA, and estimated safe harvest under changing conditions from 1987 to 2017. Thermal–optical habitat area (TOHA)—the proportion of lake area in which the optimal thermal and optical conditions for walleye overlap—was estimated using a thermodynamic simulation model of daily water temperatures and light conditions. We then used a SOS model to analyze how walleye carrying capacity and safe harvest relate to walleye thermal–optical habitat. Thermal–optical habitat area varied annually and declined over time due to increased water clarity, and maximum safe harvest estimated by the SOS model varied by nearly an order of magnitude. Maximum safe harvest levels of walleye declined with declining TOHA. Walleye harvest exceeded safe harvest estimated by the SOS model in 16 out of the 30 yr of our dataset, and walleye abundance declined following 14 of those years, suggesting that walleye harvest should be managed to accommodate changing habitat conditions. By quantifying harvest trade-offs associated with loss of walleye habitat, this study provides a framework for managing walleye in the context of ecosystem change.


INTRODUCTION
Freshwater resources are threatened by environmental change, including climate change, land-use change, invasive species, and harvest (Carpenter et al. 2011). Ecosystem responses to these drivers of global change may be non-linear, responding gradually until a tipping point or threshold is reached from which recovery can be difficult or impossible (Scheffer et al. 2015, Carpenter et al. 2017). The safe operating space (SOS) concept for managing ecosystems identifies the boundaries of acceptable conditions that are defined by interactions between continental or global scale drivers and local management (Scheffer et al. 2015). Theory suggests that by following the boundaries of the SOS, local management actions can be adjusted in response to environmental change to maintain ecosystem services and increase resilience. However, empirical tests of this important resilience concept on scales relevant to natural resource decision-making are lacking (Carpenter et al. 2017). Here, we present an empirical test of the SOS concept based on walleye (Sander vitreus), an economically and ecologically important freshwater fish species. We estimate the effect of changing habitat on sustainable harvest of walleye, identify the boundaries of the SOS for walleye harvest as a function of habitat, and show that safe harvest levels based on the SOS differ from those based solely on traditional fisheries models.
Walleye prefer low water clarity (Ryder 1977) and cool temperatures (Christie and Regier 1988). Walleye are low-light specialists due to a specialized retinal structure known as the tapetum lucidum that develops during the first year of life and allows them to successfully forage in dim conditions Anctil 1977, Vandenbyllaardt et al. 1991). For most lakes, Secchi depths of 2-3 m are optimal for walleye (Lester et al. 2004), and increasing clarity above this range reduces optical habitat area. Thermal habitat area can be positively affected by water temperature due to increased growing season duration ( Fig. 1), or negatively affected if temperatures exceed upper thermal limits. A lake's clarity, temperature, and bathymetry determine its thermaloptical habitat area (TOHA), that is, the area of a lake in which optical and thermal conditions for walleye overlap. A lake's TOHA is positively related to walleye production and catch rates at broad spatial scales (Lester et al. 2004, Tunney et al. 2018. Increasing water clarity and warming temperatures are associated with declining walleye and increasing Centrarchid (sunfishes and black bass) populations in lakes throughout North America (Robillard and Fox 2006, Hansen et al. 2015. It is assumed that temporal trends in TOHA will influence walleye carrying capacity and yield (Lester et al. 2004), but to date, empirical tests of the effects of changing thermal-optical habitat on walleye populations are lacking (but see Chu et al. 2004).
Most fisheries models and stock assessments assume that relationships between stock size and population rates are stationary, and set harvest policies accordingly (Walters 1987). However, climate change, invasive species, harvest, and other stressors can alter productivity, with important implications for sustainable harvest policies (Walters et al. 2008). Sustainable fisheries management in the 21st century must account for the effects of global change (Paukert et al. 2016), although few concrete examples exist of recreational fisheries management systems that explicitly incorporate environmental change. The SOS of a recreational fishery is defined by environmental conditions and local management, and harvest reductions may compensate for habitat loss and prevent population collapse as conditions change (Fig. 1;Carpenter et al. 2017). Under this framework, longterm harvest is increased by adapting annual harvest in response to changing environmental conditions ( Fig. 1).
In this study, we quantify relationships between water clarity and temperature, walleye habitat, and safe harvest. Our study focuses on Mille Lacs, Minnesota, USA, where walleye populations have dramatically declined since the 1990s. Due to the lake's popularity and economic importance, strong social and political pressures exist to restore walleye in Mille Lacs to support previous levels of harvest. However, if ecological changes have altered the productive capacity of the lake, harvest may need to remain low to maintain a sustainable fishery. The objectives of this study were to (1) quantify changes in walleye habitat area due to changing water clarity and temperature, and (2) quantify the effects of habitat area and predators on sustainable walleye harvest levels in an empirical test of the SOS concept.

Study area
Mille Lacs is a large (519 km 2 ), shallow (mean depth = 8.7 m), mesotrophic, polymictic lake located in central Minnesota, USA (46.233, À93.6502). Mille Lacs was historically one of Minnesota's most popular and productive walleye fisheries, but the walleye population and thus harvest has declined since the 1990s ( Fig. 2; Venturelli et al. 2014). The timing of this decline coincides with many changes, including warming temperatures (Fig. 2), changes in the fish community such as increasing smallmouth bass (Micropterus dolomieu) and fluctuating northern pike (Esox lucius) abundance (as measured by gillnet catches in standardized surveys; Fig. 2), the establishment of an Ojibwe tribal fishery in 1997 when treaty rights were reaffirmed (Minnesota v. Mille Lacs Band of Chippewa Indians 1999), and the invasion of zebra mussels in 2005 and spiny water flea in 2009 (MN DNR 2018). Notably, a marked increase in water clarity co-occurred with the onset of walleye declines, with Secchi depth changing from an average of 2.5 m from 1977 to 1996 to 3.5 m from 1997 to 2016 (Fig. 2). Other water quality data are sparse, but suggest that total phosphorus concentrations were higher in the 1970s through early 1990s than during the 2000s (

Walleye data
Walleye population size in Mille Lacs is estimated annually using a statistical catch at age (SCAA) model (Schmalz and Treml 2014). The SCAA model projects population numbers by sex, age, and length using fishery-dependent and fishery-independent data. Fishery-independent data include sex-and age-specific gillnet catch rates, catch rates of age-0 and age-1 walleye from fall electrofishing, and six independent markrecapture population estimates. Outputs from the SCAA model were used as observations in the SOS model (described below). Relevant outputs used here include annual population estimates of age-3 and older walleye and total walleye kill from 1987 to 2017 (Fig. 2). Total kill includes walleye harvested by tribal fisheries, walleye harvested by recreational anglers (estimated from non-uniform probability accessbased creel surveys; Pollock et al. 1994), and walleye killed via hooking mortality in the recreational fishery (estimated from creel surveys and a temperature-dependent statistical model; Reeves and Bruesewitz 2007). For brevity, these three sources of walleye mortality are collectively referred to as harvest throughout. The SCAA safe harvest limit was set at 24% of the biomass of walleye ≥35.6 cm (14 in) in total length from 1997 to 2014 (Minnesota v. Mille Lacs Band of Chippewa Indians 1999). Actual harvest quotas were negotiated based on these SCAA limits as well as additional information, and total walleye kill was always lower than SCAA safe harvest limits during this period (M. Treml, unpublished data). Declining populations have led to increasingly strict walleye harvest regulations in the late 2000s (Schmalz et al. 2011), and recreational harvest has been closed at various times from 2015 to 2017.

Optical, thermal, and thermal-optical habitat area
To estimate changes in walleye habitat area, we quantified optical, thermal, and TOHA for Mille Lacs for each day of 1980-2016. We used a combination of observed data, hydrodynamic modeling, and statistical modeling to reconstruct thermaloptical parameters. Water clarity was measured by Secchi depth (Appendix S1: Table S1) and converted to daily light extinction coefficients using a non-linear hierarchical model (Appendix S1). Daily water temperature profiles were estimated using an open-source hydrodynamic model (General Lake Model v2.2; Hipsey et al. 2019), modified to incorporate daily estimates of light attenuation. The temperature model was calibrated to in situ temperature data using a Nelder-Mead gradient descent algorithm, whereby the overall root-mean-squared error (RMSE) was minimized by altering model parameters (Appendix S1: Table S2; final RMSE of 1.26°C).
Daily water temperature and clarity estimates were combined to calculate daily TOHA following Lester et al. (2004), with minor modifications as described in Appendix S1 (for R code, see Winslow et al. 2017). Thermal habitat area was defined by temperatures for which simulated walleye growth rates were within 50% of the maximum from bioenergetics model simulations (11-25°C; Lester et al. 2004). We calculated thermal habitat area as lake benthic area for which water temperature fell within this range. Mille Lacs does not stratify, so habitat area calculations included the entire lake bottom. Optical habitat area was defined as lake benthic area for which estimated light levels fell between 8 and 68 lux (Ryder 1977, Lester et al. 2004. Thermal habitat area, optical habitat area, and the combination (TOHA) were expressed as the total benthic area satisfying the aforementioned criteria for each day. Daily estimates were summed to calculate annual habitat area estimates. This total annual area estimate was divided by the total  In an unstratified lake that does not usually exceed upper thermal limits for walleye such as Mille Lacs, the entire lake bottom is thermally optimal for each day that temperatures fall within the optimal range, and warming temperatures increase this duration (but can also cause water temperatures to exceed thermal limits for optimal growth on some days). Optical habitat changes daily and seasonally as a result of diurnal and seasonal sun angle and seasonality of water clarity. In Mille Lacs, increasing water clarity reduces optical habitat as conditions become too bright throughout the ❖ www.esajournals.org benthic area times the number of days in the year. Thus, habitat area is expressed as the proportion of the maximum potential habitat area averaged across the entire year (Fig. 1A). Three-year moving averages of annual TOHA were used as inputs to the walleye population model.

Walleye population model and the safe operating space
We assessed the effects of habitat area on walleye abundance and harvest using a previously studied model of fish population dynamics (Carpenter 2002). This model allows for a SOS (Carpenter et al. 2017) that is influenced by harvest and lake conditions. The general form of the population models we considered was x t is the number of age-3 and older walleye in year t, H t is the TOHA (averaged across years t À 2, t À 1, and t), P t is the catch per net night of northern pike and smallmouth bass from gillnet surveys, b i are fitted parameters, F t is the number of walleye killed by harvest (including hooking mortality), and e t is the model residual error assumed to be normally distributed, mean 0, variance r 2 estimated from the data, and uncorrelated over time . Both x t and F t were estimated by the SCAA model using observations from Mille Lacs (Schmalz and Treml 2014). Models were fit by maximum likelihood and compared using AIC (Akaike 1973). Two model forms fit the data relatively well: Both models include a linear autoregressive term (b 1 ) and a quadratic habitat term (b 2 ) analogous to a logistic equation. Model (2a) includes a predator term (b 3 ) with a linear (i.e., Lotka-Volterra) functional response. We considered more complicated functional responses (Walters and Martell 2004) and none fit as well as a linear response. In (2a) and (2b), we write f(x) instead of f(x t , H t , P t , b t ) to simplify notation.
The SOS boundary is the highest possible fishing mortality that still has a positive equilibrium growth rate for given and fixed values of H and P (Carpenter et al. 2017). Solutions to Eq. 3 with e t = 0 provide deterministic equilibrium walleye population sizes.
The edge of the SOS occurs where the flat line F is tangent to the hump-shaped relationship between population size and population growth rate as described by Eq. 3. F s is the maximum walleye mortality, and x s is the corresponding walleye population for specific fixed levels of habitat area and predator biomass. If particular values x s and F s are at the edge of the SOS, then (3) is satisfied and the first derivative of (3) is zero. Thus, x s can be found by solving the first derivative of (3) with respect to x, which is Given estimates of the parameters b i , (4) is solved numerically for x s using the uniroot() package in R (R Core Team 2017). Then, F s can be found by solving (3) using x s : water column for most hours of most days. (B) The safe operating space (SOS) is defined by both habitat and harvest. Environmental change such as increasing water clarity and reducing habitat can push the system out of the SOS (yellow dot), meaning that previously sustainable harvest levels now exceed safe harvest limits and the population collapses (red dot). Harvest reductions can move the system back to the SOS, and harvest will gradually increase to a new equilibrium (dark blue dot). If harvest is gradually decreased as conditions change (dashed orange line), total harvest over the time interval (area under the dashed orange curve) will be much larger than if the population is allowed to collapse (area under the black line). Even though the final rate of annual harvest is the same in the final year, the total harvest across all years is greater in the case where harvest adapts to changing habitat.
( Fig. 1. Continued) ❖ www.esajournals.org summer (July-August-September) air temperatures based on interpolated topoclimatic daily air temperatures: https:// catalog.data.gov/dataset/topowx-topoclimatic-daily-air-temperature-dataset-for-the-conterminous-united-states. Note that zebra mussels were discovered in 2005. Starting from observed time series of the walleye population, TOHA, and predators, we estimated the parameters b i of (2a) and (2b), as well as the standard deviation of model errors, by maximum likelihood. We examined residuals from the model fits to assess that model errors were approximately normally distributed and uncorrelated in time. We then estimated errors of the b i by bootstrapping, using 10,000 random permutations of residuals with replacement (Efron and Tibshirani 1993). We then computed a population of bootstrap estimates of the SOS for specified values of TOHA or predators.
An estimate of the SOS consists of a pair of numbers, F s and x s , representing the maximum fishing mortality with positive population growth and the corresponding population size, respectively. We estimated the boundaries of the SOS by fitting the model to observed time series, solving the model for the upper bound of the SOS, and estimating errors of model parameters and the SOS by nonparametric bootstrapping. 90% confidence intervals (CIs) for F s and x s were approximated by estimating the 95% and 5% quantiles from values estimated from the 10,000 bootstrapped parameter sets. In a small number of cases, bootstrapped parameter sets yielded negative values for the SOS. In these cases, we set the SOS to zero. Maximum safe harvest at the SOS (F s ) and its CI were compared to total walleye kill in Mille Lacs to assess whether the fishery was overharvested based on the SOS estimates. These differences in harvest were also compared to change in estimated population size in the following year to examine whether years in which walleye were overharvested were followed by population declines.

RESULTS
Thermal-optical habitat area for walleye in Mille Lacs declined over time (Fig. 3). Optical habitat area was most widespread in 1988 and 1999, when 15% of potential habitat area fell within preferred optical conditions. By contrast, optical habitat area was most restricted in 1997, when only 3% of habitat area was optically suitable. Thermal habitat area remained relatively constant over the same time period (Fig. 3). Thermal-optical habitat area declined over time, driven by changes in optical habitat area.
Thermal-optical habitat area was most restricted in 1997, and most available in 1988, 1993, and 1999 when over 7% of potential habitat area fell within preferred light and temperature ranges.
The walleye model without predators fit the data slightly better than a model containing predator abundance (AIC [predators] = 97.49; AIC [no predators] = 95.69). Here, we present results for the model without predators; see Appendix S2 for results of the model including predators. The walleye population model was able to recreate population trends (Appendix S2: Fig. S1). The autoregressive parameter and the effect of TOHA on carrying capacity were both positive (Appendix S2: Table S1). Maximum safe harvest at the SOS boundary (F s ) was positively and non-linearly related to TOHA (Fig. 4). Population size at the SOS boundary (x s ) was also positively related to TOHA (Appendix S2: Fig. S2). Small changes in habitat led to relatively large changes in safe harvest level-for example, when TOHA was 6.5% of lake area, maximum safe harvest was approximately 430,000 walleye. If habitat area declined by half (TOHA = 3.2%), maximum safe harvest dropped to less than a quarter of peak values to under 100,000 walleye. Actual harvest exceeded the maximum likelihood predicted safe harvest in 16 of 30 yr, and exceeded upper bootstrapped 90% CIs in 3 yr (1992, 1996, and 1998). In recent years (2013-2016), harvest has been well below estimated safe harvest levels. Walleye abundance declined following 14 of the 16 yr in which actual harvest exceeded the maximum safe harvest estimated by the SOS model (Fig. 5). Population increases were observed in 10 yr out of the time series, and in 8 of these years, the previous years' harvest fell below estimated safe levels.

DISCUSSION
The SOS for fisheries defines the range of conditions that maintains fish biomass and harvest at acceptable levels even as the environment changes (Carpenter et al. 2017). We identified the SOS for walleye populations in Mille Lacs as a function of habitat (TOHA) and walleye harvest in the first empirical test of the SOS concept applied to the management of harvested resources of which we are aware. Walleye habitat area in Mille Lacs declined over the past several decades, with important implications for fisheries management. The historical range of TOHA resulted in estimated safe walleye harvest levels that varied by nearly an order of magnitude. Walleye mortality in Mille Lacs exceeded mean safe levels based on TOHA in about half of the Fig. 4. Walleye harvest (in 100,000s of walleye) as a function of thermal-optical habitat area. Maximum safe harvest at the safe operating space across a theoretical gradient of habitat values estimated from maximum likelihood parameter estimates (black line) and 90% confidence intervals (CIs) estimated from 10,000 bootstrapped parameter sets (gray band) based on the SOS model. Walleye kill from 1987 to 2016 is also shown as a function of estimated thermal-optical habitat in each year (three-year moving averages). Colored points represent walleye harvest relative to the maximum safe harvest generated from the SOS model (red, harvest exceeded safe harvest; green, harvest was less than safe harvest). Filled circles represent years in which actual harvest fell outside the 90% CIs of safe harvest estimated by the SOS model. Fig. 5. Walleye harvest relative to estimated maximum safe harvest (x-axis) versus change in walleye abundance the following year (y-axis). Quadrants indicate years where the population was overharvested and declined (red), overharvested and increased (yellow), harvested within safe limits and increased (green), and harvested within safe limits and decreased (blue). Years are harvest years. Lines are 90% confidence intervals (CIs) of the x-axis values, representing uncertainty in safe harvest levels. Filled circles represent years in which actual harvest fell outside (either above or below) the 90% CIs of safe harvest estimated by the SOS model. past 30 yr, though uncertainty surrounding these safe harvest estimates is high. In the majority of cases where harvest exceeded estimated safe levels, the walleye population declined in the following year. Similarly, population increases were more likely to occur following years in which harvest fell within estimated safe levels. The major exception to these patterns occurred following harvest year 1990, when the walleye population increased in 1991 due to much higher than average recruitment to the fishery in spite of overharvest. The walleye population also declined in five years for which harvest fell within safe levels (blue quadrant in Fig. 5), although the magnitude of these decreases was small, indicating other sources of mortality not accounted for in our analysis. While TOHA and harvest cannot explain all variation in walleye populations, accounting for TOHA in harvest decisions makes harvest less risky and more likely to stay within the SOS.
In response to observed walleye population declines, walleye mortality has been well below estimated safe harvest levels in recent years (2013)(2014)(2015)(2016). Based on our population model incorporating TOHA, such precautionary management should allow the walleye population to increase, depending on future habitat availability. Note that maximum safe harvest levels and population size at SOS are equilibrium values for a fixed TOHA value (x-axis in Fig. 4) and do not account for annual dynamics of TOHA. Furthermore, our model does not account for all potential drivers of walleye population abundance, such as population, food web, or ecosystem productivity changes associated with invasive zebra mussels  or spiny water flea (Strecker et al. 2011) beyond what is captured by changing water clarity, and therefore may overestimate the walleye production that could currently be supported. Continued precautionary management and monitoring will help elucidate whether changes in TOHA are the main driver of walleye declines.
Increasing water clarity was the main driver of changing TOHA. Several mechanisms could account for this increase. Total phosphorus declined from 1992 to 2005-2013, and improvements to septic systems and land use around the lake may have played a role (Heiskary and Egge 2016). Zebra mussels increase water clarity (Higgins and Vander Zanden 2010) and hence optical habitat area (Geisler et al. 2016), but zebra mussels were discovered in Mille Lacs in 2005 and cannot explain the observed increase in water clarity in the 1990s. Water clarity peaked in 2013, suggesting that zebra mussels may have further increased water clarity once they established. Thermal habitat area was relatively unaffected by increasing temperatures; Mille Lacs water temperatures exceeded 25°C in only 12 days of the 30 yr modeled here. However, changes in thermal habitat as defined by optimal growth conditions may not correlate with changes in fish populations. For example, increasing temperatures are associated with walleye declines in inland lakes (Robillard andFox 2006, Hansen et al. 2017), despite increased walleye thermal habitat area in most lakes as climate warms (Fang et al. 2004). As the climate continues to warm, the number of days per year exceeding walleye thermal tolerance will increase and may negatively influence survival and growth of walleye in the future.
The response of TOHA to changing conditions depends on lake characteristics including morphometry, historical baseline, and stratification, although on average across all lakes TOHA is optimized at Secchi depths of 2 m (Lester et al. 2004). Mille Lacs is shallow and well-mixed lake, and historic Secchi depths were around the optimum value of 2 m. These factors increase sensitivity to increased water clarity and likelihood of impacts to walleye (Geisler et al. 2016). The trajectory of Mille Lacs in terms of ecosystem and fish community changes appears to be similar to what has been documented in Lake Oneida (New York, USA), another high-profile walleye fishery which has undergone dramatic changes in recent decades , suggesting that these dynamics may not be unique. Still, other lakes with more complex bathymetries and different trophic status may respond differently to changing conditions. Lakes with more TOHA support higher walleye biomass and harvest (Christie and Regier 1988, Lester et al. 2004, Tunney et al. 2018, although changes in TOHA over time have rarely been quantified or linked directly to walleye abundance (but see Chu et al. 2004, Jones et al. 2006. Water clarity and TOHA affect walleye populations through a number of pathways. Walleye are more active in low-light conditions (Ryder 1977), and increased water clarity is associated with a shift from feeding during the day in turbid water to crepuscular or nocturnal feeding in clear water . Recent research has also demonstrated that walleye populations in lakes with low water clarity can access multiple prey sources and achieve higher biomass compared to lakes with high water clarity (Tunney et al. 2018). As water clarity increases and TOHA decreases, walleye may be restricted to offshore and deepwater habitats throughout most daylight hours. If energy resources in these habitats are limiting (as is likely to be the case in a system invaded by zebra mussels, Higgins and Vander Zanden 2010), populations are likely to be negatively affected by increasing water clarity. Understanding the mechanisms through which changes in water clarity affects walleye behavior, growth, and population dynamics and how such responses differ among lakes is a fruitful area of future research.
Successful fisheries management requires accounting for climate change and other global and regional stressors (Paukert et al. 2016). Sustainable harvest policies can differ substantially as ecosystem productivity changes (Walters et al. 2008). Changes in habitat area can alter population growth rates such that harvest levels that were once sustainable are no longer so (Carpenter et al. 2017). Our results suggest that altering harvest in response to changing conditions may allow Mille Lacs to retain its function as a walleye fishery. Our model estimates safe harvest levels as a function of walleye stock size and TOHA averaged over the previous three years, and could be run annually to estimate safe harvest. Continued monitoring of water clarity and temperature is relatively inexpensive and is already a part of standard monitoring of Mille Lacs and many other lakes; these data can be used to adjust harvest in response to environmental changes. Of course, a management regime that adjusts harvest based on environmental changes will require flexible structures and institutions that can adapt to change (Green et al. 2017), as well as a commitment to sustaining walleye populations over the long term. Accepting harvest reductions is socially and politically difficult, and will require coordination, communication, and collaboration among stakeholders, policy makers, and scientists (i.e., adaptive governance; sensu Folke et al. 2005). However, under rapidly changing conditions, the potential for exploitation is constrained and reductions in harvest can facilitate adaptation (Roberts et al. 2017). Our results, like all models fit to empirical data, are bounded by the range of variability in our dataset and the assumptions in our model. Parameter estimates and management implications may change if critical variables move outside the range we have previously observed. Therefore, sustained monitoring is essential for managing walleye and adapting to rapid and uncertain ecosystem change.
Theories of global change suggest the need for approaches based on a SOS for living resources, whereby harvest or other local variables are adjusted to compensate for large-scale changes in climate or other drivers (Scheffer et al. 2015, Carpenter et al. 2017). This study uses field observations to estimate the SOS, explain changes in walleye stocks in relation to the SOS, and suggest local changes in harvest that could sustain a valuable walleye stock in its SOS for the long term. We thereby demonstrate empirically a general approach that could be used to sustain diverse living resources in a time of extensive long-term environmental change (USGCRP 2018).