Escape from the heat: thermal stratification in a well-mixed estuary and implications for fish species facing a changing climate

Climate change may cause organisms to seek thermal refuge from rising temperatures, either by shifting their ranges or seeking microrefugia within their existing ranges. We evaluate the potential for thermal stratification to provide refuge for two fish species in the San Francisco Estuary (Estuary): Chinook Salmon (Oncorhynchus tshawytscha Walbaum, 1792) and Delta Smelt (Hypomesus transpacificus McAllister, 1963). We compiled water temperature data from multiple monitoring programs to evaluate spatial, daily, hourly, intra-annual, and inter-annual trends in stratification using generalized additive models. We used our models to predict the locations and periods of time that the bottom of the water column could function as thermal refuge for salmon and smelt. Periods in which the bottom was cooler than surface primarily occurred during the peak of summer and during the afternoons, with more prominent stratification during warmer years. Although the Estuary is often exceedingly warm for fish species and well-mixed overall, we identified potential thermal refugia in a long and deep terminal channel for Delta Smelt, and in the periods bordering summer for Chinook Salmon. Thermal stratification may increase as the climate warms, and pockets of cooler water at depth, though limited, may become more important for at-risk fishes in the future.


2011
). Estuaries are also some of the most degraded ecosystems on the planet: tributary rivers have been dammed and diverted, food webs have been altered by invasive species, and contaminants have been introduced into the waters (Lotze et al., 2006). Many of these anthropogenic impacts are exemplified in the San Francisco Estuary of California, United States (Nichols et al., 1986;Brown & Bauer, 2010). The San Francisco Estuary has lost a large portion of its historical tidal wetlands due to diking, freshwater flow is largely controlled by water control infrastructure, and multiple native fish species are in severe decline (Kimmerer, 2004). Some species, such as the San Francisco Estuary-endemic Delta Smelt (Hypomesus transpacificus McAllister, 1963) and the anadromous Chinook Salmon (Oncorhynchus tshawytscha Walbaum, 1792) are of particular interest due to their listing under the federal and California Endangered Species Acts (Perry et al., 2016;Moyle et al., 2018;NMFS, 2019;USFWS, 2019).
Water temperature is a key variable controlling the physiology, behavior, and distribution of fishes. Temperature can influence the growth rate, metabolism, and swimming activity of individual fish (Jobling, 1997;Breau et al., 2011;Jeffries et al., 2016;Davis et al., 2019a). These effects can impact fish species at a population level through increased predation-associated mortality or reduced spawning success (Rose et al., 2013;Davis et al., 2019b;Michel et al., 2020). Shifts in water temperature can also alter the routing and timing of migration and, therefore, can change the distribution of fishes at large geographical scales (Munsch et al., 2019;Goertler et al., 2021). Fish can exploit the spatial and temporal temperature variation within a system to optimize their bioenergetics (Bevelhimer & Adams, 1993;Neverman & Wurtsbaugh, 1994;. They can also seek refuge (whether cooler or warmer) in order to avoid thermal stress or potentially lethal conditions (Torgersen et al., 1999).
Climate change is forecasted to increase water temperatures through atmospheric warming and in California, this will be exacerbated by declining snowpack in the wet season (i.e., winter) (Dettinger & Cayan, 1995). In response to warming, marine fishes have been migrating towards higher latitudes (Hickling et al., 2006) while freshwater fishes have been migrating to higher elevations to seek more suitable thermal conditions (Comte & Grenouillet, 2013).
However, some species with limited opportunities for dispersal (e.g., spring and lake-dwelling fishes, estuary-specialist species with narrow abiotic tolerances) or anadromous species that rely on thermally suitable migration corridors may have few options for avoiding warming conditions (Ray, 2005). Evidence suggests that water temperatures are already reaching stressful levels for thermally sensitive fishes in some areas of the upper San Francisco Estuary and that climate change will continue to increase water temperature over time (Brown et al., 2013Nobriga et al., 2021). Although modeling (Vroom et al., 2017) and observational  studies suggest that significant vertical differences in water temperature are likely uncommon in the upper San Francisco Estuary, the hypothesis that such thermal stratification can act as temporary refuge for aquatic species has not been well-tested. Furthermore, assessment of thermal heterogeneity in the landscape can help better understand the ecology and evolution of the biota of any system.
Here we seek to evaluate existing discrete (i.e., boat-based surveys) and continuous (i.e., sondes) records of water temperature in the San Francisco Estuary to determine the timing, frequency, duration, and magnitude of differences between near-surface and near-bottom water temperature. We compiled thousands of discrete temperature data points that have been collected since 2011 and data from four continuous monitoring stations that have been collected since 2012. These data were interpolated over space and time using empirical models to identify patterns of thermal stratification. Our study questions are as follows: 1) Where on the landscape does thermal stratification occur?, 2) When during the day and over the course of the year does thermal stratification occur?, and 3) Can thermal stratification provide refuge for fish species of concern? We specifically consider crucial life stages for two species of management concern in the San Francisco Estuary: subadult Delta Smelt and juvenile outmigrant Chinook Salmon. Delta Smelt is an example of an estuarine resident species with limited capability to disperse to other estuaries (Sommer & Mejia, 2013), while Chinook Salmon is an example of an anadromous species that relies on the San Francisco Estuary for only a part of its life cycle (Perry et al., 2016) (Fig. 1). Patterns of thermal stratification and its potential impacts to fish species as described in this study can be instructive for other major estuaries around the world inhabited by thermally sensitive and estuarydependent fishes.

