Which environmental factors control extreme thermal events in rivers? A multi-scale approach (Wallonia, Belgium)

Study area

The study focused on Wallonia, the southern region of Belgium (Europe), which covers an area of 16,901 km². It was based on 92 measurement sites, each of which delimited a watershed (Fig. 1). The watersheds had a total drainage area of 11,309 km². Elevation and topographic slope varied strongly among the studied measurement sites (Fig. 2), as well as land use and shading conditions (Fig. 3). The studied watersheds ranged widely in their size (surface area) and their physical conditions (channel sinuosity, flow regime and water level) (Fig. 2).

Model studied/temporal thermal dynamic characterization

The spatiotemporal ranges of the variables associated with an empirical model (Eq. (1)) developed by Georges et al. (2021) were studied to understand the temporal thermal dynamics over the study area. This model describes dynamics of WT around extreme climatic (temperature) events. The extreme events were identified with generalized additive models (GAMs). Further details on the temporal thermal dynamic model and the identification of extreme events are described in Georges et al. (2021).

(1) $${\displaystyle \frac{d{T}_w}{dt}={\displaystyle \frac{\left[a\left(T_{aM}-T_w\right)+b\left(T_{am}-T_w\right)\right]}{Q}+c}}$$with ${\displaystyle \frac{d{T}_w}{dt}}$ corresponding to the evolution over time of the day maximum WT (T_w). Time here is a period of +7 days and −7 days around an extreme date identified using GAM modeling (i.e., a period of 15 days). T_am = day minimum air temperature, T_aM = day maximum air temperature, Q = streamflow. a, b and c are model variables representing respectively day thermal sensitivity, night thermal sensitivity, and non-convective thermal flux (Georges et al., 2021).

The model includes three composite variables. Variable a relates to the day convective heat flux, which corresponds to a heat transfer from a hot to a cold body depending on the temperature difference. Variable a modulates the temperature difference between the day maximum air and water temperature (i.e., daytime convective flux). Variable b also relates to the convective heat flux at night, when the air temperature typically leads to a lowering of WT. In what follows, variables a and b will be respectively denoted ‘day thermal sensitivity’ and ‘night thermal sensitivity’. Thermal sensitivity is the response of WT to changes in air temperature (Kelleher et al., 2012; Chang & Psaris, 2013). Variable c (denoted ‘non-convective thermal flux’) was added to the equation to incorporate other physical processes influencing WT such as advection not specifically studied because of scant data availability. The model was fitted to extreme events covering the different seasons (37) identified by Georges et al. (2021) at each measurement site (92) during a 6-year period to estimate specific best fit of variables a, b and c. WT data collected by a dense network recording WT at pre-determined intervals of 10 minutes between 2012 and 2018 was used (Fig. 1).

In this study, a total of 1,564 values (17 extreme summer events (Table 1) × 92 measurement sites) of the three variables were studied in relation with environmental variables in order to improve our understanding of the temporal thermal dynamics around extreme events. Figure 4 shows the range of these model variables in the study area between the different summer extreme events.

Dataset: environmental variables studied

A set of environmental variables characterizing land cover, topography (channel slope, elevation), hydromorphology (channel sinuosity, water level, watershed area, baseflow index) and shade conditions, identified in previous studies as the most impactful drivers influencing WT were considered as potentially able to explain the model variables. The water flow and air temperature variables used to build the differential equation (Georges et al., 2021) were therefore not retained as explanatory variables. The explanatory variables taken are spatially and temporally distinct. Spatially, they fall into three groups: area-dependent, station scale-dependent, and stream network-dependent.

Spatial scales of environmental variables

We then investigated the role of environmental variables at different spatial scales (Fig. 5) to determine whether the local environment had a greater influence on the WT thermal dynamic than the upper reaches or the whole catchment. Five area-dependent variables were measured only at the site measurement scale: watershed area, local channel slope, elevation, water level and baseflow index. All the others were extracted from six spatial scales:

• Catchment scales: the entire upstream catchment (CA) and “circular” buffer of 1 km (buff_1) and 2 km (buff_2). The circular buffer corresponds only to the part intersecting the watershed, upstream of the WT site measurement;

