Spatial Stream Network Analysis of Skeena Watershed Stream Temperatures 2025

Stream Temperature

The data were analyzed using a Spatial Stream Network model (Ver Hoef and Peterson 2010; Peterson and Hoef 2010), with code adapted from the SSNbayes package (Santos-Fernandez et al. 2022). The necessary stream network distances and connectivity were calculated using the BC Freshwater Atlas and the SSNbler (Peterson et al. 2024) and SSN2 (Dumelle et al. 2024) R packages. Air and stream temperature data were averaged by site and week; modeling was done on this weekly time scale.

The expected stream temperatures were modeled using the 3-parameter version of the air2stream model (Toffolon and Piccolroaz 2015). The average stream temperature (in ˚C) in the first week, \(W_{s,j=1}\) was estimated by the model and assumed to be the same for all sites.

For all subsequent weeks (i.e., \(j > 1\)), the change in the stream temperature (in ˚C) between week \(j - 1\) and week \(j\) for the \(s^{th}\) site, \(\Delta W_{s,j}\), was modeled as follows:

\[\begin{equation} \Delta W_{s,j} = a1_s + a2_s A_{s,j} - a3_s W_{s,j - 1}, \end{equation}\]

where \(a1_s\), \(a2_s\), and \(a3_s\) are the parameters of the air2stream model for the \(s^{th}\) site, \(A_{s,j}\) is the air temperature (in ˚C) for the \(s^{th}\) site in the \(j^{th}\) week and \(W_{s, j - 1}\) is the expected stream temperature at the \(s^{th}\) site in the previous week.

The expected stream temperature for the \(s^{th}\) site in the \(j^{th}\) week was then calculated:

\[\begin{equation} W_{s,j} = W_{s,j - 1} + \Delta W_{s,j}. \end{equation}\]

Growing Season Degree Days (GSDD) are the accumulated thermal units (in ˚C) during the growing season based on the mean daily water temperature values, which is a useful predictor of age-0 rainbow and westslope cutthroat trout size at the beginning of winter. The start and end of the growing season were based on the definitions of Coleman and Fausch (2007):

• Start: the beginning of the first week that average stream temperatures exceeded and remained above 5˚C for the season.
• End: the last day of the first week that average stream temperature dropped below 4˚C.

GSDD were derived for each site and year by assuming that the daily stream temperatures at each site were the predicted weekly mean stream temperature for every day in the given week.

Key assumptions of the model include:

• The stream network is dendritic, not braided.
• The expected stream temperatures were set to 0˚C if they were estimated to be negative by the model.
• The stream temperature in the first week is the same for all sites.
• The parameters of the air2stream model (\(a1\), \(a2\), and \(a3\)) vary randomly by site.
• The residual variation is multivariate normally distributed.
• The covariance structure of the residual variation combines the following covariance components:
  • Nugget (allows for variation at a single location)
  • Exponential tail-down (allows for spatial dependence between flow-connected and flow-unconnected locations)

Preliminary analysis found that:

• Allowing the initial stream temperature to vary randomly by site produced unrealistic stream temperatures for January (> 10˚C).

Stream Temperature

