Spatial Stream Network Analysis of Nechako Watershed Stream Temperatures 2022b

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. 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 4-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} = \frac{1}{(\theta_{s,j})^{a4_s}}(a1_s + a2_s A_{s,j} - a3_s W_{s,j - 1}), \end{equation}\]

where \(a1_s\), \(a2_s\), \(a3_s\), and \(a4_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, \(W_{s, j - 1}\) is the expected stream temperature at the \(s^{th}\) site in the previous week, and \(\theta_{s,j}\) is the dimensionless discharge for the \(s^{th}\) site in the \(j^{th}\) week. \(\theta_{s,j}\) was calculated as follows:

\[\begin{equation} \theta_{s,j} = \frac{d_{s,j}}{\bar{d_s}} \end{equation}\]

where \(d_{s,j}\) is the discharge for the \(s^{th}\) site in the \(j^{th}\) week, and \(\bar{d_s} = \frac{\sum_{j = 1}^{J}(d_{s,j})}{J}\) is the mean discharge across all \(J\) weeks for the \(s^{th}\) site.

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\), \(a3\), and \(a4\)) 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:

• The exponential tail-down model was better at explaining the spatial correlation in the data than exponential tail-up or euclidean distance models (Ver Hoef and Peterson 2010; Peterson and Hoef 2010).
• The full 8-parameter air2stream model did not converge.
• 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
  int <lower=1> i_y_obs [N_y_obs] ;  // [N_y_obs,T]
  int <lower=1> i_y_mis [N_y_mis] ;  // [N_y_mis,T]
  vector [N_y_obs] y_obs;  // matrix[N_y_obs,1] y_obs[T];
  real discharge [nsite * nweek];
  real air_temp [nsite * nweek];
  int<lower=0> site[nsite * nweek];
  int<lower=0> week[nsite * nweek];

  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=0> s4;
  real<lower=-5, upper=15> m1;
  real<lower=-5, upper=1.5> m2;
  real<lower=-5, upper=5> m3; 
  real<lower=-1, upper=1> m4;
  real<lower=-5, upper=15> a1[nsite];
  real<lower=-5, upper=1.5> a2[nsite]; 
  real<lower=-5, upper=5> a3[nsite]; 
  real<lower=-1, upper=1> a4[nsite];

transformed parameters {
  vector[nsite * nweek] y; // long vector of y
  vector[nsite] Y[nweek]; // array of 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;
  vector[nsite] mu [nweek];
  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] = (1/(discharge[i]^a4[site[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(0, 20000) T[0, ]; // range tail-down
  bInitialTemp ~ normal(0, 0.1) T[0, ];

  s1 ~ exponential(50);
  s2 ~ exponential(50);
  s3 ~ exponential(50);
  s4 ~ exponential(50);
  m1 ~ normal(0.8, 1);
  m2 ~ normal(0.4, 1);
  m3 ~ normal(0.4, 1);
  m4 ~ normal(0.1, 1);
  a1 ~ normal(m1, s1);
  a2 ~ normal(m2, s2);
  a3 ~ normal(m3, s3);
  a4 ~ normal(m4, s4);

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

Block 1. Model description.

Stream Temperature

Table 1. Parameter descriptions.

Table 2. Model coefficients.

Table 3. Model convergence.

Table 4. Model sensitivity.

Stream Temperature