• Riparian buffer scales: 50 m wide corridors on each side of the stream along the entire upstream river network (ripa), within the circular buffer of 1 km (buff_1_ripa), and 2 km (buff_2_ripa) upstream river.

The spatial scales were designated based on previous studies (Wu et al., 2020; Chang & Psaris, 2013; Sliva & Williams, 2001; Scott et al., 2002), context of rivers studied, and field investigation.

The different environmental variables were calculated at the different spatial scales using R (R Core Team, 2018) and QGIS (QGIS Development Team, 2020).

Because some variables (riparian shade and water level) are time-dependent as stated above, and different spatial scales are studied, Table 2 summarizes the variables and their potential variation in space and time. Like the model variables, the environmental variables were also calculated for each measurement site and at each extreme event when the variable depended on it.

Statistical analysis

Results were examined in terms of seasonal timescales focusing on summer when the solar radiation is strongest; the vegetation is fully grown (maximum shade) and the ecological impact most critical for the aquatic ecosystem (Langan et al., 2001; Meisner, 1990; Miara et al., 2018). The extreme dates highlighted by Georges et al. (2021) were therefore assigned to a meteorological season (winter: December–February, spring: March–May, summer: June–August, autumn: September–November). For the summer period (Table 1), a single value of each environmental variable was retained for each measurement site (all extreme dates combined). Only shade and water level varied temporally. They were therefore aggregated by the mean to obtain a single value for each measurement site. Finally, 1,656 observations were listed (92 measurement sites × 18 variables (15 environmental variables + 3 model variables)).

Relationships between model variables (day thermal sensitivity, night thermal sensitivity, and non-convective thermal flux) and environmental variables were evaluated for the six spatial scales during summer with stepwise multiple linear regressions (MLR) to ensure non-collinearity. Prior to all statistical analysis, environmental data were standardized to zero mean and unit variance.

To understand the relationship between environmental variables and the model variables, stepwise MLRs were carried out for each combination of model variable and spatial scale (24 models: six spatial scales × three variables). The effectiveness of these 24 models was assessed based on the adjusted R-squared value (Adj. R²). In the stepwise selection, best influential environmental variables based on stepwise AIC values were selected. The significance of each variable selected from the stepwise procedure was assessed by MLR with the associated p-value. The statistical significance was conventionally set at 0.05.

The estimate value of each variable retained by the stepwise MLR was analyzed. Estimate value corresponds to a coefficient quantifying the influence of the variable taking into account the influence of all the others.

Statistical analyses were performed using R version 3.5 (R Core Team, 2018).

A file with the model variables (variables to be explained) and the environmental variables (explanatory variables) is provided for each spatial scale as a supplemental file. The data is standardized as explained in the methodology. The description of the different column (variables) is to be considered with Table 2.

Descriptive analysis

The watersheds studied in the analysis covered a wide range of environmental conditions. The result of descriptive statistical analysis is presented in Fig. 2. Figure 3 shows the different land cover classes of the watersheds studied.

Most stations were located in watersheds ranging between 56 ha and 115,207 ha and at low elevation (<200 m). Water levels measured were low in summer (<0.5 m) but could reach 1.06 m. On the other hand, most of the watersheds had a high BFI, the values being closer to 1 than 0. Watersheds were therefore mainly composed of aquifers with groundwater storage. Sinuosity ranged around 1.15 (moderately meandering), but covered a wide range of conditions from straight to strongly meandering (Horacio, 2014).

Land cover of the 92 watersheds studied showed a broad diversity. Northern watersheds were mainly composed of cropland, while forest and grassland dominated in the south and east.

Influence of environmental variables and spatial scale on the thermal sensitivity variables

To understand the influence of environmental variables and spatial scales on the three model variables, stepwise MLRs were conducted. Results of the stepwise MLR between model and environmental variables are shown in Table 3.

Influence of other variables

Our results show that at most 60% of the temporal thermal dynamic variability was explained by the environmental variables studied (Table 3). This result suggests that environmental variables other than those studied also influenced thermal dynamics around extreme events. For example, thermal sensitivity is also influenced by human activities and thermal pollution (from industries), deforestation, flow modification and global change. Abstraction of water from a river reduces water volume and so thermal capacity. These modifications lead to an increased maximum WT and a decreased minimum WT (Dymond, 1984). Dymond (1984) analyzed the WT changes associated with abstraction of water by modeling the downstream WT from upstream data (upstream WT, net energy flux, discharge) and reusing this model, varying the discharge to quantify the change in WT due to water abstraction.