Study system
The San Francisco Estuary (Estuary) is the largest estuary on the Pacific Coast of the United States, stretching from the tidal saline San Francisco Bay to the tidal freshwater Sacramento-San Joaquin Delta (Delta) (Fig. 2). This system has been highly altered by habitat modification, construction of water conveyance infrastructure, and species invasions (Nichols et al., 1986). Once composed of largely dendritic tidal wetlands and seasonal floodplains, the Delta of today is a network of diked channels bordered by levees built to protect agricultural tracts (Whipple et al., 2012). The Delta is managed to be mostly freshwater yearround, with some brackish water (1-2 parts per thousand salinity) intruding in the western Delta during the summer and fall of dry years. Salinity is managed through upstream reservoir releases, as well as adjustments to water exports at the pumping facilities in the southwestern Delta near Old and Middle rivers. The Estuary downstream of the Delta is dynamic in salinity and the low salinity zone between 1 and 6 ppt is of particular interest to managers due to the association of the endangered Delta Smelt with this habitat (Sommer & Mejia, 2013). The Estuary has a Mediterranean climate with warm and dry conditions in the summer and fall (June-October) and cold and wet conditions in winter and spring (November-May). The amount of precipitation and flow in the Estuary Fig. 1 Simplified life cycle and associated presumed thermal maxima of Chinook Salmon and Delta Smelt in the San Francisco Estuary. Note that both species demonstrate a diversity of life history types and the life cycles depicted here are broad generalizations for illustrative purposes. For Chinook Salmon, we focused on the dominant races (spring and fall-run). Addi-tionally, thermal maxima listed above were gathered from several studies and generalized to aid the interpretation of results and not meant to be definitive (Marine & Cech, 2004;Myrick & Cech, 2004;Nobriga et al., 2008Nobriga et al., , 2021Marston et al., 2012;Komoroske et al., 2014;Damon et al., 2016;Lewis et al., 2021) varies considerably from year to year. This high interannual variability plays a significant societal and ecological role in California (Dettinger et al., 1998;Dettinger & Cayan, 2003). Temperature in the Estuary is typically highest in July and lowest in January.