data {
  int nsite;
  int nweek;

  int <lower=0> N_y_obs; // number observed values
  int <lower=0> N_y_mis; // number missing values
  array[N_y_obs] int<lower=1> i_y_obs;
  array[N_y_mis] int<lower=1> i_y_mis;
  vector [N_y_obs] y_obs;  // matrix[N_y_obs,1] y_obs[T];
  array[nsite * nweek] real air_temp;
  array[nsite * nweek] int<lower=0> site;
  array[nsite * nweek] int<lower=0> week;

  matrix [nsite, nsite] D;
  matrix [nsite, nsite] I;
  matrix [nsite, nsite] H;
  matrix [nsite, nsite] flow_con_mat;

parameters {
  vector<lower=0, upper=30>[N_y_mis] y_mis; // declaring the missing y

  real<lower=0> sigma_nug; // sd of nugget effect
  real<lower=0> sigma_td; // sd of tail-down
  real<lower=0> alpha_td; // range of the tail-down model
  real<lower=0> bInitialTemp;
  real<lower=0> s1;
  real<lower=0> s2;
  real<lower=0> s3;
  real<lower=-5, upper=15> m1;
  real<lower=-5, upper=1.5> m2;
  real<lower=-5, upper=5> m3; 
  array[nsite] real<lower=-5, upper=15> a1;
  array[nsite] real<lower=-5, upper=1.5> a2;
  array[nsite] real<lower=-5, upper=5> a3;

transformed parameters {
  vector[nsite * nweek] y; // long vector of y
  array[nweek] vector[nsite] Y;
  matrix[nsite, nsite] C_td; // tail-down cov
  real <lower=0> var_nug; // nugget
  real <lower=0> var_td; // partial sill tail-down

  vector[nsite * nweek] eTempDiff;
  vector<lower=0, upper=30>[nsite * nweek] eTemp; 
  array[nweek] vector[nsite] mu;
  y[i_y_obs] = y_obs;
  y[i_y_mis] = y_mis;
  var_nug = sigma_nug^2; // variance nugget
  var_td = sigma_td^2; // variance tail-down
  // Place observations into matrices
  for (t in 1:nweek){
    Y[t] = y[((t - 1) * nsite + 1):(t * nsite)];
  }
  eTemp[1:nsite] = rep_vector(bInitialTemp, nsite);
  for (i in (nsite + 1):(nweek * nsite)) {
    eTempDiff[i] = (a1[site[i]] + a2[site[i]] * air_temp[i] - a3[site[i]] * eTemp[i - nsite]);
    
    eTemp[i] = eTemp[i - nsite] + eTempDiff[i];
    if (eTemp[i] < 0) {
      eTemp[i] = 0.0;
    }
  }
  // Define 1st mu
  mu[1] = eTemp[1:nsite];
  // Define rest of mu; ----
  for (t in 2:nweek){
    mu[t] = eTemp[((t - 1) * nsite + 1):(t * nsite)];
  }
  // Covariance matrices ----
  // Tail-down exponential model
    for (i in 1:nsite) {
    for (j in 1:nsite) {
      if (flow_con_mat[i, j] == 1) { // if points are flow connected
        C_td[i, j] = var_td * exp(-3 * H[i, j] / alpha_td);
      }
      else{ // if points are flow unconnected
        C_td[i, j] = var_td * exp(-3 * (D[i, j] + D[j, i]) / alpha_td);
      }
    }
  }

model {
  sigma_nug ~ exponential(0.05); // sd nugget
  sigma_td ~ exponential(2); // sd tail-down
  alpha_td ~ normal(20000, 20000); // range tail-down
  bInitialTemp ~ normal(1, 1);

  s1 ~ exponential(1);
  s2 ~ exponential(1);
  s3 ~ exponential(1);
  m1 ~ normal(0.8, 1);
  m2 ~ normal(0.4, 1);
  m3 ~ normal(0.4, 1);
  a1 ~ normal(m1, s1);
  a2 ~ normal(m2, s2);
  a3 ~ normal(m3, s3);

  for (t in 1:nweek) {
    target += multi_normal_cholesky_lpdf(Y[t] | mu[t], cholesky_decompose(C_td + var_nug * I + 1e-6));
  }

generated quantities {
  array[nweek * nsite] real ePrediction;
  array[nweek] real log_lik;
  real lprior;
  for (t in 1:nweek) {
    log_lik[t] = multi_normal_cholesky_lpdf(Y[t] | mu[t], cholesky_decompose(C_td + var_nug * I + 1e-6));
    ePrediction[((t - 1) * nsite + 1):(t * nsite)] = to_array_1d(multi_normal_cholesky_rng(mu[t], cholesky_decompose(C_td + var_nug * I + 1e-6)));
  }
  lprior = exponential_lpdf(sigma_nug | 0.05) +
           exponential_lpdf(sigma_td | 2) +
           (normal_lpdf(alpha_td | 20000, 20000) - normal_lccdf(0 | 20000, 20000)) + // truncated normal
           (normal_lpdf(bInitialTemp | 1, 1) - normal_lccdf(0 | 1, 1)) + // truncated normal
           normal_lpdf(m1 | 0.8, 1) +
           normal_lpdf(m2 | 0.4, 1) +
           normal_lpdf(m3 | 0.4, 1) +
           exponential_lpdf(s1 | 1) +
           exponential_lpdf(s2 | 1) +
           exponential_lpdf(s3 | 1);

Block 1. Model description.

Locations

Table 1. Stream temperature site locations.

Stream Temperature

Table 2. Parameter descriptions.

Table 3. Model coefficients.

Table 4. Model convergence.

Table 5. Model sensitivity.

Table 6. Model sensitivity by parameter.

Table 7. Mean GSDD by year (with 95% prediction intervals).

Table 8. Mean GSDD by site (with 95% prediction intervals).

Stream Temperature