Warming of the Willamette River, 1850–present: the effects of climate change and river system alterations

Abstract. Using archival research methods, we recovered and combined data from multiple sources to produce a unique, 140-year record of daily water
temperature (Tw) in the lower Willamette River, Oregon (1881–1890, 1941–present). Additional daily weather and river flow records from the 1850s onwards are used to develop and validate a statistical regression model of Tw for 1850–2020. The model simulates the time-lagged response of Tw to air temperature and river flow and is calibrated for three distinct time periods: the late 19th, mid-20th, and early 21st centuries. Results show that Tw has trended upwards at 1.1 ∘C per century since the mid-19th century, with the largest shift in January and February (1.3 ∘C per century) and the smallest in May and June (∼ 0.8 ∘C per century). The duration that the river exceeds the ecologically important threshold of 20 ∘C has increased by about 20 d since the 1800s, to about 60 d yr−1. Moreover, cold-water days below 2 ∘C have virtually disappeared, and the river no longer freezes. Since 1900, changes are primarily correlated with increases
in air temperature (Tw increase of 0.81 ± 0.25 ∘C) but also occur due to alterations in the river system such as depth increases from reservoirs (0.34 ± 0.12 ∘C). Managed release of water affects Tw seasonally, with an average reduction of up to 0.56 ∘C estimated for September. River system changes have decreased variability (σ) in daily minimum Tw by 0.44 ∘C, increased thermal memory, reduced interannual variability, and reduced the response to short-term meteorological forcing (e.g., heat waves). These changes fundamentally alter the response of Tw to climate change, posing additional stressors on fauna.


Willamette River Basin, the data we found and use, and the processes that may influence water temperature.