Discrete dataset
We integrated water temperature data collected at discrete time points in the Estuary from Regions are shown by dark blue outlines and grey outline denotes the legal Delta boundary. Note that region cutoffs were used only for data processing and visualization purposes (not for analyses) various boat-based monitoring programs and studies (Bashevkin, 2021). Further details on the data integration process and access to the data can be found in Bashevkin (2021). For this study, we only retained data with temperature measurements at both the surface and near-bottom of the water column (typically 0.5 m below the surface and 0.5 m above the bottom). Data containing both surface and bottom temperature information has limited spatial coverage prior to 2011. Therefore, we only used data from 2011 to 2019. We removed a few data points which appear to have been entered incorrectly (i.e., 0 °C recorded, data with spatial coordinates that are outside of the Estuary's water boundaries, 10 °C gap between surface and bottom). We also excluded data from regions with limited sample size (i.e., regions with no station containing 25 or more sampling occasions; see Fig. 2 for regional boundaries). To reduce spatiotemporal autocorrelation in the dataset, in cases where multiple samples were taken at a single location in a day, only a single sample was retained (the data point closest to 12 P.M. was chosen). The final discrete temperature dataset contained 9,463 data points from seven different studies in the Estuary: the Fall Midwater Trawl survey (Stevens & Miller, 1983), the San Francisco Bay study (Armor & Herrgesell, 1985), the Summer Townet survey (Turner & Chadwick, 1972) (Bashevkin, 2021), and the U.S. Geological Survey's San Francisco Bay survey .

Continuous dataset
We integrated continuous surface and bottom water temperature data from four stations monitored by the California Department of Water Resources. Stations included: Antioch (ANH) in the San Joaquin River, Mallard Island (MAL) in the confluence area of the Sacramento and San Joaquin rivers, Martinez (MRZ) at Carquinez Strait, and Rough and Ready Island (RRI) in the San Joaquin River at the Port of Stockton (Fig. 2). All bottom sensors were approximately 1 m from the bottom of the riverbed, and all surface sensors were near the surface of the water column at a depth of approximately 3-4 m. Surface temperature data were filtered for relevant stations from IEP et al.
(2020). Bottom temperature data were standardized, and quality checked as described by IEP et al. (2020), including filtering the 15-min event data to hourly, and applying range, missing values, repeating values, anomaly, spike, and rate of change filters. Lastly, we filtered the dataset to a time period for which all stations had existing and paired surface and bottom temperature data. The final dataset contained 252,711 datapoints from November 2012 to October 2019 (Nelson & Pien, 2021).

Discrete dataset
We used generalized additive models (GAM) to analyze general spatiotemporal patterns of temperature stratification in the discrete dataset. We calculated the bottom-surface temperature difference (bottom temperature minus surface temperature) and used it as a response variable. As such, negative values correspond to cooler temperatures at the bottom and positive values correspond to warmer temperatures at the bottom. Covariates used include geographical coordinates (x, y) and day of year (1 for January 1 st and 365 or 366 for December 31 st dependent on whether it is a leap year). Because temperature stratification can be affected by changes in air temperature at the daily or interannual scale, we also included a "temperature anomaly" covariate that represents the variability of surface temperature relative to the expected daily temperature based on day of year and location. To calculate the temperature anomaly value, we first constructed a GAM to predict surface temperature based on day of year and geographical coordinates at a coarse level: where S is the surface temperature (°C), te is tensor product smooth modeling the interaction between variables, x and y are the geographical coordinates.
Smooth terms were thin-plate regression spline and cyclic cubic regression spline for spatial coordinates and day of year terms, respectively. The GAM parameter k, the basis dimension of each smoother, was set to 5 for day of year and 10 for the spatial coordinates. These basis dimensions were set fairly low for the model because it was used mainly to remove S ∼ te (x, y), (day of year) collinearity between temperature, space, and time. To calculate the temperature anomaly value, we subtracted the observed temperature in our dataset by the predicted temperature from the temperature anomaly model. As such, the final temperature anomaly value provides us with the deviation of the surface temperature measurements from the general expectation based on their location and season within the Estuary. Additional variables we expect to influence stratification such as tide and wind were not included because they were not available at the proper spatiotemporal scale.
We set up four candidate models with smoother terms based on hypotheses (see models 1-4 in Table 1). A conventional thin-plate regression spline GAM with geographical coordinates as predictors would extrapolate the model across space without regards to geographical boundaries (e.g. prediction of water temperature over land). To minimize potential bleedover of information between distinct water bodies that are within a short distance of each other (e.g., eastern end of Suisun Marsh and western end of Cache and Lindsey Sloughs) (Fig. 2), we applied a fourth model with a soap-film smoother (Wood et al., 2008) using the same covariates as the best fitting model out of the original four ( Table 1). Because of the number of islands and channels that exist in the upper Estuary, and the inability of the soap-film smoother to handle exceedingly complex boundaries, we ran our soap-film smoother model with a simplified boundary that represents just the outline of our study area (i.e., without islands) (Supplementary Information Fig. S3). Given that water temperature in the upper Estuary is primarily driven by air temperature (Wagner et al., 2011;Vroom et al., 2017), a relatively simple boundary should have provided accurate results.
Prior to model selection, we conducted a preliminary analysis to select the best basis dimensions (k) for each covariate. We selected k of 5 for the day of year term based on visual inspection of the temperature anomaly calculation model and the general seasonal pattern of the system. For the spatial coordinates term (x,y), we constructed GAMs with bottom-surface temperature difference as the response variable and the tensor product smooth of the spatial coordinates as predictors with varying k. We constructed models from k = 10 to k = 50 at increments of 5, resulting in a total of nine models. We evaluated changes in adjusted R 2 and Akaike information criterion for limited sample size (AICc) across the different k values to select the best k. We selected k = 20 for the spatial coordinates term based on the incremental change with subsequent increase of k ( Supplementary  Information Fig. S1). We followed a similar procedure to select k for the temperature anomaly term, but the GAMs for this were constructed with the tensor product smooth of both spatial coordinates and temperature anomaly. A total of eight models were constructed with varying k for the temperature anomaly term from k = 3 to k = 10, with k for the spatial coordinates term set at 20. We selected a low k at 3 for the temperature anomaly term based on visual inspection of changes in adjusted R 2 and AICc, and to also ease the interpretation subsequent model output (Supplementary Information Fig. S2).
We assessed relative model fit of the four candidate models using adjusted R 2 and AICc (Table 1). Table 1 Description and justification of the four candidate Generalized Additive Models used to evaluate patterns in thermal stratification using the discrete dataset. Models are ordered by increasing complexity. T is the bottom-surface temperature difference, te is tensor product smooth, x and y are the spatial coordinates, day is day of year, and ta is temperature anomaly. For smooth terms, tp is thin-plate regression spline, cc is cyclic cubic regression spline, and so is soap-film smooth Candidate model number

Model structure
Smooth terms Interpretation and hypothesis 1 T ~ te(x, y) x, y = tp Thermal stratification varies across space, but not by season 2 T ~ te((x, y),(day)) x, y = tp, day = cc Thermal stratification varies by space and season 3 T ~ te((x, y),(day),(ta)) x, y = tp, day = cc, ta = tp Thermal stratification varies by space, season, and surface temperature anomaly 4 T ~ te((x, y),(day),(ta)) x, y = so, day = cc, ta = tp Thermal stratification varies by space, season, and surface temperature anomaly. Water boundaries are also important in understanding stratification patterns To further evaluate the performance of the models, we conducted tenfold cross validations with randomly assigned folds to the dataset (but same across models). We assessed the results of this tenfold cross-validation process by calculating the root mean square error (RMSE) and Pearson's correlation coefficient (r) between the out-of-sample predictions and observations. Global RMSE across folds was calculated by dividing the sum of squares of RMSE from the various folds by the number of folds (10), and then taking the square root of this value. We also calculated the proportion of predictions from the tenfold cross-validation that had the same sign as the observed value (e.g., whether the model predicted a negative value when the observation has a negative value and vice versa). For this calculation, actual observations with identical bottom and surface temperatures (zero difference) were removed. To assess the degree of residual spatiotemporal autocorrelation, we calculated the spatiotemporal variogram on the best model's residuals (Pebesma, 2004).
To visualize our results, we first generated a 150 × 150 grid of cells over the spatial extent of our dataset. We then removed any points outside of the Estuary's water boundaries, and subsequently generated predictions from our best fitting model based on the specified covariate values. All analyses for the discrete dataset were conducted in R version 3.8.2 (R Core Team, 2021) using the "mgcv" package for GAM construction (Wood, 2011), "cvTools" for k-fold cross-validation process (Alfons, 2012), "gstat" for spatiotemporal variogram (Gräler et al., 2016), and "ggplot2", "ggpubr", and "sf" package for plotting and visualization (Wickham, 2016;Pebesma, 2018;Kassambara, 2020).

Continuous dataset
Our discrete dataset was collected only during daytime and each site was rarely sampled more than a few times in a given month. As such, we used the continuously collected data in our study to better understand the changes in temperature stratification throughout a diel cycle and across months. We calculated bottom-surface temperature difference in a similar manner as the discrete dataset (bottom temperature minus surface temperature for each time point). To evaluate how bottom-surface temperature differences varied across a diel cycle, we organized time of day into 3 categories based on hour (Pacific Standard Time): Early (00:00-07:59), Middle (08:00-15:59), and Late (16:00-23:59). To assess changes across months we plotted bottomsurface temperature differences by day of the year with gam smooths constructed in ggplot in R with the formula: where y represents the bottom-surface temperature differences, k represents the GAM basis dimension, and bs = "cc" defines a cyclical smooth type. We observed that substantial differences occurred only at two of our stations: RRI and MRZ (Fig. 3). For this reason, we focused on RRI and MRZ stations for subsequent analysis.
For models of RRI and MRZ, we first calculated a surface temperature anomaly, similar to that described for the discrete dataset, to provide context about how stratification varied during cooler or warmer days. For each station, a GAM examining the interactive effect of day of year and hour on surface temperature was run. The basis dimension of each smoother, k, was set to a low value of 5 for both day of year and hour, and we used a cyclic cubic regression spline for both hour and day of year due to the cyclic nature of day of year and hour: The anomaly was then extracted from the response residuals of the model. Next, we ran separate models for RRI and MRZ examining the interactive effect of hour, day of year, and surface temperature anomaly on bottom-surface temperature difference (T) using a tensor product smooth model. We recognized that there was temporal autocorrelation not accounted for in these models; however, we used the models primarily to visualize observed patterns. Basis dimension (k) values were selected based on our expectations and uses of this GAM model for visualization. We chose lower basis dimensions for hour and anomaly (7 and 6, respectively) because we expected simple curves and were only interested in the broader relationships. We chose a k-value of 14 for day of year to roughly reflect the monthly time-scale we were most interested in. We used a cyclic cubic regression spline for hour and day of year, and a thin-plate regression spline for temperature anomaly (ta): y ∼ s day of year, k = 5, bs = "cc �� S ∼ te(day of year, hour) Model fit was assessed by visualizing residual patterns. We visualized model results by predicting model results for all combinations of hour (1-24), day of year (1-365 or 366), and surface temperature anomalies of − 1.5 °C, 0 °C, and + 1.5 °C, and plotting rasters of the results for each anomaly. Analysis for the continuous dataset was done in R version 4.0.1 (R Core Team, 2021). GAMs were run using the "bam" function in the "mgcv" package (Wood, 2011) and visualizations were displayed using the "ggplot2," "ggpubr", "grid" and "gridExtra" packages (Wickham, 2016;Auguie, 2017;Pebesma, 2018).

Study species and temperature thresholds
To put our results in the context of species' habitat, we evaluated our data and model output using key temperature thresholds based on the existing literature for Delta Smelt and Chinook Salmon, two fish species of high conservation concern in the Estuary. Delta Smelt is an endemic, annual forage fish species that spends the entirety of its life cycle in the upper Estuary (i.e., estuary-dependent) (Fig. 1). Delta Smelt T ∼ te(day of year, hour, ta) spawn in the cold and wet winter-spring months in the freshwater portion of the Estuary, followed by a downstream dispersal to low salinity brackish habitat where they rear in the dry summer-fall months (Sommer et al., 2011). High water temperatures in summertime can negatively impact Delta Smelt survival (Mac Nally et al., 2010;Polansky et al., 2020). Catch of Delta Smelt generally declines starting at 20 °C and they are rarely seen when water temperature reaches > 25 °C (Nobriga et al., 2008). Temperature of 25 °C and above is generally believed to cause high Delta Smelt mortality (Swanson et al., 2000;Brown et al., 2013Brown et al., , 2016, and therefore, we chose 25 °C as our temperature cutoff for the species. Chinook Salmon is an anadromous fish species of high commercial and recreational interest that is native to the Pacific coast of North America. Adult Chinook Salmon spawn in Sacramento and San Joaquin river tributaries over a hundred of kilometers away from the Delta. Juvenile Chinook Salmon mostly rear for a few months in streams or the Estuary before making their way into the ocean (Fig. 1). Chinook Salmon display a diversity of life history types. Winter-run and spring-run, the two Chinook Salmon runs federally listed under the Endangered Species Act, spawn in rivers during the dry months between April and  (Moyle, 2002). Most juvenile Chinook Salmon in the Estuary face rising water temperature as they outmigrate, which can lead to increased predation risk. Growth for juvenile Chinook Salmon in the laboratory appear to decline above 20 °C (Marine & Cech, 2004;Myrick & Cech, 2004) and mortality in the field can be near 100% above this temperature due to predation (Nobriga et al., 2021) so we chose 20 °C as our temperature metric for this species.
With the final models we constructed from our discrete dataset, we created surface and bottom temperature predictions that were converted to different suitability categories depending on the species (i.e., suitable temperatures at both surface and bottom, suitable only in the bottom of the water column, unsuitable overall). For Delta Smelt, we plotted temperatures during a typical condition (temperature anomaly = 0) for July 15th, as this would generally be the hottest time of the year. For Chinook Salmon, we plotted temperatures during a typical condition (temperature anomaly = 0) for June 15th because this is the tail-end of the juvenile Chinook Salmon outmigration season and when they would experience the warmest conditions. We also constructed suitability plots with the continuous dataset for the two stations with substantial amount of stratification: MRZ and RRI. Continuous data allow us to better understand the true daily and seasonal pattern of temperature stratification that these at-risk fish species experience. We selected the year 2015 for the plots as this was the warmest average year in the dataset, representative of future conditions, with a relatively complete dataset for both MRZ and RRI stations. We calculated the proportion of each day that was below a threshold of 25 °C for Delta Smelt and below 20 °C for Chinook Salmon. We also subtracted the proportion of suitable hours at the surface by the proportion of suitable hours at the bottom to show when the bottom might provide refuge from the surface at each threshold temperature. While we acknowledge that these thresholds are somewhat arbitrary and that in truth a continuum of suitable temperatures exist for each species, a broad illustration of the thermal landscapes for each species during their warmest months can assist management agencies with identifying areas for potential recovery actions (e.g., hatchery fish release, restoration, etc.).

Discrete dataset
In the discrete temperature dataset, surface temperatures ranged from 5.6 to 29.6 °C with a mean of 18.1 °C and median of 19.2 °C, while bottom temperatures ranged from 5.6 to 29.9 °C with a mean of 17.9 °C and median of 19 °C. Temperature difference (bottom temperature minus surface temperature) in the discrete dataset ranged from -7.4 to 6.5 °C with a mean of − 0.2 °C (standard deviation of 0.6) and median of − 0.1 °C. About 15% of the data had zero difference between bottom and surface temperatures. The surface temperature GAM constructed to calculate the temperature anomaly value followed our general expectations, with surface temperatures being highest in July and August and lowest in January ( Supplementary Information Figs. S4-S6). Although there was some variability in the spatial pattern across seasons, surface temperatures were generally higher towards the Sacramento River Ship Channel to the north and the southern portion of the Delta (Figs. 2,  4).
Of the four models we constructed with the discrete dataset (Table 1), model 4 with the soap-film smoother had the best fit overall with the best performance across all metrics ( Table 2). The temperature anomaly, which could be due to either diel or interannual changes in air temperature, appears to be an important variable for predicting bottom-surface temperature difference based on the substantial jump in R 2 , AICc, RMSE, and r between model 2 and model 3. The addition of the soap-film smoother in model 4 provided only a modest improvement from model 3 based on the same metrics and negligible change in model accuracy based on % of correctly predicted directions of bottom-surface temperature difference (+ or -). The spatiotemporal variogram of model 4 showed some autocorrelation for sites within 1 km of each other sampled on the same day, but no consistent patterns for other time frames or distance within a two-week span ( Supplementary Information Fig S7). There was a slight increase in autocorrelation for data  Table 1) are shown on the right panel. Rows are the dates used to represent the different seasons. Columns within each panel represent three different scenarios based on temperature anomaly values: cool conditions (-1.5 °C from expected surface temperatures), average conditions (0 temperature anomaly, predicted surface temperatures), and warm conditions (+ 1.5 °C from expected surface temperatures). To avoid extrapolation, only regions with data that match the categories (season and temperature anomaly group) are shown for each plot. Note that the color-temperature scale changes among seasons  Fig S8), likely due to the nature of the long-term monitoring programs that produced the dataset where sampling occurs roughly once a month. Both raw data and results from the soap-film model indicated that considerable differences (> 0.5 °C) between surface and bottom temperatures are relatively rare (Figs. 4, 5). Bottom-surface temperature differences tend to be negative (cooler in the bottom of the water column relative to surface), though conditions where bottom temperature is warmer than surface do occur in winter and spring. We observed the highest temperature differences between surface and bottom at the warmest conditions: during the peak of summer and as the temperature anomaly became more positive (Fig. 4, Supplementary Information Figs. S4-S6). Although the fall season saw lower magnitude temperature differences overall relative to summer, cooler temperatures in the bottom are more widespread than other seasons (Fig. 5). Thermal stratification is the most prominent and consistent in the Sacramento River Ship Channel towards the northern end of the Delta, especially towards the northern tip of this dead-end channel. The Sacramento Ship Channel is where we observed the maximum temperature difference between surface and bottom in the discrete dataset at -7.4 °C, and where the soap-film model predicted the highest difference overall with over 3 °C cooler temperatures in the bottom during the summer months. Although infrequent, the soap-film model indicates that warmer bottom temperatures can be observed relatively often at the western end of Suisun Bay during the coldest months and during more negative temperature anomaly conditions.

