Ecosystem services in connected catchment to coast ecosystems: Monitoring to detect emerging trends

. We ﬁ rst quanti ﬁ ed covariation among three ecological indicators that had time-series data: pasture biomass, vegetation greenness and barramundi catch per unit effort. Higher values of all indicators oc-curred in years with greater river ﬂ ow. We then predicted the emergence times for each indicator, as the time taken for a trend in an indicator to emerge from the background of natural variation. Emergence times were > 10 years in all cases, quanti ﬁ ed at 80 % and higher con ﬁ dence levels. Past trends and current status of ecosystem service ﬂ ows are often used by decision makers to directly inform near-term actions, particularly for provisioning services (such as barramundi catch) due to their important contribution to regional economies. We found that ecological indicators could be used to assess historical performance over decadal timespans, but not as short-term indicators of recent change. More generally, we offer an approach to performance testing of indicators. This approach could be useful for quantifying timescales of ecosystem response in systems where cross-ecosystem connections are important.


Introduction
The supply of ecosystem services can depend on ecological processes that operate across multiple ecosystems (Barbier et al., 2011;Lowe et al., 2022).Interconnections among ecosystems are facilitated by the movement of plant propagules, migration of animals and the transport of physical matter in air and water.In particular, the flow of water from catchments to coast creates many challenges for the management of human activities that may impact upon the supply of ecosystem services (Alvarez-Romero et al., 2011).For example, the capture of water in dams for agricultural, power or municipal uses upstream interferes with the downstream supply of ecosystem services all the way to the ocean, potentially affecting ecosystem inputs to agriculture on floodplains, the productivity of estuarine and coastal marine fisheries (Gillanders et al., 2011) and the capture of carbon in estuarine habitats (Macreadie et al., 2017).Therefore, monitoring of multiple ecological indicators that together span interconnected systems is important to characterize the processes that support ecosystem services in catchments and inform the management of human activities.This monitoring of ecological processes plays multiple roles in management decisions.First, it informs decision makers about trends in downstream ecosystem processes so they can respond proactively to address changes in hydrology upstream.Second, it helps to detect change in ecosystems and with attribution to management actions.Finally, monitoring can hold decision makers accountable for actions that affect the integrity of ecosystems downstream (Stevenson et al., 2021).
A challenge for monitoring programs that span a whole catchment is to attribute trends in ecological processes to changes in water flow, because the causes of changes in water flow are separated in time and space from their impacts (Alvarez-Romero et al., 2011;Brown et al., 2019).Interannual variation in flow linked to long-term climate cycles (e.g.El Nino southern oscillation) is common in many systems, and this variation can mask long-term trends caused by climate change or human interventions.Monitoring design, such as 'before/after, control/impact' approaches can be used to resolve the attribution problem (Underwood, 1992).However, in whole of catchment monitoring programs, commonly there are no controls and monitoring programs are not designed with a specific intervention, such as water extraction for irrigation, as the focus.In such cases we need to know how long underlying trends in hydrology/flow will take to manifest as trends in the ecological indicators.
The testing of indicator performance, where the ability to detect trends in indicators is explored through simulation models, is a popular method in fisheries science for evaluating indicators (Fulton et al., 2005;Blanchard et al., 2014;Nickols et al., 2019), but it is only just beginning to be applied to ecosystem indicators (Collen and Nicholson, 2014;Watermeyer et al., 2021).Performance testing involves constructing a system model that represents the fundamental system dynamics, and then running simulations that represent the stages of underlying trends, data collection and interpretation of data by decision-makers (Fulton et al., 2005).Simulations are run multiple times so that we can define error bounds on the probability of a trend of a given magnitude being detected with a given level of confidence.This testing can be evaluated through univariate indicators, such as binary indicators of whether a threshold value is reached, or via power analyses that consider the sensitivity and specificity of indicators for different decision thresholds (i.e.risk of missed detections versus false alarms, see Mapstone, 1995).
The development of performance testing methods in fisheries science and biodiversity monitoring has run parallel to the development of the emergence time concept in climate studies (Hawkins and Sutton, 2012).Climate models can be used to quantify when human-induced trends in climate caused by greenhouse gas emissions emerge from background variation in climate.Thus, models of trends underpin the attribution of climate trends to human causes and provide an easily communicable statistic about lags between human actions and detection of changes in the system under study.The emergence time concept could also be applied to the testing of indicators for interconnected catchment ecosystems, where high interannual variation in river flows may mask long-term trends in multiple downstream ecological processes.However, development of such performance testing approaches is often held back by the requirement to develop complex system models.Developing system models for catchment to coast systems is a particularly complex task because the modelling needs to span multiple disciplines (Brown et al., 2019).Catchment to coast systems models are therefore rare.Simpler approaches to performance testing would facilitate more widespread testing of indicators and inform how managers can use indicators.
The suitability of indicators for their intended purpose is a widespread issue that affects conservation monitoring at all scales from global targets in the Convention on Biological Diversity (Nicholson et al., 2021) to local scale decision making (e.g.Jakobsson et al., 2021).Here we develop an approach to quantifying emergence time of multivariate catchment processes.We use a case-study of a riverine-estuarine-coastal marine system in the Gulf of Carpentaria in tropical northern Australia.The Mitchell catchment is earmarked for water resource development, which may provide agricultural benefits, but will also disrupt freshwater flows (CSIRO, 2018) and impact the downstream ecosystems and fisheries that depend on them (Turschwell et al., 2019;Broadley et al., 2020).Similarly, climate change is causing long-term shifts in weather patterns and flows in the region (Karim et al., 2015) We aimed to test commonly used monitoring indicators for their ability to detect change in the interconnected ecosystems of the Mitchell catchment.We sought to fit a process model that quantified the relationship among indicators and that allowed the effect of uncertainty within those indicators to be assessed (Objective 1).From this model we inferred natural bounds of variation for each indicator (Objective 2).We finally quantified the sensitivity of indicators to long-term trends in water flow via their emergence time from background variation (Objective 3).The emergence time for each indicator represents when the impacts of change will become apparent to human observers (Hawkins and Sutton, 2012).

