Restored off-channel pond habitats create thermal regime diversity and refuges within a Mediterranean-climate watershed

Cool-water habitats provide increasingly vital refuges for cold-water fish living on the margins of their historical ranges; consequently, efforts to enhance or create cool-water habitat are becoming a major focus of river restoration practices. However, the effectiveness of restoration projects for providing thermal refuge and creating diverse temperature regimes at the watershed scale remains unclear. In the Klamath River in Northern California, the Karuk Tribe Fisheries Program, the Mid-Klamath Watershed Council, and the U.S. Forest Service constructed a series of off-channel ponds that recreate floodplain habitat and support juvenile coho salmon (Oncorhynchus kisutch) and steelhead (Oncorhynchus mykiss) along the Klamath River and its tributaries. We instrumented these ponds and applied multivariate auto-regressive time series models of fine-scale temperature data from ponds, tributaries, and the mainstem Klamath River to assess how off-channel ponds contributed to thermal regime diversity and thermal refuge habitat in the Klamath riverscape. Our analysis demonstrated that ponds provide diverse thermal habitats that are significantly cooler than creek or mainstem river habitats, even during severe drought. Wavelet analysis of long-term (10 years) temperature data indicated that thermal buffering (i.e. dampening of diel variation) increased over time but was disrupted by drought conditions in 2021. Our analysis demonstrates that in certain situations, human-made off-channel ponds can increase thermal diversity in modified riverscapes even during drought conditions, potentially benefiting floodplain-dependent cold-water species. Restoration actions that create and maintain thermal regime diversity and thermal refuges will become an essential tool to conserve biodiversity in climate-sensitive watersheds. 

Methods
Data Collection
In July 2020, we deployed 30 temperature sensors (HOBO MX2201, Onset Corporation, Massachusetts) programmed to measure temperature every 15 minutes in ponds and creeks. We placed 1-4 sensors in each pond to capture local-scale temperature variation. Sensors were installed at approximately ⅓ the water depth (at time of placement), except for two sensors in Goodman Pond, one in Upper Lawrence Pond, and one in Lower Lawrence Pond, where sensors were placed on the bottom of the pond. We chose these locations to capture within-pond variation in thermal habitat, to maximize access and safety, and to facilitate future monitoring. We also placed 1 sensor in the creek upstream of the outlet of each pond. We placed sensors between 7 and 13 July 2020 and read them out between 11 and 13 July 2021. We removed incomplete sensor time series (n = 6 pond sensors and n = 3 creek sensors) resulting either from sensor malfunction or sensors that were no longer submerged because of drought-related decreases in water level. In ponds and creeks with multiple sensors, we averaged the remaining sensor readings to obtain an average time series per site. In the 5 ponds with only one sensor, we used that sensor’s time series. We calculated and modeled daily temperature means (instead of using sub-daily data) to avoid having to account for diel periodicity in the MAR models (Hampton et al. 2013; Holmes et al. 2023), which would have made these models unnecessarily complex.

Water levels in the pond fluctuated throughout the year, leading to different depths for the sensors throughout the study period, which could influence temperatures. We removed from analysis sensors that were completely out of the water (thus, recording air temperature rather than water temperature) because of depth fluctuations. To understand how well the remaining sensors represent thermal habitats in the ponds, we took post hoc temperature-depth profiles in June 2023 at several locations in each pond. 

Klamath River temperature data were collected by the Karuk Tribe and accessed with permission from the Karuk Tribe Water Quality Department (Accessed 27 September 2022). We used data between May 2020 and February 2021. We used a combination of data readings from the Seiad Valley station as well as interpolated data using a linear regression from the Orleans station when Seiad Valley data was unavailable (5.5% of Seiad Valley data was interpolated). Additionally, we obtained air temperature time series from the National Oceanic and Atmospheric Administration’s Climate Data Online database for Siskiyou County, CA (NOAA 2020). We used the Slater Butte air sensor, located relatively close (13 km) to our study sites in Seiad Creek. Although these two sites differ in elevation (1423 m vs. 430 m), we expected fluctuations in air temperature at these two locations to be correlated, and we note that our models quantify the effects of fluctuations around the mean rather than absolute values of air temperature (see next section). Also, we measured DO in a single location in each pond over several days in July 2020 (Figure S3), and we took post hoc DO and temperature-depth profiles in each pond in June 2023 (Figure S2). 

Thermal diversity
To analyze variation in thermal regimes across the riverscape, we used multivariate auto-regressive model (MAR) models. The MAR model is a time series model that takes advantage of temporal correlation in environmental variables to estimate the effects of a particular driver, while also accounting for stochastic process error (Holmes et al. 2014; Ives et al. 2003; Ruhí et al. 2015). A MAR model in the matrix form can be expressed as follows:

Xt = BXt-1 +  Cct-1 + wt,  where wt ~ MVN(0, Q) 

where temperature at a given day (Xt) is a function of temperature the previous day (Xt-1) plus sensitivity to a covariate, here variation in air temperature (Cct-1); and process error (wt). As a covariate (ct-1), we used a time series of air temperature with a one-day time lag, after examining support for other lags (results not shown); and the C matrix captured site-specific sensitivity to air temperature. In turn, process error (wt) was drawn from a multivariate normal distribution, with mean zero and covariance matrix Q. In our case, Q captured stochasticity in water temperature (i.e., temporal variation in water temperature that was unrelated to air temperature). B is an interaction matrix that can model the effect of each state on itself (diagonal parameters) and on each other (off-diagonal parameters). In our case, we set off-diagonal parameters to zero (as we did not expect sites to interact with each other) and estimated the diagonal parameters, often used to capture “density-dependence” in population processes, or pull-back to mean. 

To test our prediction that off-channel ponds have significantly different thermal regimes compared to creeks and the mainstem, we developed four MAR model hypotheses that represent different levels of complexity in thermal regimes (as in Leathers et al. 2022). Each hypothesis was tested by manipulating the matrices of the MAR model, capturing stochastic or ‘unexplained’ variation (Q matrix), and deterministic or covariate-explained variation (C matrix). This strategy allowed modeling mean daily temperatures among pond, creek, and river habitats in different ways. The first hypothesis was that all sites had different levels of stochastic and deterministic variability (i.e., as many thermal regimes as sites). The second hypothesis was that each habitat type (pond, creek, and river) had some typical level of stochastic and deterministic variability, but sites within the same habitat type did not differ from each other. The third hypothesis predicted that stochastic and deterministic variability depended on the watershed (Horse Creek vs. Seiad Creek vs. Klamath River), but not the specific site or habitat type. The fourth hypothesis predicted that all sites would have the same level of stochastic and deterministic variability (i.e., a single, watershed-level thermal regime). We used Akaike’s Information Criterion corrected for small sample size (AICc) to compare support for the different hypotheses. All data and covariate data was z-scored, and model outputs were examined for normality and autocorrelation of residuals via the Autocorrelation Function (ACF). We used the MARSS package version 3.11.3 (Holmes et al. 2021) in R (R Development Core Team 2021).

Thermal buffering
To quantify thermal buffering of ponds (relative to creeks), we compared daily maximum temperatures (averaged across all sensors in a site) in each pond and creek during the three hottest months of 2020 (15 July – 15 September), and then ran a one-way ANOVA of temperature ~ site. We repeated the same process for the winter, focusing on daily minimum temperatures during the three coldest months (15 December 2020 – 15 February 2021). We ensured that model residuals met assumptions of normality and homogeneity of variances.

We also assessed daily thermal buffering capacity of ponds and creeks by calculating the coefficient of variation (CV) for each day, using 15-minute temperature data. We then averaged daily CVs for each site over the yearlong study period. We used mean CV values to calculate the ratio of creek to pond CV for each pond/tributary pairing. If the creek:pond CV ratio was equal or less than 1, that suggested no significant buffering took place. If the ratio was greater than 1, we considered the pond to “buffer” thermal fluctuations compared to the creek.

Thermal stabilization over time
We used wavelet analysis to examine thermal regimes in the frequency and time domains and to determine whether some scales of variation strengthened over time. Although the wavelet method does not require pre-specifying a frequency of interest, here we focused on temperature variation at diel (24-hour) and seasonal scales (12-months), and asked whether diel and seasonal variation changed over the years. We interpolated missing values in the historical temperature datasets for Alexander and Stender Ponds (3.3% and 3.9% of days, respectively) via an autoregressive integrated moving average model (ARIMA) and a Kalman filter. An ARIMA model is generally expressed as ARIMA(p, d, q), where p is the order of the autoregressive model, i.e. the dependence of the model on prior values; d is the order of non-seasonal differences, i.e. degree of differencing of raw observations; and q is the order of the moving average, i.e. the model’s dependence on longer-term values and stochastic “shocks”. After identifying the best-fit ARIMA model, we used the Kalman filter to interpolate missing data. We then ran wavelets on the complete time series, using the WaveletComp package in R (Roesch & Schmidbauer 2018). We used the Morlet wavelet function and compared observed power to a null background generated with red noise (i.e. temporally autocorrelated data).