Continuous dataset
For the entire continuous dataset, surface temperatures ranged from 6.6 to 28.8 °C with a mean of 17.1 °C and median of 17.6 °C, while bottom temperatures ranged from 6.8 to 27.3 °C with a mean of 17.0 °C and a median of 17.5 °C. ANH, MAL and MRZ had similar surface temperature ranges, while RRI had the greatest range and reached the highest temperatures of any station by several degrees (Table 3). Stations ANH and MAL had more moderate bottom-surface differences, while MRZ and RRI exhibited greater negative values and variability in their bottom-surface differences ( Fig. 3; Table 3).
For all stations, the highest temperature differences between surface and bottom usually occurred in the summer (June, July, August). Furthermore, median maximum temperatures occurred either during the 15th or 16th hour (3:00 or 4:00 PM in the afternoon).
Results from the GAMs further indicate that for RRI and MRZ, bottom temperatures can be several degrees cooler than surface temperatures, especially during the summer (July -early September), and in the afternoon, (14:00 -19:00) (Fig. 3, 6). These differences, where bottom is cooler than surface, occur for 35% and 23% of predicted values in RRI and MRZ, respectively, but differences more negative than − 1 °C occur in only 3.4% and 4.7% of predicted values in RRI and MRZ, respectively. During cooler times of year, and the early and late hours of the day, we often observe no bottom-surface temperature difference, or the surface being slightly warmer than the bottom (up to 0.3 °C at RRI and 0.7 °C at MRZ).
While RRI and MRZ exhibit similar patterns by day of year and hour, we observed differences in patterns in relation to the surface temperature anomaly. Station MRZ exhibits similar trends as those observed in the discrete dataset: bottom-surface temperature differences are greater during warmer than expected surface temperatures (greater anomaly), reaching differences of -2.3 °C in the summer at a positive anomaly of + 1.5 °C (Fig. 6). For RRI, we observe the greatest bottom-surface temperature difference (-2.5 °C) during cooler than expected surface temperatures (anomaly of -1.5 °C), followed by warmer than expected surface temperatures (-1.4 °C difference at an anomaly of + 1.5 °C) (Fig. 6).