Prediction of WT under reduced flow, due for example to irrigation, hydroelectric power stations, or water withdrawal, was also analyzed by Bartholow (1991), Morin, Nzakimuena & Sochanski (1994), and Hockey, Owens & Tapper (1982).

Water releases (reservoirs, dams, industrial effluent) have also been shown to influence thermal downstream conditions (Lowney, 2000).

Day thermal sensitivity a

Watershed area (AREA) was a very highly significant variable influencing day thermal sensitivity at all spatial scales (Table 3). Hilderbrand, Kashiwagi & Prochaska (2014) found no relationship between thermal sensitivity and basin size. Nevertheless, the increase in day thermal sensitivity with the size of the watersheds is explained by stream size. Larger watersheds have gentler hillslopes and larger streams, with therefore less forest shading (Chang & Psaris, 2013; Culler et al., 2018). In addition, thermal sensitivity increased with AREA because the time spent in the river was longer, implying a longer equilibration of WT with air temperature (Beaufort et al., 2020). Time spent determines the amount of exposure to radiation (Wehrly, Wiley & Seelbach, 2006).

Shade controls day thermal sensitivity during summer. Studied with the other environmental variables, the correlation between shade and day thermal sensitivity was negative. A correlation analysis showed that only forest land cover was significantly correlated with shade (correlation = 0.41; p = 0). The other variables had no influence on shade.

Most prior research studied the influence of shade on maximum stream WT (Caissie, 2006; Rutherford et al., 2004; Li et al., 1994; Abell & Allan, 2002; Hrachowitz et al., 2010). However, few studies have addressed the influence of riparian shade on thermal sensitivity as here (Chang & Psaris, 2013; Beaufort et al., 2020). Beaufort et al. (2020) also demonstrated that the thermal sensitivity of a stream decreased with the riparian vegetation, which captures solar radiation.

Wawrzyniak et al. (2017) modeled the shade cast by the vegetation on the river surface using LiDAR data to study its thermal effect. Their results show that one tenth of solar radiation was intercepted by riparian vegetation. This corresponds to a temperature decrease of 0.26 ± 0.12 °C (summer 2010) and 0.31 ± 0.18 °C (summer 2011). Riparian shade retains direct solar inputs and heating (Hilderbrand, Kashiwagi & Prochaska, 2014), so the watercourse heats up less quickly. Net radiation can be five times lower in a stream under shading riparian forest than in an unshaded stream (Moore, Spittlehouse & Story, 2005). Our study confirms that thermal sensitivity is controlled by shade mainly during summer because solar heating (solar radiation) is more intense then (Needham & Jones, 1959; Beschta, 1997; Wawrzyniak et al., 2017).

Our first hypothesis is thus supported by our stepwise MLR showing the influence of different environmental variables on day thermal sensitivity a.

Night thermal sensitivity b

The proportion of storage area, aquifers and reservoirs compared to watersheds with low storage represented by the BFI index in Table 3 influenced night thermal sensitivity negatively. Decreased BFI is associated with increased night thermal sensitivity, meaning that watersheds with small storage areas are more sensitive. In other studies, Kelleher et al. (2012), Chang & Psaris (2013), and Beaufort et al. (2020) analyzed BFI as a measure of groundwater contribution, in relation to thermal sensitivity. They also demonstrated that thermal sensitivity increased with decreasing groundwater contribution.

Other studies specifically analyzed the relationship between groundwater inflows and thermal sensitivity (Luce et al., 2014). Erickson & Stefan (2000) demonstrated that groundwater inflows limited the influence of air temperature on WT, lowering thermal sensitivity. In their study on the thermal sensitivity of Pennsylvania streams, Kelleher et al. (2012) showed that narrow annual variation (about 3 °C) was observed in some sites with high groundwater inputs. These observations can be explained by the fact that groundwater inflow supplies water at constant temperature, which exerts a stabilizing influence on stream WT (Erickson & Stefan, 2000). WT change in response to air temperature change is reduced, causing a decrease in thermal sensitivity and an increase in thermal inertia (O’Driscoll & DeWalle, 2004).