Indicators for the Mitchell catchment
The Mitchell catchment covers an area of over 72,000 km 2 and discharges into the Gulf of Carpentaria (Fig. 1).This catchment is subject to highly seasonal flow regimes.During the dry season many of the ephemeral rivers become disconnected water holes and in the wet season a large area of the floodplain becomes inundated.There is high inter-annual variability in rainfall, which drives inter-annual variation in flow and floodplain productivity (Ndehedehe et al., 2021).
We focused on indicators of capacity to supply ecosystem services that included (i) fish biomass that supports barramundi (Lates calcarifer) catch in a commercial fishery and subsistence harvesting by Indigenous Australians, and (ii) pasture biomass that supports commercial livestock production.In particular, the catadromous barramundi life-cycle exemplifies cross-system dependence (Barramundi migrate from freshwater to estuaries to spawn, Milton and Chenery, 2005).Variability in barramundi biomass may be correlated to variability in pasture growth on flood plains which supports livestock production, because both are dependent on regional cycles in river flow.
Time-series data (Table S1) were collected for the ecosystem service capacity indicators.Stream flow, a hydrological indicator, was included as a driver of ecosystem services.Stream flow data was sourced from the Northern Australia Water Resource Assessment.
Yearly data on barramundi (Lates calcarifer) catch and commercial fishing effort were used to derive an indicator of the barramundi population size.Data were sourced from the Queensland Fisheries 'QFish' data cube (State of Queensland Department of Agriculture Fisheries and Forestry, 2018) for the Mitchell region.Barramundi are an important fishery species in northern Australia, and they are caught in commercial, Indigenous, recreational and charter fisheries.They also play an important cultural role within the ecosystem (Jackson et al., 2012).Barramundi spawn in estuaries, with juveniles then moving upstream to freshwater habitats, returning to saltwater as adults (Robins et al., 2005).Seasonal floodplains provide connections to off-stream freshwater habitats as well as highly productive, short-term habitats.Individuals residing in upstream habitats return to the estuary and saltwater to spawn as adults.Barramundi are impacted by flow regimes through changes in connectivity throughout the catchment which benefit growth and survival (Robins et al., 2005;Balston, 2009;Roberts et al., 2019;Leahy and Robins, 2021).Barramundi opportunistically exploit multiple aquatic habitats at different stages of their life cycle.Therefore, if connectivity between these different ecosystem assets is diminished, this can impact their harvestable biomass in following seasons.Barramundi catch per unit effort (CPUE) is a suitable indicator of the interlinked aquatic ecosystems, because barramundi are a near-apex predator in the system.Their production aggregates multiple trophic levels, and they depend on persistent freshwater in lagoons and waterholes to survive dry periods (Robins et al., 2005).
Normalized difference vegetation index (NDVI) is an indicator of vegetation greenness (Pettorelli et al., 2005).We included NDVI because it may be an alternative indicator for the ecosystem service of supply of forage.NDVI was obtained from the Bureau of Meteorology (Table S1).Pasture biomass was included as an indicator of supply of forage for grazing.Pasture biomass was obtained as the predicted output from Queensland Government's AussieGRASS model (State of Queensland Department of Environment and Science, 2021).AussieGRASS is a pasture production and water balance model that covers the Australian continent (Carter et al., 2000) and is run with daily climate data (Jeffrey et al., 2001).AussieGRASS also uses stream flow as input data and NDVI as one of several calibration datasets; however, since other variables (i.e.rainfall, soil type, temperature and local grazing pressure) are also inputs to AussieGRASS, the relationship between pasture biomass and stream flow is not deterministic.