Species thermal landscape
Based on our surface temperature model and soapfilm model results from the discrete dataset, the Sacramento River Ship Channel and the southern portion of the Delta are two regions that can pose thermal risk to both juvenile Chinook Salmon and Delta Smelt for the average June and July conditions (Fig. 7). Parts of the Sacramento River (Middle Sacramento River region) towards the northeastern portion of our study region also seem to be relatively warm; however, we have low sample sizes for June and July for this area (N = 9). Water temperatures were generally cooler at the bottom and may provide refuge in portions of the Sacramento Deepwater Ship Channel and Model results for the discrete dataset showing areas with significant bottom temperature difference from surface (p < 0.01) to aid interpretation of results. Red areas indicate warmer temperature in the bottom relative to surface, while blue areas indicate cooler temperature in the bottom relative to surface. Rows are the dates that are used to represent the different seasons. Columns within each panel represent three different scenarios: cool conditions (-1.5 °C from expected surface temperatures), average conditions (0 temperature anomaly, predicted surface temperatures), and warm conditions (+ 1.5 °C from expected surface temperatures). Only regions with data that match the categories (season and temperature anomaly group) are shown for each plot. The predicted mean magnitude of difference can be found in Fig. 4 the upper Sacramento River, representing habitat for Delta Smelt and a key migration corridor for Chinook Salmon, respectively. Meanwhile, the southern edge of the Delta seems to offer little reprieve for both species with high temperatures above our thresholds for both the surface and bottom of the water column during June and July of the average year.
The continuous dataset provided a similar picture, where the bottom of the water column can be cooler and more suitable for both at-risk species at different times of year (Fig. 8). While both the surface and bottom are at temperatures likely unsuitable for both species throughout July-September for MRZ, and May-October for RRI, the bottom may provide periods of refuge during the months bordering these hottest months, such as in June and October for MRZ, and May and November for RRI (Fig. 8). Juvenile Chinook Salmon will likely be downstream of Suisun Bay by July (MRZ), so the species would not be too affected by the warm temperatures in July-September. However, Chinook Salmon may be affected by warmer temperatures in April in RRI, as this overlaps with their outmigration period from the San Joaquin River. As such, Chinook Salmon could