In summer, hyporheic exchanges supply groundwater that is generally cooler than mainstream water (Mayer, 2012). Thermal and water exchanges occur in the stream’s internal structure because the aquifer is hydraulically connected to the main channel (Boano et al., 2006). A back-and-forth transfer therefore occurs between the stream and the aquifer (Gibert, Danielopol & Stanford, 1994; Cristea & Janisch, 2007). Other studies have found that hydraulic properties of the aquifer are a driver of thermal sensitivity in streams (Kelleher et al., 2012; Mayer, 2012; Cristea & Janisch, 2007). The total amount of groundwater contained in aquifers depends on rainfall, the intensity of extraction (agriculture, factories, etc.), and on the nature of the subsoil (SPW Agriculture, Ressources naturelles et Environnement, 2020). For example, in a region with a permeable subsoil, such as chalk, the proportion of precipitation reaching the groundwater is very high, as opposed to a schist region (SPW Agriculture, Ressources naturelles et Environnement, 2020). Johnson (2004) found that maximum WT in the upstream bedrock reach was up to 8.6 °C higher than maximum WT in the downstream alluvial reach and 3.4 °C lower for minimum WT (distance between measurement stations: 350 m). However, sedimentation of the channel, for example through forest harvesting, can decrease the influence of hyporheic exchange through stream bed clogging (Moore & Wondzell, 2005).

Although significant control of BFI is observed for night thermal sensitivity, groundwater inflows also influence day thermal sensitivity. The studies cited above refer to day WT (Kelleher et al., 2012; Chang & Psaris, 2013; Beaufort et al., 2020). No significant effect of the BFI on day thermal sensitivity was observed in our study, probably because our watersheds were large enough to allow greater stream interaction with the surrounding environment (especially shade) during the day (Erickson & Stefan, 2000), masking the influence of daytime groundwater inflows.

Our first hypothesis is thus supported by our stepwise MLR showing the influence of different environmental variables on night thermal sensitivity.

Management insights

The analysis of model variables reflecting thermal sensitivity makes it possible to highlight the most influential environmental variables. The variables studied are not likely to be all equally important in the different watersheds and spatial scales. Determining which environmental variables most influence WT and at which spatial scales is necessary to provide management advice for counteracting heating due to climate change.

Shade, watershed area, water level and BFI were identified as the most influential variables. Among these variables, shade and water level (and more generally the river morphology) are the only relevant variables on which a manager can easily act. O’Briain, Shephard & Coghlan (2017) found that shallower sites showed greater WT extremes, together with sites with little riparian tree cover. Our results are consistent with those of O’Briain, Shephard & Coghlan (2017), suggesting that the thermal resilience of rivers can be improved by riparian vegetation (shade) and restored channel morphology.

Typical stream morphology restoration focuses on reinstating geomorphic features (riffle/pool, meanders) and modifying width-to-depth ratio (Malavoi & Bravard, 2011). O’Briain, Shephard & Coghlan (2017) showed that restoration of riparian tree cover is a more sustainable approach because it allows natural recruitment of large wood debris, recognized as a catalyst and a restoration tool impacting channel morphology within a few years (Osei, Harvey & Gurnell, 2015). Riparian tree cover also provides shade, mitigating WT rise. Moreover, acting on riparian cover is commonly recognized as an accessible management tool in terms of application, cost and time for environment managers (O’Briain, Shephard & Coghlan, 2017).

In restored wood-laden rivers (large wood debris), Osei, Harvey & Gurnell (2015) have demonstrated that the channel does not remain wide. Fine sediment accumulates in and around wood jams and flow velocities increase, leading to the mobilization of fine sediment. Pool formation and sediment bar deposition also occur in wood-restored streams (Elosegi et al., 2017; Pilotto et al., 2016).