Historical Changes in the Willamette River Basin
Humans have influenced the Willamette River for millennia, and thus the mid-19 th century reflects a watershed that was already influenced by people. For thousands of years the Willamette Basin was inhabited by Native Americans, who influenced the watershed in many ways, including through controlled burns and small-scale fish dams (Boyd, 1999, Taylor, 1999. European and American settlement began in the early 1800s; Portland, founded in 1843, became the largest city in Oregon by 1860 (US Census, 1866). Thus, the beginning our analysis (1850) approximately coincides with the beginning of increasing American influence and settlement (Oregon became a US state in 1858). A number of references detail what is known about the mid-19 th century condition. For example, in the 1850s, approximately 97,500ha of the Willamette Valley was mapped by the Government Land Survey Office as riparian and wetland forest, and was dominated by tree species such as Quercus garryana (Oregon white oak), Fraxinus latifolia (Oregon ash), Acer macrophyllum (bigleaf maple), Alnus rubra (red alder), and Populus trichocarpa (black cottonwood) (Christy and Alverson 2011). Other references are described in the main text.

Notes on Data
We found many previously unused and forgotten records for this study. From 1881-1890, the US Signal Service (USSS) measured top-and bottom Tw at Portland at 11:00 (local time) every day. The successor to the USSS, the US Weather Bureau (USWB) measured Tw from 1941-1961 between 6:30 am and 7:30 am daily (local standard time). We digitized and quality assured the previously unanalyzed USSS and USWB records, which were obtained from the National Centers for Environmental Information (NCEI). A spot-check of US Army Corps of Engineers records from Willamette Rkm 10.5 from 1941-42 (Moore, 1968) showed a general consistency with USWB measurements, to within 1 o C. Measurements of Tw are available from the US Geological Survey (USGS) since 1961, with ~26 station years available in the Portland metropolitan area since 1971 (Table 1). These federal records are supplemented by additional state and local records. Intermittent Grab-sample measurements of Tw are available from the State of Oregon Department of Water Quality, particularly during summer (1949, 1953-present; obtained from the City of Portland). Nearly continuous daily measurements of Tw at the Willamette Falls fish ladder from 1985-2020 were obtained from the Oregon Department of The following supplement provides more context for the chronology of changes in the Supplement Fish and Wildlife. Finally, a long, continuous record has been made available by the City of Portland at half-hourly increments from 1992-1999 and 1997-2015 at the Saint Johns Bridge and the St Johns Railroad Bridge, respectively (see also Annear et al., 2003).
Within Portland, air temperature may be influenced by the urban heat Island effect (e.g., Voelkel et al., 2018). However, the city has been relatively urbanized (cleared of forest) since the beginning of the time series, and Ta measurements have primarily occurred by either the Willamette or Columbia River, both reasons that changes in temperature bias caused by infrastructure may be relatively small. Moreover, the distribution of air temperatures around the climatological mean has remained virtually unchanged ( Figure 5). Given the long history of Portland and later the Airport as the primary regional measurement station, and the consistency of trends with the regional average (e.g., Mote et al., 2019), we conclude that the Ta measurements are reasonably representative of regional climate patterns.
Average air temperatures during the 1881-1890 calibration period (during the Signal Service Tw measurements) are only 0.4 o C cooler than the 2000-2015 calibration period (Figure 11d), markedly lower than the 1.3 o C difference between the 30y climatological averages ( Figure 11b). A possible reason is that pre-1888 measurements may not have been properly sheltered (Mote 2003). However, comparison with Tw measurements (compare Figure 11c with 11e) suggests that air and water temperature patterns during this decade were similar and warmer than previous and subsequent decades. For example, both springtime Ta and Tw measurements in the 1880s were higher than instrumental measurements from the 2000-2015 period. The correspondence between Ta and Tw measurements in the 1880s increases confidence that measurements indicate a real climate signal, possibly caused by decadal fluctuations in climate (e.g., Peterson & Kinkel, 2001), rather than an instrumental artifact.

Notes on the Advection-Diffusion Equation
To develop our statistical model approach, understand its limitations, and motivate its form, we consider here the underlying physical dynamics. So that the discussion is coherent, some repetition with the main text is found.
Heating and cooling of river water is governed by the Advection-Diffusion equation (ADE; e.g., Fischer et al., 1979). When vertical and cross-sectional variations in Tw are neglected, the 1-D ADE for Tw as a function of time t and along-channel coordinate x (positive downstream) reads: where K is a horizontal diffusion coefficient, u is river velocity, H is the sum of heat flux into or out of the system, d is the cross-sectionally averaged depth, and cp is the heat capacity of water, and is approximately constant to within 1% for typical variations in Tw. This simple ADE does not consider groundwater flow, which cools the off-channel alcoves of the Willamette River during summer (Faulkner et al., 2020).
Scaling provides insight into the relative importance of the advection, diffusion, and heating terms, relative to the time rate of change . Over a 12 hour time scale during the day, temperatures in summer are observed to vary by ~0.5 o C, yielding ( )~10 -5 o C/s. Over a month, larger changes of order 5 o C are observed, yielding ( ) ℎ~2 ×10 -6 o C/s. The time rate of change for daily and monthly time scales must be balanced by the terms on the right hand side of Equation (1). An evaluation of measurements suggests that:  The diffusive term is negligible. Over most of the year, the monthly average of daily is << 10 -5 o C/m, except from July-September when a monthly-averaged increase of 1-2 o C per 100km is observed (Figure 2b). Using 100km as a typical length scaleand K ~1000 m 2 /s for the diffusive term, the ( ) term is generally < 10 -7 o C/s, much less than .
 The nonlinear advective term is likely influential during summer, due to a positive ( Figure 2b). During other seasons, river discharge can either cool or warm Portland water because of the presence of both negative and positive ( Figure 2). Therefore, the net influence of the advective term on monthly averaged temperatures is likely small, though it may matter during weather events (such as a rain-on-snow event).  Seasonal variations in discharge (Figure 2a) influence the magnitude of the advective term. During early summertime (June) conditions, Lee (1995) measured velocities of ~0.8 m/s in the upper Willamette; tidally averaged currents are typically 0.05-0.1 m/s during the same period in Portland (USGS Gauge 14211720). Since discharge is smallest during August/September, the decrease in u counteracts the increase in in the advective term . Overall, considering typical magnitudes of u and , we find that the advective term scales as 10 -5 o C/s to 10 -6 o C/s during the summer, depending on location.
 Based on the considerations above, the heating term is usually the leading order term that drives the time rate of Tw, as also found, for example, by Wagner et al., (2011).
When advection and diffusion are unimportant, the non-linear heating term ( ) governs the time rate of change of temperature, . The term can be linearized using a number of assumptions, enabling use of a linear regression approach in which Tw is a function of Ta and river discharge Q. The details, described briefly below, reveal some inherent limitations. See Mohseni & Stefan (1999) for a more detailed discussion of linearization assumptions.
First, we make the approximation that the reciprocal of depth, 1/d, is a function of Q: where a1 and a2 are constants. The negative sign reflects the observation that 1/d decreases (depth increases) as discharge Q increases.
Further, the heat flux term is a function of at least 5 different terms (e.g., Fischer et al., 1979): The sensible heat flux is proportional to the difference between air temperature and Tw (both measured in Celsius): where k1 is a constant that depends on air density and several empirical coefficients, and w is the wind speed at 10m. The energy loss because of evaporative heat flux, He, depends on wind speed, the latent heat of evaporation, and atmospheric conditions, and is generally small in winter but potentially significant in summer (Wagner et al., 2011). The third term, the heat input from radiation from water vapor, is Where , is a constant that depends on cloud cover. When ∆ is small relative to (273.15 + ), such as occurs in the Willamette, Equation 5 is approximately linear with respect to . Similarly, heat loss due to long-wave radiation is modeled as where the power term is approximately linear in for temperature differences < 20 degrees Celsius (see also Mohseni & Stefan, 1999). Finally, the heat input from incoming shortwave radiation, , is a function of sun angle, albedo, and atmospheric effects. Wagner et al. (2011) used the climatologically averaged insolation as a basis function in their Tw model, but most models implicitly assume that is proportional to Ta, (Benyahya et al., 2007).
Combining Equations 3 to 6, and neglecting the evaporation term, we find that H can be linearized as follows: where b1, b2, and b3 are constants.
Combining Equation 7 and Equation 2, the heating term can be approximated by: Where is the approximation error and c1, c2, c3, and c4 are coefficients. Equation 8 shows that even after many simplifications and approximations, there are still nonlinear interactions between terms such as air temperature and river flow (i.e., the term). In practice, it is found or assumed that air temperature is the most important factor in heating, and only the dependence is retained (e.g., Erickson & Stefan, 2000, Webb et al., 2003. Most statistical models implicitly start with this assumption, though some non-linear regression approaches have been applied (see review by Benyahya et al., 2007). For our purposes here, we note that simplifying heating to be a linear function of works best during periods of relatively constant water temperatures and river discharge (see also Mohseni & Stefan, 1999). This is one reason why models calibrated to a specific season such as summer often works better than a model fit to an entire year (see below).
The advection term in Equation 1 can similarly be linearized by assuming that either or Q is constant or slowly varying, relative to the other. This yields either a regression term in Q or in . Removing nonlinear terms, the following linearized basis function emerges: where , , and are coefficients and the minus sign indicates that river flow reduces water temperature. Using the approximation ≈ − −1 ∆ , we find that Tw at time step n is equal to the Tw at the previous time step n-1, plus a correction that is a function of Ta and Q: This is an autoregressive (AR1) process. Hence, at time n-1, Tw is a function of the at time n-2, and the at n-2 depends on at n-3. If we develop and then substitute the solutions for −1 , −2 , …. into Equation 10, we find that where and are regression coefficients at some time lag , C is a constant of regression, and the time period j is chosen to be long enough that the coefficients and effectively become negligible and/or statistically insignificant. The coefficients and can be modeled using an exponential filter approach (e.g., Al-Murib et al., 2019); here, as explained below, we estimate the coefficients directly. At a large time lag, the influence of the time-lagged temperature term in Equation 10 becomes negligible and drops out; hence Equation 11 effectively represents Tw as a function of time lagged Ta and river discharge.
The discussion above suggests that linear regression models have a basis in the underlying physical dynamics (see also Mohseni & Stefan, 1999). However, a number of assumptions and approximations must be made to represent the 1D ADE as a linear model. Factors such as wind, evaporation, time or spatial variation in parameters and heating terms, and alterations in depth are only approximately represented by Tw and Q. Moreover, depending on conditions, different terms (e.g., depth, heat flux, and velocity) may contribute in varying degrees to the overall heat balance. Thus, a linearized representation of average conditions during a particular season may work less well under unusual or extreme conditions.