Discussion
Evidence of climate change and its impacts on biodiversity continues to accumulate. However, these impacts will not be uniform across space and time, and the spatial and temporal variability in temperature stratification within estuaries may provide an important thermal refuge to key fish species. Our study shows that the upper Estuary, while mostly wellmixed, demonstrates episodic thermal stratification at certain locations, seasons, and times of day that may provide thermal refugia for at-risk species. Cooler water towards the bottom can be observed more often when overall temperatures are the highest: summer and fall months, warmer than average years, and the warmest time of the day. Our study also provides some insights into thermal stratification patterns that could be occurring in other estuarine systems with similar climate and geomorphology: during winter when freshwater outflow is high, warmer waters towards the bottom could occur (possibly related to salinity stratification; see Vroom et al., 2017), while deep terminal channels may result in a fairly consistent thermally stratified water body throughout the year. Moreover, we took a statistical and empirical approach in this study that not only allows for a better understanding of the variability of thermal stratification, but also provides support and validation for mechanistic physical models.

Water temperature dynamics
Although it is generally well-understood that overall water temperatures are rising due to climate change, the thermal landscape may be complex . Results from our spatially intensive discrete dataset indicate that thermal stratification varies considerably across space and time. When thermal stratification occurs, it often results in cooler temperatures in the bottom layer of the water column relative to the surface (Figs. 5, 6), as would be predicted from studies of other estuaries (Strange, 2013;Weinke & Biddanda, 2019). Warmer water is less dense, so it rises in the water column. In addition, surface water is exposed to more sunlight and warm air temperatures compared to deeper water (Simpson et al., 1990). The cooler bottom layer of water is the most consistent and widespread in the warm summer-fall months when freshwater input to the system is typically at its lowest. In highly stratified estuaries with a salt wedge, freshwater flow can intensify stratification by forming a buoyant layer on top of higher-density and colder salt water (Simpson et al., 1990). However, in shallower, fresher estuaries such as the Delta, freshwater flow may increase mixing. This pattern has also been seen in the Chesapeake Bay, where high freshwater flows can cause mixing and breakdown of stratification (Xu et al., 2012).
The bottom layer of the water column also becomes cooler relative to the surface as overall surface conditions become warmer (as proxied by the positive temperature anomaly value), possibly because thermoclines tend to form more frequently during periods of low wind, which often coincide with warmer air temperature. Surface water temperatures are also generally warmer in spring and summer during droughts, which can potentially increase thermal stratification . Conversely, warmer conditions towards the bottom of the water column can occur during cool conditions. The Estuary typically experiences higher freshwater outflow in winter-spring relative to summer-fall, and longer exposure of the freshwater surface layer to cold air during wintertime may explain this phenomenon. Observations of warmer water in the bottom relative to surface are most common at the most downstream location of our study area during wintertime of colder years (Fig. 5), likely a byproduct of salinity stratification that can take place in this more marineinfluenced region (Vroom et al., 2017).
In concordance with previous studies (Kimmerer, 2004;Vroom et al., 2017), the magnitude of temperature difference from thermal stratification in the Estuary is generally low (within 0.5 °C). The sole exception to this is the Sacramento River Ship Channel, a man-made channel built in 1963 to allow large vessel access to the port of West Sacramento, located towards the northern end of the Delta. The Sacramento River Ship Channel is a deep (~ 11 m) and long (~ 40 km) dead-end slough with high residence time and low water exchange with the rest of the system (Gross et al., 2019). However, the Sacramento River Ship Channel is also tidally influenced towards its downstream portion, resulting in unique hydrodynamics relative to the rest of the Delta and a strong thermocline throughout most of the year (roughly between 0.6 and 3.5 °C cooler towards the bottom, with a maximum of 7.4 °C based on our dataset). It is important to note that our discrete dataset can only demonstrate the broad spatiotemporal patterns and cannot capture the short-term variability in thermal stratification due to tides and winds. Thermal stratification in the Sacramento River Ship Channel is more prevalent during ebb tides when warm water from the northern end moves south at the surface layer level (RMA, 2021). Meanwhile, wind and freighter traffic can break down thermal stratification due to vertical mixing at the Sacramento River Ship Channel (Lenoch et al., 2021). Changes to wind patterns due to climate change could impact where and when we see thermal stratification in the future. It is possible that the declining wind speed observed in the estuary (Bever et al., 2018) may have further amplified the thermocline seen at the Sacramento River Ship Channel. Wind has been well-documented in breaking down stratification in other estuaries, including the Chesapeake Bay (Xu et al., 2012), Great Lakes (Weinke & Biddanda, 2019), and the York estuary (Scully et al., 2005).
The temporally intensive, but spatially limited continuous water temperature data provides an important counterpoint to our spatially intensive discrete data set. We found the potential for greater average temperature differences at two of the continuous sites than the discrete data would suggest, possibly because it captures more data from the warmest part of the day (~ 14:00-16:00). These stations also gave us a window into nighttime temperature patterns, which are not captured during discrete water quality sampling runs. Higher stratification during warmer times of day is not surprising, but the difference between the four continuous stations was unexpected, and these differences provide insight into the range of potential thermal landscapes. Station MRZ was furthest to the west, more subject to ocean influence, and had the highest range of difference between surface and bottom temperatures. There were even periods of early mornings during the winter when bottom temperatures were warmer than surface temperatures, potentially due to an interaction between thermal stratification and salinity stratification . RRI was furthest to the east, but in a much smaller channel. High stratification was observed at this site as well, particularly in the summer in the afternoon. The moderate amount of thermal stratification observed at RRI may due it being a deep channel with low water velocity for large parts of the year (Vroom et al., 2017). ANH and MAL were midway between these two end points, and these two stations had very little thermal stratification. We do not have sufficient sample size to speculate on which conditions are controlling these patterns, but the contrast between these four stations indicates more paired surface/bottom continuous temperature sensors can help better understand the thermal landscape of the estuary.
Another interesting observation from the continuous temperature sensors was that while bottomsurface temperature difference increased with temperature anomaly at MRZ, as expected from surface warming, we did not observe this pattern at RRI. The increased stratification during cooler years may be related to inflow and water year type, as the coolest years observed at this station also happen to be the wettest years (2017,2019), and similarly, the warmest years also happen to be some of the driest years in our study period (2014-2016) (Supplementary Information Figs. S12, S13). Meanwhile, the decreased stratification during warmer periods may be due to the interaction between inflow and the exceptional longitudinal dispersion observed in the San Joaquin River where RRI station is located (Monismith et al., 2009). This phenomenon may also be related to the Stockton Deep Water Ship Demonstration Dissolved Oxygen Aeration Facility (Aeration Facility) pumping that occurs at RRI during warmer years and during times of lower levels of dissolved oxygen, which may increase vertical mixing. The intake for the Aeration Facility is at the same location as the RRI water quality station and the diffuser location is approximately 0.2 miles downstream of RRI (ICF International, 2010). Pumps or "bubblers" are commonly used to mitigate local occurrences of low dissolved oxygen levels and have been found to disrupt stratification of temperature and salinity in other systems (e.g., Hamilton et al., 2001).