Objective 1: quantify interrelations among condition indicators
We combined two common modelling approaches to develop a model of indicator interrelationships (Fig. 2).The first was a state-space model of barramundi catch per unit effort (CPUE) dynamics (Millar and Meyer, 2000).State space models are widely applied for modelling fisheries dynamics, informing on reference points for fisheries management and for estimating environmental effects on fishery dynamics (Aeberhard et al., 2018), but have rarely been applied to model whole ecosystem change (Auger-Méthé et al., 2021).The state-space model was used to fit catches (a measure of ecosystem service flow), then estimate latent states for stock biomass (a measure of ecosystem capacity to deliver services) and variation in population growth (Fig. 2).
The second modelling approach was a structural equation model that related the observed indicator variables to a latent (unobserved) state, which we termed the 'catchment condition index' (Fig. 2).Structural equation models can represent correlations among observed indicators and a latent state.The observed variables are then interpreted as outcomes of an underlying causal process, the latent state (Grace et al., 2010).We followed earlier literature in interpreting this latent state as the condition of the interlinked ecosystems of the catchment (Brown and Hamilton, 2018).We linked the two modelling approaches by assuming that the observed indicators of NDVI, pasture biomass, and CPUE were all dependent on an unobserved catchment condition index.Catchment condition, in turn, was assumed to be driven by freshwater flow (Fig. 2).
The state space model had the form: where B t was the latent estimate of barramundi biomass, t were yearly timesteps, r was the intrinsic growth parameter, K was the carrying capacity, C t was catches, q was the catchability coefficient, u t were process errors and σ CPUE was the standard deviation for CPUE errors (we use superscripts to demarcate parameters relating to different variables, subscripts to demarcate indices).B 1 was assumed a fixed fraction of 20 % of the estimate of K (Millar and Meyer, 2000;Streipert et al., 2019).We tested the effects of the assumption that B 1 was 20 % of K and that q was constant by repeating analyses with q increasing by 1 % per annum and with B 1 = 50 and 80 % of K.
The state space model was linked to the structural equation model through a dependence of the process errors on a latent variable that we interpreted as the catchment condition index: where β u was a regression coefficient and ν t was the latent condition index.
The process errors were assumed to be log-normally distributed so they had a multiplicative effect on biomass.The natural log formulation is consistent with population theory because changes in survival have multiplicative effects on population size (Hilborn and Walters, 2013).
The state-space model assumed a one-year lag between flow and changes in barramundi biomass.Flow likely impacts barramundi catches and biomass at a number of different time-scales, as described above, but we assumed a one-year lag for simplicity.To verify that the data did not suggest multi-year lags we confirmed that residuals from the modelled CPUE relationship were not autocorrelated.
NDVI and pasture biomass were related directly to the condition index: where y t refered to either NDVI or pasture biomass and a y is the intercept.We used a normal (not log-normal) distribution for NDVI and pasture biomass because this model assumption best fit the data and we were not modelling population growth for these indicators.NDVI and pasture biomass were mean centred and standardized prior to analysis.We note that it would be theoretically possible to integrate the AussieGRASS model directly into the SEM, however for our purposes this was unnecessary because there were no dynamic feedbacks from pasture production to the other SEM variables.
Finally, catchment condition was estimated as a regression on freshwater flow: where F t was freshwater flow and the standard deviation was arbitrarily set at 1, because ν t is scale invariant.We used the natural log of flow to reduce the leverage that extremely high flow years had on condition.Representing flow in natural log also reflects a hypothesis that the sensitivity of ecological condition to flow decreases at high flows.
The NDVI data started in 1992, so the condition index was estimated only from flow, CPUE and pasture biomass in 1990 and 1991.The model also interpolated the missing NDVI values.We presented results as marginal predictions for each data type and the latent states for barramundi biomass and catchment condition.
We obtained parameter and prediction probabilities within a Bayesian estimation framework.The model was estimated with the STAN programming language, implemented from the R program with the rstan package (Stan Development Team, 2020).We set priors based on prior information for each parameter and prior expectations for the variance in NDVI and CPUE.It was also important to set weakly-informative priors to aid Fig. 2. Overview of the approach showing the measured variables, the latent variables and the scope of the model.algorithm convergence and identification of parameters.In particular, the latent state ν t and the β coefficients can 'sign-switch', e.g.−ν t and +β are exchangeable with +ν t and -β (Hui et al., 2015).Sign-switching makes interpretation of the posteriors challenging, because the interpretation of one parameter depends on the sign of another parameter.Priors and algorithm settings were chosen consciously, rather than using software defaults (Table S2, supplementary material sections 1.1 and 1.2).
The model was verified for its fits to the data and through simulation testing.The model assumptions were checked by confirming independent normally distributed residuals, where residuals were defined as the observed values of CPUE, NDVI and pasture biomass minus their posterior predicted median value.Residual independence was confirmed by plotting the auto-correlation function with the R function 'acf' (base R, R Development Core Team, 2019).Further, we checked the conditional independence assumptions implied by the SEM by calculating the crosscorrelation of residuals across indicators to confirm the model had captured all dependencies between indicator values (Grace et al., 2012).Finally, since we had some missing data, we verified the robustness of parameter estimates to missing data with a simulation study.We created 50 randomised flow time-series; simulated values of NDVI, pasture biomass and CPUE from these timeseries; and then refit the model.We repeated these steps where the first 0, 5, 11 and 15 years of the NDVI data were missing in the model fitting.We could then compare the estimated parameter values to the true values from the simulation (supplementary material section 1.4).
We determined the strength of the relationship between each indicator variable, catchment condition, and flow as the direct effect of condition on the pasture biomass and NDVI indicators (simply the marginal posterior distributions for β y ) and the total effect of flow (in natural log) on the condition index =β ν β y .For barramundi we calculated the direct and total effect sizes on surplus production, where surplus production was defined as the sum of change in biomass and catch (i.e.SP t = B t − B t−1 + C t ).We then plotted the total effect of a 1 S.D. increase in the log of flow (=5.6 times increase in flow) on each indicator and the direct effect of a one S.D. change in catchment condition on each indicator (formulas for effects were derived in supplementary material section 1.5).We also calculated one-tailed probabilities that the total effect was greater than zero.We note that given the prior choices (Table S2) the prior probability that β ν .β y >0 was 0.5.