Besides wood debris production, riparian vegetation produces shade. Shade cast by trees reduces energy inputs to the river and so counteracts thermal extremes. This observation was also made in other studies (Broadmeadow et al., 2011; Hannah et al., 2008; Garner et al., 2015). Future river management needs to preserve (Feld et al., 2018) and include riparian planting (Brown et al., 2010; Garner et al., 2015; Johnson & Wilby, 2015). For example, Feld et al. (2018) review the use of fences to avoid riparian degradation by livestock. Romasco-Kelly (2018) reports in the Walla Walla basin of southeastern Washington (USA) that local organizations set up riparian restoration projects to use shade from trees to diminish thermal pollution.

Merely advocating the plantation and upkeep of riparian vegetation oversimplifies the actions to be implemented by managers. Using vegetation as a tool to reduce heat energy input into the water requires considering different ‘vegetation variables’ influencing the shade provided (Wondzell, Diabat & Haggerty, 2019). Rutherford, Meleason & Davies-Colley (2018) show the influence of the ratio of tree height to stream width on stream shade. The shade effect is more pronounced on narrow streams (<5 m wide) (Whitledge et al., 2006). According to Feld et al. (2018), length, width, and density should also be considered, as should species (LeBlanc, Brown & FitzGibbon, 1997; Lenane, 2012). Collins et al. (2006) made a literature search on riparian functional width, and showed that the average recommended or observed maximum riparian width for the water body cool function was 40 m. The literature review of Sweeney & Newbold (2014) concludes that buffers ≥30 m wide are needed for full protection against thermal change. Chang & Psaris (2013) found that riparian shading was only significantly negatively correlated with thermal sensitivity 1 km upstream relative contributing area scale.

The effect of riparian vegetation on WT can also vary depending on external factors: stream flow, channel morphology, anthropogenic action (grazing, harvest), and hydraulic structures (reservoirs, dams) (Justice et al., 2017; Wondzell, Diabat & Haggerty, 2019). The orientation of the watercourse (included in our shade calculation) will also influence the amount of shade on the river. As demonstrated by Pluhowski (1970), a north-south watercourse receives less sunshine than an east-west. In the Northern Hemisphere, south bank vegetation maximizes shade on the stream (Larson & Larson, 1996). Streams moderated by groundwater inflows are less sensitive to riparian shade (Johnson & Wilby, 2015).

Besides mitigating thermal pollution, acting on riparian vegetation improves the connectivity between rivers and their floodplains (Opperman et al., 2010), and reduces, retains (Booth & Leavitt, 1999), and filters urban runoff (Welsch, 1991).

While the riparian zone influences thermal sensitivity, among other things, specific management features must be considered to make strategic choices (Beaufort et al., 2020). It is not necessary to maintain or plant a riparian buffer around all streams. Managers need to plan their investments (finance and staffing workforce) strategically. Adding vegetation will be ineffective against thermal heating if other factors are more impactful, such as groundwater flows or stream configuration (width, orientation, topographical conditions, etc.) (Johnson & Wilby, 2015). However, trees can be planted as a priority where there is no vegetation and where the BFI is low (no storage).

Spatial scale

Concerning the spatial scale, Sliva & Williams (2001) demonstrate that variables at the catchment scale have slightly greater influence on water quality than those at the buffer scale (100 m). Pratt & Chang (2012) show that the watershed scale has a greater effect on water quality than the buffer scale (100 m riparian) except for WT. In their review on relationships between land use, spatial scale and stream macroinvertebrate communities, Sponseller, Benfield & Valett (2001) found that land cover influenced mean and maximum WT within respectively riparian corridors and riparian subcorridors (200 m).

The studies cited above show the difficulty defining the most influential spatial scale for management. Although choosing the spatial scale is controversial, as observed by Chang & Psaris (2013), our study shows that managers can focus on relatively small geographical areas.

Although no one spatial scale appears clearly more influential than any other in our analyses (Adj. R² in Table 3), the near stream area appears sufficient to impact thermal sensitivity. The environment across the entire riparian or the entire watershed spatial scales does not significantly explain temporal thermal dynamics around extreme events any better than the local scale, meaning that large scale management will not be more efficient. Because the financial and time/labor costs need to be considered in a management project, the results of our study suggest that managers should focus on small spatial scales. However, the too-small number of environmental variables (three: landcover, sinuosity and shade) at different spatial scales cannot validate our second hypothesis.