Implications for species of concern
Identifying thermal refugia for key species under a warming climate is a critical issue for conservation biology (Keppel et al., 2012). In the Estuary, two declining fish species play a substantial role in water management and are of high interest: Chinook Salmon and Delta Smelt. Chinook Salmon is an ecologically and economically important species throughout the Pacific coast of North America, with two federally listed runs under the Endangered Species Act within the Estuary (NMFS, 2019). Meanwhile, the Delta Smelt is an annual, forage pelagic fish endemic to the Estuary that acts as a sentinel species and is nearly extinct in the wild (Hobbs et al., 2017). Both fish species are considered sensitive to warm water temperatures (> 20 °C) often seen in the Estuary (Myrick & Cech, 2004;Moyle et al., 2016).
Based on our results, the northern and southern edges of the Delta already pose thermal risk in late-spring and summer for both Chinook Salmon and Delta Smelt based on the average conditions from the past decade (Fig. 7). However, the impacts of these high temperature conditions and stratification differ considerably between the two species given their distinct life histories.
Chinook Salmon are highly impacted by temperature stress at their juvenile life stage, as they migrate from their spawning grounds in the upper watershed through the Delta on their way to the ocean (Hance et al., 2021). This stage in their migration is generally characterized by high mortality, even during good conditions (Buchanan et al., 2018). If they reach the Delta when temperatures are too warm, they may experience increased predation, physiological stress (Marine & Cech, 2004), low dissolved oxygen (Jassby & Van Nieuwenhuyse, 2005), synergistic impacts of toxins (Dietrich et al., 2014), and other barriers to successful migration. Salmon have been shown to use differences in water temperature to optimize their growth and feeding in other systems (Tiffan et al., 2006;. Therefore, juvenile Chinook Salmon may make use of thermal stratification as they migrate from the relatively warm Sacramento River (Fig. 7a) and San Joaquin River (Fig. 8a, see RRI station on the San Joaquin River) towards the cooler downstream parts of the Estuary, whether it is based on season, location, or time of day (Figs. 4,6). Juvenile Chinook Salmon are almost never found in the frequently stratified Sacramento River Ship Channel ; however, they can often be routed towards the large water export facilities at the southwest end of the Delta where water can be unsuitably warm in June (Kimmerer, 2008). The unsuitably warm waters around the Delta water export facilities suggest that during warmer parts of the years, there may be high mortality rates for juvenile Chinook Salmon prior to their entrainment into the facilities (Jahn & Kier, 2020). We found the most consistent thermal stratification in the warmer months of June-September (Supplementary Information Figs. S9-S11) when few salmon are migrating through the region. Additionally, because the ability for juvenile Chinook Salmon to migrate vertically varies by size (Brown et al., 2009) and therefore dependent the environmental factors they experienced prior to their estuarine entry, our interpretation of the impacts of thermal stratification is rather complicated. Nevertheless, increased temperature due to climate change may extend the period of thermal stratification earlier into the spring or later into the fall based on our results, possibly making cooler depths more important for early or late migrants. Further studies can be done to examine the depth preference of juvenile Chinook Salmon in relation to temperature, as well as extending our analysis to the adult life stage that can also experience thermal stress during migration (Marston et al., 2012).
In contrast to Chinook Salmon, the endangered Delta Smelt reside within the Estuary year-round. Delta Smelt spawn in freshwater and largely migrate towards the low salinity zone around summertime, but parts of their population can reside in freshwater year-round in the northern portion of the Delta that includes the Sacramento River Ship Channel . Semi-anadromous Delta Smelt that migrate into the low salinity zone can experience relatively cool waters in the Suisun Bay in the summer (Fig. 4); however, freshwater resident Delta Smelt have to contend with high and occasionally lethal temperatures in the north Delta (Mahardja et al., 2019;Young et al., 2020). Given the sensitivity of Delta Smelt to high temperature and the considerable amount of temperature stratification in the Sacramento River Ship Channel, it seems likely that Delta Smelt would take advantage of the cooler, deeper water at this location. However, we note that further research is needed to better understand Delta Smelt movement behavior and swimming performance relative to temperature regulation and depth. The relatively high turbidity and food density in the Sacramento River Ship Channel preferred by Delta Smelt help explain why this area hosts the last few remaining Delta Smelt populations over the past few years (USFWS et al., 2020). Brown et al. (2016) predicted water temperature in the Sacramento River Ship Channel under several climate change scenarios and found the region may be too warm for Delta Smelt by mid-century, however Brown et al.'s study was based on surface water temperature. If refugia are available in deeper areas of the channel as our results suggest and Delta Smelt is indeed capable of conducting vertical migration, then the Sacramento River Ship Channel may remain a refuge into the future and this suggests that recovery efforts, such as supplementation of hatchery Delta Smelt into the wild, should not exclude the Sacramento River Ship Channel as a release site.
An important caveat to our results is that the thermal limits of both species are likely to be much more dynamic than what we presented here. For example, high temperatures can negatively impact Delta Smelt at a population level in the summer and fall (Mac Nally et al., 2010;Polansky et al., 2020), their spawning timing , suband whole-organism physiology (Komoroske et al., 2015;Jeffries et al., 2016), and behavior (Davis et al., 2019b). The 25 °C threshold we presented for Delta Smelt is likely too conservative even for the sub-adult life stage, given that at above 20 °C, juvenile Delta Smelt behavior starts to change as they exhibit sublethal stress (Davis et al., 2019b) and their growth rates start to decline steeply (Lewis et al., 2021). We also have concentrated on surface-depth variation in the context of providing thermal refugia solely for fishes. However, temperature stratification may have other implications yet to be investigated in this system. Many lakes and estuaries have found interactions between zooplankton vertical migration and the depth of the thermocline, impacting peak densities of both phytoplankton and zooplankton (Berger et al., 2010;Carstensen et al., 2015;Leach et al., 2018). Stratification may also impact rates of nutrient transport and transformation (Sharp et al., 1986;Testa & Kemp, 2008), as well as dissolved oxygen levels (Murphy et al., 2011).
Our results suggest that the occasional cooler deeper water resulting from thermal stratification are indeed available for fishes, but questions remain about the ability of these species to take advantage of them as temperature refuge and trade-offs associated with the bioenergetic cost. Further field experiment and simulation studies should be used to evaluate the swimming capability and behavior of Chinook Salmon and Delta Smelt, as well as the energetic costs that may be associated with vertical movement. Moreover, we focused mainly on large-scale stratification patterns, only used data collected at two depths (bottom and surface), and did not consider islands that are present in the system. Management actions revolving around commonly stratified regions would benefit from having mechanistic physics-based models that can display the variability of stratification at a finer time-scale (daily or hourly). Stratification can form and break down multiple times within a day due to factors not considered in this study (e.g., tidal cycle, wind), and more mechanistic models would be necessary to anticipate future actions that can alter stratification patterns in the long-term (e.g., dredging, installation of aeration facilities).

Conclusions
Our investigation confirmed previous findings that the Estuary is a well-mixed water body lacking widespread thermal stratification. However, we found that thermal stratification associated with cooler waters towards the bottom can become stronger, more consistent, and spatially extensive when overall conditions are the warmest (i.e., summer-fall period, warmer than average years, warmest time of the day). We also found that deep, terminal channels and other areas with low water velocity show more consistent thermal stratification. We expect that thermal stratification will become more common as temperatures continue to rise, the low flow period shifts to earlier in the year, and droughts become more frequent due to climate change. Our study serves as an important step in describing the thermal landscape of two key Estuary fish species and sheds light into the areas and time at which stratification can potentially offer some respite from temperature stress. We identified one location where frequent thermal stratification overlap with an endangered fish species' hotspot (i.e., Sacramento River Ship Channel), which suggests that such a region can be targeted for future recovery actions such as hatchery fish supplementation or habitat restoration. Warming associated with climate change is an intractable global issue that cannot be easily resolved at a regional level. Thermal stratification is one potentially cost-effective path that can help local management by highlighting areas for recovery actions and/or providing a blueprint to replicate stratification effects. However, many uncertainties remain, and we recommend further studies on the extent to which fish species can take advantage of vertical temperature difference prior to undertaking any major investment for a specific management action.
Van Nieuwenhuyse, and Nicole Kwan for their insightful comments on earlier version of this manuscript. We also thank the anonymous reviewers for their insightful reviews of the manuscript. The findings and conclusions of this study are those of the authors and do not necessarily represent the views of their respective agencies.
Author contributions All authors contributed to the study conception and design. Data preparation and analysis were performed by BM, SMB, CP, and MN. All authors contributed to writing, reviewing, and editing the manuscript. All authors have read and approved the final manuscript.
Funding No specific funds, grants, or other support were received during the preparation of this manuscript.

Data availability
The discrete temperature dataset analyzed during the current study is available in the Environmental Data Initiative repository: Bashevkin (2021). The continuous temperature dataset analyzed during the current study is available in the Environmental Data Initiative repository: Nelson & Pien (2021).

Conflict of interest
The authors have no relevant financial or non-financial interests to disclose.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.