A hierarchical model of daily stream temperature using air-water temperature synchronization, autocorrelation, and time lags

Descriptive statistics

As a coarse comparison of daily water temperatures, we calculated correlations among sites. We also explored patterns in water temperature over time and among sites by comparing cumulative residuals from a spline fit to all the data (function gam() in R, Fig. 2). We calculated residuals for each water temperature data point and then developed empirical cumulative curves over days of the year for each year and site combination.

Breakpoints

The goal is to develop a robust model for the relationship between mean daily water and mean daily air temperature. A key limitation in developing this relationship is that lower water temperatures in the winter are bounded near 0 °C while air temperatures are not. This means that water and air temperatures can become decoupled when air temperatures are cold resulting in only a weak relationship, at best, between water and air temperature. In contrast, as air temperatures warm in the spring and before they get too cold in the autumn, daily water and air temperatures can be synchronized (Figs. 3A and 3B), suggesting the possibility of a strong relationship between daily water and air temperature during the synchronized portion of each year.

The key to the approach is identifying a breakpoint in the spring when daily water and daily air temperature become synchronized and a breakpoint in the autumn when daily temperatures become desynchronized. To identify the synchronization breakpoints we calculated a simple index (waterT−airT)/waterT (waterT > 0), where waterT was mean daily water temperature and airT was mean daily air temperature (Fig. 3). This temperature index (tempIndex) approaches 0 when water and air temperature are similar and is very different from 0 when temperatures diverge (Figs. 3C and 3D). While water and air temperatures are synchronized, tempIndex flattens out (Figs. 3C and 3D), providing the opportunity to identify the beginning and end (breakpoints) of the flat period.

To identify the spring and autumn breakpoints, we used a runs analysis that determined the first (spring) and last (autumn) day of the year that the tempIndex was consistently within the flat period (Figs. 3C and 3D). We established the range of tempIndex values that comprised the flat period by calculating the 99.9% confidence interval (CI) for tempIndex using the middle 150 days of the year (late April to mid-September). The middle 150 days of the year were always within the flat period based on visual observation of tempIndex plots. Separate CI values were calculated for each year and stream. For the breakpoint estimation, we used a moving average for tempIndex with a centered 10-day window to help stabilize tempIndex values near the breakpoints. Temperatures were considered synchronized when 10 consecutive days of the moving average fell within the 99.9% CI. Beginning on day 1 and moving towards day 150, the first time 10 consecutive days were synchronized was used as the spring breakpoint and we moved from the end of the year to day 150 to establish the fall breakpoint. Numbers of days in the synchronized period for each stream and year are shown in Table 2.

We evaluated trends in fall and spring breakpoints by running three linear models with breakpoint day of the year as the dependent variable and year alone or year + stream or year * stream as independent variables. We estimated AIC to determine the most parsimonious model.

Water temperature model description

With breakpoints established for each year and site, we modeled the relationship between water temperature and air temperature for the synchronized period using a hierarchical linear autoregressive model with a cubic trend across days within a year and covariation among sites. We fit the model using a Bayesian approach.

Observed water temperature (t_s,d,y) for each site (s; s₁ = WB, s₂ = OL, s₃ = OS, s₄ = IL), day of year (d) and year (y) was assumed to derive from a normal distribution with mean μ_s,d,y and standard deviation sd (residual model error): (1)${\displaystyle t_{s,d,y}\sim N\left({\mu }_{s,d,y},sd\right).}$ We used a non-informative uniform prior [0, 10] for sd. We modeled the mean with a linear trend (ω_s,d,y) adjusted by an AR(1) autoregressive coefficient (δ_s) on the residual error from the previous day: (2)${\displaystyle {\mu }_{s,d,y}={\omega }_{s,d,y}+{\delta }_s\left(t_{s,d-1,y}-{\omega }_{s,d-1,y}\right).}$ We placed a hierarchical structure on δ_s: (3)${\displaystyle {\delta }_s\sim N\left({\mu \delta }_s,sd\delta \right)T\left(-1,1\right)}$ where site-specific δ_s were drawn from a truncated normal distribution with mean μδ_s and standard deviation sdδ. Values for δ_s were truncated to keep them within the admissible range for a correlation. Priors for the mean and standard deviation were non-informative; μδ_s ∼ U(−1, 1), and sdδ ∼ U(0, 2) (an upper limit of 2 for sdδ is non-informative for the truncated data).

When observed temperature data were not available for the previous day (beginning of a series or following a break in the series) we modeled the mean without the autoregressive component: (4)${\displaystyle {\mu }_{s,d,y}={\omega }_{s,d,y}.}$ We modeled the linear component with a combination of fixed and random effects: (5)${\displaystyle {\omega }_{s,d,y}=\alpha +{{\beta }_{1}T}_{s,d,y}+{{\beta }_{2}T}_{s,d-1,y}+{{\beta }_{3}T}_{s,d-2,y}+{\beta }_4{F}_{s,d,y}+{\beta }_5{T}_{s,d,y}\cdot F_{s,d,y}+{\beta }_{6:8}s+{\beta }_{9:11}s\cdot T_{s,d,y}+Y_y}$ where α is the overall intercept, the β are the coefficients for the fixed effects (T is mean daily air temperature, F is mean daily stream flow, s is site) and Y_y represents random effects among years. Priors for the β_1:11 were independent and non-informative, N(0, 100).

Y_y represented random effect temporal trends (cubic) across years where: (6)${\displaystyle Y_y={\alpha }_y+{{\beta }_{12,y}D}_{s,y,d}+{\beta }_{13,y}{D}_{s,y,d}^2+{\beta }_{14,y}{D}_{s,y,d}^3.}$ For convenience, this equation can be written in matrix notation as (7)${\displaystyle Y_y={\mathrm{B}}_y{X}_{s,d,y}}$ where X is a data matrix with l columns (l = 4; the number of year-level predictors) with the first column a vector of 1’s for the intercept and B_y is the y × l matrix of year-level regression coefficients. Priors for the mean were non-informative, with (8)${\displaystyle {\mathrm{B}}_y\sim MVN\left(M_l,\Sigma \right)}$ where M_l = (μ_α, μ_β12, μ_β13, μ_β14) is a vector of length l, representing the mean of the distribution of intercept and slopes. The l × l covariance matrix is represented by Σ where the variance of each regression coefficient is on the diagonal and the covariance on the off-diagonals. The hyperprior for the means were non-informative with μ_α, = 0 and μ_β12, μ_β13, μ_β14 ∼ N(0, 100). Standard deviation priors were also non-informative and were drawn from an inv-Wishart distribution: (9)${\displaystyle \Sigma \sim \mathrm{inv}-\mathrm{wish}\left(\mathrm{diag}\left(l\right)l+1\right).}$

Goodness of fit and prediction

We assessed goodness of fit in two ways. First, we compared observed and predicted values for the complete dataset. Second, we ran a series of cross validation tests where we randomly left out a portion of the water temperature data, estimated parameters with the remaining (training) data and compared predictions of the left out (testing) data to original values. This involved leave-p-out cross-validation where we randomly left out a proportion (p) of the data, where p = 0, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8. We ran 10 replicates for each value of p. For each condition, we calculated the root mean square error (RMSE) of the residuals for the training and the test data sets.