Objective 2: bounds of natural variation
We used the fitted model to define the bounds of natural variation for each indicator.The bounds were defined to represent the typical ranges of indicator values given historical variation in flow.The historical condition of flow was characterized by fitting an autoregressive moving average model to the historical flow series , then simulating 1000 random flow series from the fitted model.The optimal model for simulating flow was chosen with the AIC and had an autoregressive order = 6 and moving average order = 3. Flow models were fit with the 'Arima' function and simulated with the 'simulate' function from the forecast' R package (Hyndman and Khandakar, 2008).To define bounds for the indicators we then simulated their values from the fitted SEM, where the SEM took each of the 1000 flow series as an input.For NDVI and pasture biomass, the bounds were defined as the marginal probability quantiles integrating over measurement noise (also called predictive intervals) and the flow simulations.The bounds for barramundi CPUE also depended on biomass, so for simulations of barramundi we assumed a biomass fixed at the biomass that gives maximum sustainable yield (B 1 = 50 % of B 0 ).This is a common reference point used in fishery management (Mace, 1994).For the simulations of bounds, we also assumed that biomass was initially at the biomass that gives maximum sustainable yield, doing so avoided transient dynamics affecting the simulations.

Objective 3: emergence time for indicators
We then estimated each indicator's emergence time, given an assumed long-term trend in flow.Environmental change is causing long-term changes in flows in the region (Karim et al., 2015), which may affect interconnected catchment ecosystem services, including barramundi catches and pasture growth.Thus, our scenarios represent the ability to detect this long-term change under an assumed trend magnitude.We simulated flow series and added a linear declining trend in flow, as an example of a situation where environmental change or water resource development could increasingly constrain ecosystem services.We simulated 1000 iterations of 30-year flow time-series for each trend in flow and then predicted values of all indicators.We then defined emergence times as the time taken for the mean indicator value to cross the lower 5 %, 20 % and 40 % probability quantiles.The different quantiles represent different degrees of confidence for emergence, the lower the quantile the greater the confidence the trend has emerged from background variation.This quantilesbased approach is analogous to a power analysis in frequentist statistics, where a trend is refitted to simulated data to test a model's ability to detect change.In our case we used a Bayesian analysis, which allowed us to propagate uncertainty from parameter estimates through to the posterior, so refitting was not required.
Code and data to repeat the analyses is available at: https://github.com/cbrown5/ecological-condition-latent-model.

