Simple Bayesian reconstruction and forecasting of stream water temperature for ecologists—A tool using air temperature, optionally flow, in a time series decomposition approach

Mitigating the impacts of global warming on wildlife entails four practical steps. First, we need to study how processes of interest vary with temperature. Second, we need to build good temperature scenarios. Third, processes can be forecast accordingly. Only then can we perform the fourth step, testing mitigating measures. While having good temperature data is essential, this is not straightforward for stream ecologists and managers. Water temperature (WT) data are often short and incomplete and future projections are currently not routinely available. There is a need for generic models which address this data gap with good resolution and current models are partly lacking. Here, we expand a previously published hierarchical Bayesian model that was driven by air temperature (AT) and flow (Q) as a second covariate. The new model can hindcast and forecast WT time series at a daily time step. It also allows a better appraisal of real uncertainties in the warming of water temperatures in rivers compared to the previous version, stemming from its hybrid structure between time series decomposition and regression. This model decomposes all-time series using seasonal sinusoidal periodic signals and time varying means and amplitudes. It then links the contrasted frequency signals of WT (daily and six month) through regressions to that of AT and optionally Q for better resolution. We apply this model to two contrasting case study rivers. For one case study, AT only is available as a covariate. This expanded model further improves the already good fitting and predictive capabilities of its earlier version while additionally highlighting warming uncertainties. The code is available online and can easily be run for other temperate rivers.

Understanding and hence mitigating climate change impacts on biological communities involves four steps.First, we need to quantify the interactions between the process of interest and temperature.Second, we need to predict temperature under future greenhouse gas (GHG) scenarios.Third, we need to forecast our process under these future climate scenarios.Then, as the fourth step, we can work towards mitigating impacts, preferentially through an adaptive management framework [15,16].Robust ecological and temperature data are essential cornerstones of this process.
Managing freshwater species conservation in the face of climate change is particularly challenging.First, such species are highly susceptible to climate change as they are mostly ectothermic and thus sensitive to temperature [17][18][19].Shifts in spatial distribution are also constrained by the river network and habitat fragmentation [20].Second, water temperature time series are limited compared to those of air temperature.Indeed, water temperature time series are often quite short and prone to missing data, and unlike air temperature, are not generally included as outputs of climate change models.Generic models to complete and hindcast/forecast water temperature based on covariates such air temperature and flow are needed by ecologists and managers.
Different types of stream water temperature models exist but they present caveats to the ecologist [21].Process based models, e.g.energy budget models [22][23][24] are realistic but complex.They require a lot of environmental and geological data, and are also location specific.This makes them of limited use for long term and large-scale forecasting.Statistical models on the other hand are less data demanding and rather simple [25].Examples of this model type include regular time series models (AR, ARIMA), periodic autoregressive models, K nearest neighbors methods and neural network approaches [26,27].These methods perform well for short term forecasting and filling short gaps.On the contrary, they perform poorly for long term predictions [25].The last category usually encountered is regressions models [28][29][30].While they may seem adequate, they can lead to significant biases in more long-term forecasting.Indeed, they do not disentangle long term trends, seasonality and short-term variations.An extensive review of stream water temperature models can be found in [31] and allows for a better appraisal of the limitations exposed.
As a workaround, [21] proposed a Bayesian hierarchical approach which was a hybrid between time series decomposition and regression.This approach required water temperature together with air temperature and flow time series as covariates.It separated out long term, seasonal and short terms components of the time series and linked these components through correlation.This method has several advantages over regression, in that it outperformed correlation in both fit and forecasting capabilities, and appeared to be unbiased for long term forecasting.This hybrid approach also highlighted the true, and high uncertainty, in future river temperature warming.However, this approach lacked fine-scale modelling of deviations (i.e.daily) around long term and seasonal components, and required additional refinement.Specifically, short-term deviations were modeled using a first order autoregressive (AR1) process which limits fine hindcasting/forecasting abilities, and the model worked on 5-day average temperature which is rather coarse.In this paper, we improve on the previous model by i) working on daily temperature and ii) modeling daily deviations as depending on that of air temperature (AT) and flow (Q).Lastly, we offer the possibility of running this model without water flow data, as these are often not collected routinely.These new features, coupled with openly accessible code should be helpful to ecologists and stakeholders.
In this paper, we describe the fit and forecasting performances of this extended model and compare it to that of [21] on two rivers with contrasting bio-geographical conditions and sizes.The Rough River is a small Irish spate river while the Scorff river is of medium size and located in France.For the Rough river, air temperature is the only covariate.As in [21], we use a fully hierarchical Bayesian framework.This provides a probabilistic rationale to quantify uncertainty in inferences and forecasting which further ease out ecological analysis and risk management [32].