Objective 1: interrelationships among indicators
Verification of the model fits indicated that there was no autocorrelation in the residuals for the fitted values for any of the indicators (Fig. S1).Overall, the R 2 was highest for log CPUE (0.91), moderate for NDVI (0.43) and lowest for pasture biomass (0.09) (Fig. S2).Sensitivity analyses in which we varied the assumptions regarding constant catchability, the initial biomass ratio and higher fishing effort indicated the main results were robust to alternative assumptions (Figs.S2, S3 and S4).Consequently, below we present only the main results assuming stable catchability and an initial biomass ratio of 20 %.
Model fits indicated an overall increase in CPUE from 1990 to 2012, and then a slight decrease from 2012 to 2015.NDVI and pasture biomass had no long-term trends, but varied year-to-year (Fig. 3).NDVI and barramundi surplus production were strongly related to river flow (Fig. 4A, pr(β ν β y .)> 0 = 0.99 and 0.96 respectively), whereas pasture biomass had a somewhat weaker relationship with river flow (Fig. 4A, pr(β ν β y .)> 0 = 0.93).Therefore, the catchment condition variable primarily captured the correlation between flow, NDVI and barramundi biomass change (Fig. 3E).Catchment condition showed peaks in 1991, 1998-99 and around 2010 that were driven by increased flow in those years.But overall catchment condition did not vary strongly from the long-term average (95 % CIs generally overlap zero).Barramundi biomass was predicted to have increased from 0.2 to 0.5 as a fraction of B 0 over 1990-2012, and to have declined to ~0.3 by 2015.The decline in 2012-2015 corresponded to a decline in catchment condition and a decline in flow.
A simulation study where we varied the number of years of missing NDVI data suggested the model's predictions for the variables and catchment condition were not biased by missing data, but the predictions for NDVI were less accurate with a greater amount of missing data (Fig. S5).Coverage by the 95 % C.I.s was always >0.95 (Table S3).The parameter estimate for the effect of catchment condition on NDVI was most affected by missing data, and the parameter estimates for pasture biomass and barramundi biomass were least affected.