Material and method
The new model version presented here is developed with reference to that of [21].For this reason, we refer to the model in the previous publication as M 1 (for model 1).Our newest version described thereafter is called M 2 .

Model M 2 structure
Model M 2 has a fully Bayesian hierarchical structures and produces estimates of water temperature (WT) based on air temperature (AT) and discharge (Q).This is achieved using three fully integrated modules.Module 1 desegregates WT, AT and Q time series into long term variations, seasonal fluctuations and short-term variations.Module 2 links WT components to that of the AT and Q (also referred to as covariates).Module 3 then provides estimates of all unknown WT based on modules 1 and 2. Estimated WT can range anywhere from a few missing values up to several years in case of hindcasting or forecasting.The difference between model M 1 [21] and M 2 lies in module 2 and also in the way daily variation in WT are modeled.In case of model M 1 , daily variation was assumed to follow an autoregressive process.Here we expand module 2 to link these to those of the covariates.
Module 1 of the model M 2 desegregates WT and covariates time series (X t in the following equation) as follow: α y and β y are the mean and amplitude of a time window.n is the number of time steps per year.In this article n is equal to the number of days within a year (365 or 366) while we set the time window y to 6 months as it offers a good tradeoff (see [21]).t 0 sets the seasonal signal on the year.For AT and Q, we modeled � t using first order autoregressive (AR1) processes centered on 0 (with autocorrelation coefficient ρ and standard deviation σ).[21] modeled WT daily deviations the exact same way for model M 1 .This is where this new model differ and we detail this as part of the second module below.Module 2 links WT time series to that AT and Q (if available).It comprises two sub-modules (2a and 2b).Sub-module 2a corresponds to the entirety of module 2 from [21].In this first sub-module, we link parameter α y and β y of the WT time series to that of the covariates time series.To do so, we first reparametrized the sine signal within the time series using max y = α y + β y and min y = α y − β y .We then model WT maxima max WT y and minima min WT y as depending upon those of the covariates.Corresponding equations are: We kept this parametrization from [21] because of its conceptual sense.Let's assume that θ 1 is positive, θ 2 is negative and both θ' 1 and θ' 2 are positive.This corresponds to assuming that i) warm AT leads to warm WT, ii) high Q means warms WT in winter and a cooling effect in summer.
In addition, module 2b now links the short-term variations of WT to that of the covariates.We parametrized this addition to the original model M 1 as follows: WT short term deviations thus are directly proportional to those of AT, although usually buffered.For Q, we included a sinusoidal component.This allows the impact of discharge deviations to change along the year.It shifts from positive in winter to negative in summer with a smooth transition along the year.We assume it as strictly antiphasing the WT seasonal signal, hence the inclusion of n\2.Because we use the log discharge, small changes in summer flow have an impact more comparable to that of bigger changes in winter.We modeled the remaining errors (z WT t ) using a white noise process.Lastly, module 3 corresponds to the estimate of unknown temperatures.As our model is fully hierarchical, this module depends upon and integrates with modules 1 and 2. Consequently, the model propagates all uncertainties when estimating unknown WT based on covariates.In Bayesian software such BUGS, JAGS, STAN or NIMBLE ( [33] and references therein), this is performed simultaneously to the fit.AT and Q time series just have to include hindcasting and/or forecasting parts if any.Module 1 also decomposes these 'extra' data while module 2 produces the WT estimates.This process is often referred to as deriving posterior predictive distributions.