Objective 2: bounds of natural variation
The bounds for each indicator had broad intervals (Fig. 5).For CPUE the 90 % probability interval ranged from 17 to 70 t per fishing year, compared to a historical range of 19-50 t per fishing year.For NDVI the 90 % probability interval ranged from 0.19 to 0.26 units, compared to a historical range of 0.20-0.30units.For pasture biomass the 90 % probability interval ranged from 1970 to 2998 units, compared to a historical range of 2000-4134 units.This broad variation in the probability bounds was caused by the high variation in year-to-year flow, to which CPUE was particularly sensitive.

Objective 3: emergence times
Emergence times were greatest for pasture biomass and lowest for CPUE (Figs. 5 and 6).For instance, if flow decreased by multiples of 0.9 (10 % per year), the time required for mean CPUE to fall below the 5 % quantile (90 % probability intervals) of historical variation was 16 years, but >30 years for pasture biomass and NDVI.In comparison, emergence time for the 5 % quantile for a 10 % per year trend in stream flow was 20 years (Fig. S6).In general flow had slower emergence times than CPUE, but faster emergence times than pasture biomass and NDVI.Emergence times were lower for higher quantiles, reflecting the reduced level of confidence that the trend had emerged from background variation.The differences in emergence times among indicators became smaller as the quantiles were

Relationships among indicators
We identified relationships among all indicators and flow, and particularly for barramundi catch per unit effort and NDVI.There is high interannual variation in the Mitchell catchment, driven by interannual variation in rainfall.Wet years drive increased vegetation growth and barramundi production, which was reflected in correlations among river flow, NDVI and barramundi catch.Thus, catchment condition, defined here as covariation in the indicators, largely reflected year-to-year variation in flow.This variation underscores the importance of quantifying interconnections between assets via river flow and among ecosystem services.The implication is that change in one variable (flow) can cause changes in multiple processes, and consequently there would be correlated changes in the annual capacity of the system to deliver different ecosystem services.
By defining the catchment condition index as a latent variable we were able to partition out environmental effects shared by all indicators (captured in the catchment condition index) from those that were unique to specific indicators (variation in indicators that is not explained by the condition index).The ability to identify shared trends may be desirable in an indicator, because it facilitates the attribution of changes that affect only specific indicators versus those that impact the overall condition of the ecosystem (Sutherland et al., 2016;Kupschus et al., 2016).Over time there was an increasing trend in barramundi CPUE from the 1990s to 2015.This was driven by high flows 2005-2015 and management changes in the fishery in the late 1990s and early 2000s, which reduced fishing effort and put increased restrictions on fishing gear (Streipert et al., 2019).NDVI had a strong correlation with flow and CPUE, but in contrast to CPUE, NDVI did not show as strong a long-term trend, because it was not affected by fishery management changes.Pasture biomass had the weakest correlation with the other indicators, perhaps because it was more strongly driven by rainfall that is local to a region and less by inundation caused by increased river flow.
The usual approach to defining joint indices of multiple ecosystem components has been spatial and thematic aggregation by summing spatiallyspecific counts or by averaging and then normalising spatially-specific indicators (Maes et al., 2020).The aggregation procedure requires the use of some form(s) of aggregation functions and weights.Our model-based approach to defining a condition index can also be interpreted as a weighted sum of the indicators, where the weights were estimated empirically.Each approach has its own advantages and disadvantages.Our approach requires time-series data, which are not available for many indicators or in some regions (Dvarskas, 2018).Expert weights do not require historical data and are simple to compute but can also be arbitrary.Arbitrary weights are not repeatable across different groups of experts and may mis-represent management priorities (Game et al., 2013).Further, our process model separates the effects of flow on CPUE (modelled as changes in population growth) from the effects of fishery management on CPUE (modelled as changes in fishing effort).Simply summing weighted CPUE with other indicators would tend to overinflate the influence of barramundi fishery management on catchment condition.Thus, a key advantage of our approach is that it can partition changes unique to individual indicators, versus changes that are shared across the ecosystem.Identifying redundancy and complementarities among different indicators is an important step in ensuring that monitoring programs can focus on a parsimonious set of indicators that is comprehensive of the ecosystem's dynamics (Czúcz et al., 2021).

Comparison of our approach to other methods
Methods for evaluating indicator performance have been developed for fisheries and biodiversity indicators (Fulton et al., 2005;Rowland et al., 2020;Watermeyer et al., 2021).Earlier approaches tended to use mechanistic ecosystem models, which are calibrated to be representative of ecosystem dynamics.These ecosystem models represent multiple variables, so performance of indicators representing different parts of the ecosystem can be assessed.Indicator evaluation has also been conducted through pattern-based approaches to modelling trends, where time-series with realistic properties are simulated and models are refit to determine the ability of the model to detect trends in noisy time-series data (e.g.Knape, 2016).Our approach is a blend of statistical (the structural equation model) and process-based models (the state-space model).We fitted the model directly to data, which enabled us to quantify stochasticity resulting from different sourcesin this case inter-annual variation in flow as opposed to observational errors in indicators.Thus, our approach combines strengths of both statistical trend detection and process-based modelling.This type of model warrants further development for assessing ecosystem level change, particularly where there is a need to differentiate the different sources of stochasticity in observed time-series (Auger-Méthé et al., 2021).
Several modelling toolboxes have been previously developed to quantify connections among the condition of ecosystems and ecosystem service flows, but none of these directly address the challenge of effectively detecting trends in indicators.Examples include the soil and water assessment tool (Bieger et al., 2017), the Integrated valuation of ecosystem services tradeoff (INVEST) tool (Guerry et al., 2012) and globally applicable artificial intelligence tools (Martínez-López et al., 2019).These tools have been applied for measuring connected ecosystem services across catchments, in particular regulating ecosystem services relating to water quantity, water quality and sediment erosion (Bieger et al., 2017;Hamel et al., 2017).All share a common approach, based on parameterising models of soil erosion and run-off and then analysing predicted service flows spatially in a geographic information system.These toolboxes are designed for making service flow projections, given sufficient information exists to parameterize their models.Mechanistic ecosystem service models could be used in a similar way to our approach of evaluating indicator performance.This would additionally require incorporating stochasticity into the ecosystem service model, such as variation in flow, so that detection of trends could be studied relative to background environmental variation.

Caveats and limitations
There are some important caveats to consider in our analysis.First, there was missing data in some of the indicator time-series.Missing data in time-series is a common problem in modelling and will reduce the power to detect long-term change (Dvarskas, 2018).The latent variable approach enabled interpolating missing data, on the strength of correlated variables (Hui et al., 2015).A simulation study suggested accuracy of predictions for a variable was not sensitive to missing data in other variables.Thus, we suggest our approach could have value more generally for creating indices of catchment condition when there are gaps in time-series.
The analysis of barramundi used a simplified stock assessment model.For many fisheries additional data may be available, such as data on age structure, and this could improve the accuracy of the modelling.Our simple stock assessment did compare well with a more quantitatively complex stock assessment, with the stock status estimates being similar (Streipert et al., 2019) this is the consequence of both assessment models being driven primarily by trends in CPUE and the assumption of an depletion of 0.2 B 0 in 1990.More importantly, the availability of barramundi to estuarine fisheries, as well as productivity, is also likely to be affected by flow (Robins et al., 2005).An effect of flow on availability to the fishery would mean flow affects the catchability parameter instead of, or as well as, population growth.Future modelling could potentially separate the different types of flow effects, but this would require additional data types, because productivity and catchability parameters are close to being interchangeable in unstructured surplus production models.
An important caveat is that our indicators represent only one aspect each of the processes supporting capacity to supply ecosystem services in the Mitchell catchment.We focussed on assets that support commercial fish harvest and commercial livestock grazing, because there were sufficient time-series data for these indicators to assess emergence times.As an example of what the indicators miss, Indigenous uses of ecosystems are also an important component of the ecosystem services in Northern Australian catchments (Scheepers and Jackson, 2012;Jackson et al., 2012Jackson et al., , 2014)), which may be particularly sensitive to development (Stoeckl et al., 2013).Our modelling did not include condition directly relating to Indigenous services, due to lack of time-series data.It is likely that our conclusions regarding emergence time would also apply to those assets supporting traditional uses, such as water holes, given the interconnections in this catchment.Future modelling studies could integrate Indigenous knowledge to enhance the predictive power of ecosystem service models.For example, Indigenous knowledge on the spatial distribution of barramundi (Scheepers and Jackson, 2012) could be incorporated into ecosystem service models (e.g.Guerry et al., 2012).This would enhance our ability to predict change in ecosystem services, as well as quantify their value in a broader range cultural contexts.