Model application
To assess the performances of the model (M 2 ), we compared it to model (M 1 ) on two study sites with contrasting sizes and environmental conditions and tested both fit and forecasting capacities.
Assessing models' performances.We checked the consistency of both models M 1 and M 2 a posteriori with the data using the χ 2 discrepancy statistics (Gelman et al. 2015).In particular, this posterior check allows to assess whether replicated data are similar to the original data.Its formula is: Var WT y;t jy where E(WT y,t |θ) and Var(WT y,t |θ) are the expected mean and variance of WT conditionally upon all parameters (θ).For each set of parameters drawn from the joint posterior distribution, we computed replicated and observed WT χ 2 values based on the equation.Then we computed a Bayesian p-value giving the probability that χ 2 (WT rep |θ) > χ 2 (WT obs |θ).Bayesian p-values close to 0.5 suggest congruency between the model and the data.Very high or very low (about 0.95 or 0.05) values provide serious inconsistency warnings.
Secondly, we also compared the fitting performances of models M 1 and M 2 using the Deviance Information Criterion (DIC, [34]).DIC allows the comparison of goodness of fit while penalizing complexity in a way very similar to that of the Akaike Information Criterion.The smaller the DIC the better.A five points reduction is usually considered as indicative of a significant improvement.
Lastly, we checked the predictive performances of the models.For this we used cross a validation approach.Specifically, we used 2/3 of historical data as training set and forecast the last third.We used the Root Mean Square Error (RMSE) of the estimates to compare performances.
Bayesian computations details.We implemented models M 1 and M 2 in JAGS [35] through R [36,37].This software approximates the parameters posterior distributions using Monte Carlo Markov Chains (MCMC) algorithms.Specifically, we ran three MCMC in parallel for each model fit.We kept 10000 draws for each one after both a thinning of 25 and an initial burning period of 10000 draws.All diagnostics suggested the MCMC converged successfully.Priors used were all uninformative relative to the data (see S1 Table in S1 File).
The code developed together with one of the data-sets are available on GitHub (https:// github.com/GuillaumeBal/2023.bayes.stream.temperature).Readers can change the input data file with their own and run the model.Doing so does not require any in depth knowledge of both R and JAGS.
Study sites.WT, AT and Q time series come from the Rough (Mayo, Ireland) and Scorff (Brittany, France) rivers (Fig 1).Both rivers support long-term environmental monitoring and fisheries related research [5,43,44].They contrast in size, location, as well as available data (Table 1).In particular, reliable Q time series are lacking on the Rough river.This allowed for testing the performances of our approach with only one covariate.
The Scorff River flows into the Atlantic Ocean (Fig 1).INRA staff members measure daily WT at the Moulin des Princes station with Tidbit data loggers (precision of 0.2˚C, [45]).Daily AT come from the Lorient (LannBihouee ´airport meteorological station operated by the French weather services (https://meteofrance.com/), 9 km away from the water temperature measurement station.Q records come from 8 km upstream of the Moulin des Princes station where a flow measuring station is operated by the French hydrological services (https://hydro.eaufrance.fr).The three time series are 13 years long (1995 to 2007) and 9.6% data are missing on average.This river is part of a larger research observatory where diadromous fish population are extensively monitored and studied (ORE DiaPFC).
The Rough River (also called the Srahrevagh River) is a tributary of the Srahmore River in the Burrishoole Catchment.The Srahmore flows into Lough Feeagh which is connected to the saline coastal lagoon Lough Furnace.The entire system discharges into Clew Bay on the Atlantic Ocean (1.Data cover the period 2002-2016.Daily WT was measured using a StowAway TidbiT temperature data logger from Onset (TB132-05+37).Air temperature data was recorded at the Newport (Furnace) manual weather station (jointly run by Met E ´ireann and the Marine Institute, station 833).Maximum and minimum temperatures are recorded using mercury thermometers, and average daily temperature is taken as the average of these two readings.Data were downloaded from www.met.ie,licensed under a Creative Commons Attribution 4.0 International (CC-By 4.0) License.
Some more details about the catchments are provided in Table 1.Raw time series are shown within the supporting information file (S1 Fig in S1 File).
Forecasting under climate change.We used full joint posteriors from model M 1 and M 2 to project WT temperature from AT warming projections.In particular, we picked a range of AT warming of 0 to 5˚C according to the latest IPCC estimates [46], excluding the most extreme warming scenario.We did this for several reasons.First, [21] revealed WT warming should be lower than that of AT and quite uncertain.We wanted to check whether this result was robust to the updated time step and model structure.We also wanted to see whether two rather contrasted rivers would differ in their response.

Results
The results section is composed of four parts.First, we show differences in fitted series of α and β with both models.Then, we focus on quantitative results showing how the updated approach of modelling the water temperature residuals as depending upon those of air temperature and flow outperformed the version with a simple AR1 process, both for forecasting and fitting.Third, we look at parameters of model M 2 and highlight how those shared with the first model evolved to better understand model differences.Lastly, we showcase forecasting and missing value estimates from model M 2 .

Seasonal variation in the time series
Comparing time series of fitted α and β between models M 1 and M 2 provides basic insights into differences in their fit and behavior.The ranges of values are greater for model M 2 than for model M 1 whatever the parameter series and river considered (Fig 2).This results in statistically significant changes in the variance of median values.For instance, the variance of the median estimates of the α WT posterior series on the Scorff river was 0.47 for model M 2 versus only 0.17 with model M 2 (p.value = 0.01, Fischer test).While series ranges differed, mean values of median estimates for all series were equivalent.The same was observed for average 95 and 50% credible intervals.
In case of AT and Q, α and β time series derived using both models were quite similar (S2, S3 Figs within the S1 File).Mean values, variances and 95 and 50% credible intervals of the time series deviated from each other only fractionally.Only a few posteriors, such as the last couple for the Rough river α AT , were obviously different.

Fit and predictive performances comparison
Overall, results suggest better predictive performances for model M 2 on both rivers.In addition, fit statistics that indicate good predictive performances were more apparent for model M 2 .χ 2 discrepancies and associated p-values were better on both rivers with model M 2 but pvalues obtained with model M 1 did not depart enough from 0.5 to be indicative of an inconsistency between this model and the data (Table 2).
The only exception to the general superior performances of model M 2 were DIC values (Table 2) for the Scorff River.Indeed, differences of several hundred of points in deviance in   favor of model M 1 were highlighted by the analysis on that river.The same was consequently observed for the DIC value.This better fit performance of model M 1 on the Scorff river resulted from its highest fit flexibility when a lot of data are missing that is linked to modeling the residuals with an AR1 process.We also noted that the lack of covariate data tends to increase pD values in case of model M 2 as observed for the Rough river.Discrepancies between 6-month averages of replicated data and observed water temperature were clearly in favor of the more complex model M 2 (Fig 3).Indeed, related posterior distribution for M 2 appeared both better centered on zero but also exhibited a much narrower credible interval.For the Scorff river, the average 95% interval dropped from 1.63 to 0.48˚C.For the Rough River, the corresponding values decreased from 0.99 to 0.29˚C.We can also see that the distributions for the second half of years 2009 and 2010 on the Rough River are far RMSE values clearly highlighted the better forecasting performances of model M 2 .Indeed, incorporating the link between residuals when modeling water temperature allowed for a reduction in RMSE of 41 and 35% for the Rough and Scorff River respectively (Table 2).The forecasting capabilities of model M 2 are further highlighted on Fig 4 .Here, we selected a small portion of historical water temperature data from each river, containing both low and high daily variations.On both rivers, observed and predicted daily variations follow each other closely.High frequency changes are well captured by model M 2 , as can be seen for early summer (days 75-100) for both rivers (Fig 4).The average 50% credible interval around predictions is about 0.75˚C for both rivers, while the 95% credible interval spans slightly more than 4.5˚C (Fig 4).

Posteriors of parameters linking water temperature to the covariates
Modeling the link between the residuals of the different time series resulted in significant differences in θ posterior distributions (Fig 5).Indeed, three main types of changes were noticeable.First, posteriors of θ parameters tightened.This was particularly significant on the Scorff River with, for instance, reduction of about 50% in the 95% credible of parameters y 1 0 , θ 2 and y 2 0 .Second, posteriors of θ 1 and y 1 0 shifted toward stronger positive values.With model M 2 , the credible intervals of parameters linked to LFL were well centered on 0 compared to that of model M 1 .Meanwhile, the uncertainty around the regression appears greater with model M 2 than with model M 1 (see σ posteriors, Fig 5).Parameters δ and γ linking daily variations in water temperature to that of the covariates were estimated as positive together with tight posterior distributions on both rivers (Fig 6).All posteriors were centered on values between 0.5 and 0.7 suggesting changes in water temperature to be buffered when compared to those of covariates.The posterior of parameter γ appeared about 10 times wider than that of δ for the Scorff river (0.28 vs 0.03 95% credible interval), suggesting less precision into the impact of flow.

Forecast warming
Model M 2 predicted higher warming of WT than model M 1 on both rivers (Fig 7).These differences are directly attributable to the differences in θ posteriors detailed above (Fig 5).Consequently, the difference between models was greater on the Scorff river.The WT median warming predicted by model M 2 for the Scorff closely followed that of AT to peak at 4.59˚C for an AT warming of 5˚C (median value), indicating that AT warming is still likely to be a bit faster than WT warming there.Meanwhile, model M 1 median warming prediction was only 1.33 and its 95% credible interval barely overlapped that of model M 2 .For the Rough river, the WT warming predictions overlapped more between models but model M 2 suggested it could likely be faster there too.Lastly, WT warming predictions coming from model M 2 were more uncertain than with model M 1 .The average 95% credible of predicted WT warming were of 1.88 and 2.01˚C on the Scorff river with model M 1 and model M 2 respectively for a 5˚C AT warming.On the Rough river, corresponding values were of 1.23 and 1.61˚C.

Discussion
Expanding the time series decomposition-based approach of [21] to better model daily variation in WT in addition to seasonal variations lead to greatly improved WT predictive abilities.In doing so, RMSE decreased by about 41% on the Rough river and about 35% on the Scorff River.Discrepancies between 6-month averages of replicated data and observed water temperature were also lower with the new M 2 model than with the previous M 1 model.This is indicative of a better fit replicated data generated under model M 2 are more similar to observed data.It also means forecast or hindcast 6-month means should be more reliable.The use of an AR1 structure to model WT residuals as in model M 1 is very flexible and allows for capturing variability in fit.However, this can be considered as an over-fit when considering our goals.Indeed, model M 1 was less successful than M 2 at producing long term hindcasts or predictions.
Interestingly, the new model suggested that daily variation in flow may be more important to WT than seasonal variation.Indeed, in case of our new model, the θ parameters linking the seasonal component of WT to that of flow on the Scorff river appeared centered on 0, meaning that they may not be statistically significant.This appeared surprising as AT and Q six months minima and maxima time series were almost identical between models (see S1 File).In addition, the correlations between Q and AT time series entering the set of Eq (2) remained rather weak, suggesting that the signals were contrasted between covariates.Correlation were even slightly decreasing when going from model M 1 to M 2 , changing from 0.48 to 0.4 for min AT and max Q , and from -0.39 to -0.31 for max AT and min Q .It would be interesting to see if this behavior is also observed on other rivers.Still, the impact of daily flow variations on WT has practical consequences.It indicates that daily management of flow to avoid extremely low and high flow might be useful in the context of climate change adaptation [47].In particular, ensuring adequate flow in summer will serve to alleviate thermal stress, as flow can mitigate warming.The main issue is that both rainfall predictions by climate models and the subsequent rainfallrunoff models may lack enough precision [48].
The results from model M 2 also indicated that WT warming should track that of AT closer than predicted with model M 1 .In case of the Scorff river, WT potential warming appeared to be at the same scale as that of the AT warming scenario.On the Rough River, the WT forecast appeared less likely to reach 5˚C and results suggested that it might be overall slower than on the Scorff river.The more pronounced changes in predictions observed on the Scorff River with our updated model is likely linked to two key points.First it appears that modeling the residuals as done with model M 2 allows for a more precise estimate of the long-term component of WT time series than when using model M 1 .The AR1 structure used in M 1 probably absorbed a part of the signal in 6-month values.The Scorff river WT time series were also relatively short and prone to a lot of missing values.Those two things together probably led to an even greater parameter confusion while fitting the signal decomposition.Thus, we observed even more uncertainty in the 6 months min and max estimates from model M 1 on the Scorff river compared to those derived from model M 2 .Lastly, a probable more important WT warming on the Scorff river is also congruent with the energy budget theory [31].The Scorff river is much longer than the Rough river (75 and 12.3 km respectively).As such, thermal exchanges between air and water are more important and some homogenization is to be expected.Headwaters may thus be way less impacted by climate change as observed on the Rough river and in the USA [49] and their role in mitigating the global warming impact could be important.We note, however, that the opposite observation about the response of headwaters to warming has also been made [50] and warming predictions still overlap significantly in our cases.Further analyses on a range of rivers is necessary to better understand the role of different types of rivers in mitigating climate change.
The overall uncertainty associated with WT warming is high.This will have to be accounted for in decision making.Indeed, whatever the model, the 95% credible interval in WT warming was about 2˚C around a given temperature.IPCC scenarios highlight uncertainties in AT warming in 2100 close to 3.5˚C [46].The combined uncertainty is very high even over the medium term.Integrating it in a probabilistic framework is important for future ecological studies and decisions making [32].The approach we developed makes this possible.
Further extensions of our model are possible.now, the model is limited to one station per river, but may be extended to several.This has already been done for instance with models based on linear [50] and nonlinear regressions [51].The application could also encompass several rivers as done for instance with regression or neural network [52].Such an approach could include a hierarchical structure to test for an increased link between AT and WT as the stations are further downstream.The number and length of the time series would have to be chosen carefully as Bayesian fit can be slow when dealing with a large amount of data.Other covariates may also be included in further research.Shading by riparian vegetation can significantly alter WT and vegetation may change due to climate or land use management decisions [24,53].The shading factor could be included within the modelling of residuals as reducing extreme values through an additional exponential multiplier.
We offer stream ecologists and stakeholders a generic, parsimonious and effective way to complete, hindcast and predict WT time series in a temperate context.Using only AT and optionally Q time series, it predicts daily WT.The Bayesian framework allows for a full propagation of uncertainties while doing so.A probabilistic rational can then be used for management decision.The code is available on https://github.com/GuillaumeBal/2023.bayes.stream.temperature so readers can perform their own fits.Providing correctly formatted data (locally observed WT, AT and optionally Q data) and the specification a few model options is all that is required to run both models M 1 and M 2 .

Fig 2 .
Fig 2. WT α and β parameters posteriors for model fit.In case of model M 2 WT residuals are linked to those of AT and Q. https://doi.org/10.1371/journal.pone.0291239.g002

Fig 3 .
Fig 3. Boxplot of discrepancies between observed and fitted mean six months temperature.Model M 1 features an AR1 modeling of WT residuals.In case of model M 2 , WT residuals are linked to those of AT and optionally Q. https://doi.org/10.1371/journal.pone.0291239.g003

Fig 4 .
Fig 4. Observed versus forecasted temperatures coming from the cross validation.Yellow lines correspond to observed WT temperatures.Blue lines are the median value of predictive posteriors.Grey areas indicate 95% and 50% intervals of the forecasts.https://doi.org/10.1371/journal.pone.0291239.g004

Fig 5 .
Fig 5. Comparison of the posteriors of θ parameters from equation set (2) updated with both model M 1 and M 2 .The left-hand side of the figure are for parameters linking max WT y to max AT y and optionaly min Q y .The left-hand side of the figure are for parameters linking min WT y to min AT y and optionaly max Q y .https://doi.org/10.1371/journal.pone.0291239.g005

Fig 7 .
Fig 7. Warming of stream water temperature predicted across a range of air temperature warming.The dotted blue line is the 1:1 line.https://doi.org/10.1371/journal.pone.0291239.g007 srv/eng/catalog.search#/metadata/ie.marine.data:dataset.5122.In case of the Scorff river, data are accessible free of charge as well.Water temperature are accessible through a DOI too.Flow data are hosted onto the French hydrology services database and are downloadable at https://hydro.

Table 2 . Fitting and forecasting performances summary statistics. River Model RMSE χ 2 p-value Deviance pD DIC
for the posterior checking tests; RMSE: root mean square errors used to quantify the predictive performance.https://doi.org/10.1371/journal.pone.0291239.t002