Implications for management
Our results suggest what are appropriate and inappropriate uses of these indicators by management.CPUE had the shortest emergence times, so will be the most sensitive indicator of long-term trends in flow.CPUE had shorter emergence times than even flow itself.The multi-year life-history of barramundi smooths interannual variation in flow, so trends were slightly easier to detect.However, the high interannual variation in flow and catchment condition will create challenges for detecting change in the Mitchell catchment.Our estimates of the interannual variation for each indicator were broad and thus our model predicted that emergence times across all indicators at typical confidence bounds (10 % or 5 %) were >10 years and sometimes >30 years for all but the most extreme trends in flow.
These long emergence times mean the indicators will not be useful for managers aiming to use the model to pre-emptively respond to change, because the changes will not be detected on a management relevant timespan (Stevenson et al., 2021).Rather, the indicators will capture long-term lowbandwidth changes in the ecosystem.For comparison, proposed water developments in the Mitchell catchment are expected to cause much smaller reductions in flow, in order of 1-5 % per year (Petheram et al., 2018), than most of the trends we explored here.Such water developments however are predicted to impact estuarine ecosystem services (Broadley et al., 2020;Leahy and Robins, 2021).Thus, our findings suggest management planning should not rely on indicators of catchment condition.Instead, pressure indicators may better pre-empt change in catchment condition (Czúcz et al., 2021).Predictive modelling of the response to ecosystem services to future changes in land-use and water flow is also needed to inform pre-emptive management in the Mitchell catchment (McMahon et al., 2022).

Conclusion
We analyzed change in indicators for interconnected ecosystem assets and found a high level of variation.This high level of variation means long time-series are required to detect change in the capacity of the ecosystem to supply ecosystem services.Our quantification of timescales to detect change provides important context for policy-makers when they are considering indicators.For instance, we urge caution when using indicators like the ones analyzed here in cost-benefit analysis when evaluating water resource development proposals.Economic variables may respond rapidly and positively to development, but the consequent degradation of catchment condition and ecosystem's capacity to supply ecosystem services may take much longer to be noticed.These temporal differences in response times may not be adequately captured in cost benefit analysis conducted prior to approval, which may lead to decisions that are biased towards development opportunities that bring short-term benefits but have long-term environmental costs that are discounted to arrive at the present value.Our approach could be used to support development of indicator sets for attribution of trends in human impacts, catchment condition accounts and monitoring of human impacts.

Fig. 3 .
Fig. 3. Model fits to data for CPUE (A), NDVI (B); pasture biomass (C); flow (D) and the estimated latent variables for catchment condition (variable ν t in the model) (E) and barramundi biomass relative to its carrying capacity as estimated by the model (F).For model predictions, lines show median and shaded areas show 95 % CIs.Darker areas in E-F indicate 50 % CIs.The solid horizontal line in (E) indicates the mean value, in (F) it indicates 50 % carrying capacity.

Fig. 4 .
Fig. 4. Strength of the relationships between each indicator and (A) a 1 S.D. increase in log of flow; and (B) a 1 S.D. increase in the catchment condition index.Note that NDVI and pasture are differences on a scale standardized to 1 S.D., whereas surplus production is on a multiplicative scale.Horizontal bars show median, thick solid lines show 75 % probability quantiles and thin grey lines show 95 % probability quantiles.

Fig. 5 .
Fig. 5. Bounds for catch per unit effort (A) NDVI (B) and pasture biomass (C).Grey zones show 40 % (dark grey), 60 % (medium grey) and 90 % (light grey) probability intervals.Examples of trend lines shown in blue shades, as proportion of the mean.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)