Introduction

Water temperature has a significant impact on several water quality issues, such as dissolved oxygen concentration, the toxicity of chemical substances, and the rate of self-purification processes (Caissie 2006). However, most often studies on water temperature are related with stream fauna, as this abiotic variable has been considered as one of the main axes of the environmental (ecological) niche (Allan and Castillo 2007), mainly due to its key role in regulating the physiological processes of aquatic organisms (Gillooly et al. 2001). The biological effects of water temperature, discussed particularly in the case of fish, were presented in numerous studies and concern issues primarily related to their survival (Selong et al. 2001), growth rate (Handeland et al. 2008), gamete maturation (Lahnsteiner and Kletzl 2012), and egg development (Kupren et al. 2008). Thus, various properties of the stream thermal regime have a clear direct effect on their life cycles (Malcolm et al. 2008). In that way, temporal water temperature fluctuations were proven to determine the bioenergetics and growth rate of fish, the timing and duration of the thermal events are important in the context of reproduction success and migrations of fish, while maximum temperatures influence their occurrence and survival (Dunham et al. 2003; Benjamin et al. 2016). In consequence, water temperature studies are of great practical importance, as optimization and appropriate management of fisheries requires determining both spatial and temporal water temperature dynamics, particularly in streams inhabited by ecologically valuable fish species (Malcolm et al. 2008; Broadmeadow et al. 2011). Apart from issues associated with fisheries and aquatic ecology water temperature was commonly used in recent years as a climate change indicator (Isaak et al. 2012; Ficklin et al. 2013; Ul Islam et al. 2019), as well as a tracer of groundwater flow patterns (Anderson 2015). Finally, temperature monitoring seems to be highly valuable in the case of rivers and streams flowing in the proximity of sprawling urbanized areas, where changes in thermal characteristics can be used in the evaluation of future anthropopressure.

In recent years, the thermal characteristics of the streams and rivers, mainly due to the development of new measurement devices and techniques, such as digital temperature recorders and airborne thermal imagery, were broadly recognized (Webb et al. 2008; Fullertone et al. 2015; Jones et al. 2017). Studies concentrated mainly on the spatial heterogeneity of water temperature, from the microhabitat (Surfleet and Louen 2018; Dzara et al. 2019) to the catchment scale (Hrachowitz et al. 2010; Bladon et al. 2018). Considerable attention was paid to the impact of vegetation cover and its influence on thermal characteristics (Kristensen et al. 2015; Dugdale et al. 2018), as well as to the relations between landscape properties and selected thermal regime parameters (Chang and Psaris 2013; Watson and Chang 2018). Anthropogenic transformations of the thermal regime due to reservoir operations, weirs, and wastewater inflows were also discussed (Bae et al. 2016; Maheu et al. 2016). However, despite the significant interest in water temperature, studies related to its spatial variability concentrated mostly on foothill and mountain streams and rivers, inhabited by salmonids (Hrachowitz et al. 2010; Imholt et al. 2013a; Dick et al. 2015; Kędra and Wiejaczka 2018; Żelazny et al. 2018). Lowland temperate watercourses have received definitely less attention, despite the presence of such large areas not only in Poland, but also in another countries in Central and East Europe. Moreover, most of the existing studies on stream temperature in lowland landscape, conducted particularly in Poland, were based on spot, short-term measurement campaigns (Rajwa-Kuligiewicz et al. 2015) or government temperature monitoring, which only covered selected gauging stations (Marszelewski and Pius 2015; Ptak 2017; Graf 2018); in consequence, they concerned large rivers with catchment areas counted in thousands of square kilometers. Therefore, to date, small, ecologically valuable lowland watercourses have not received a thorough study, based on reliable, high-resolution measurements conducted with data loggers.

The main objective of this field-base study was to present the overall tendency of spatial variability of water temperature in lowland temperate watercourses. The specific objectives of the study were to: (1) describe similarities and differences in the basic thermal parameters of the watercourses in two spatial scales, namely main rivers and tributaries; (2) identify main catchment and channel properties suitable for characterizing spatial variability of water temperature parameters; (3) establish linear regression between air and water temperature, which could indicate the relative meteorological control of the thermal regime of watercourses.

Materials and methods

Study area

Water temperature was recorded in 20 catchments located in central Poland in the proximity of the Warsaw agglomeration area (Fig. 1). All measurement sites were distributed across the Mazovian Lowland and partly the Southern Podlachia Lowland, shaped by the presence of the Scandinavian ice sheet during the Ice Age. The landscape of the catchments is dominated by post-glacial flat plains formed as the result of erosion and deposition processes, surrounded by vast uplands, built mainly of sandy loam and boulder clays (Nowakowski 1981). In some parts, sands and gravels formed dune terraces, whereas well stratified sandy-gravel and silt deposits fill ice marginal valleys. The elevation of the catchments ranges from approximately 79 m a.s.l. in site R9 on the Utrata River to only 225 m a.s.l. near the springs of the Rządza River. The catchments are located in a transitional zone of warm temperate climate (from marine to continental), and due to the small differences in elevation, as well as the similar dominant direction of air mass inflow, spatial diversity of precipitation and air temperature is relatively small; the sum of long-term annual precipitation reaches approximately 500–550 mm, while the mean annual air temperature ranges from 8 to 9 °C (Somorowska and Łaszewski 2019). The hydrological regime of the watercourses can be considered to be nival, with the highest streamflow rates observed in March and April in result of snowmelt, and the lowest usually occurring during the hot summer, mainly throughout July, August, and September (Somorowska and Łaszewski 2017). Most of the investigated reaches have retained a quasi-natural character; in such reaches, riparian zones are covered with lines of medium-height trees and shrubs, such as white willows (Salix alba), common aspens (Populus tremula), and black alders (Alnus glutinosa). Due to intensive fish stocking, the ichthyofauna of the watercourses is represented by over 20 species, among which the most important are rheophilic cyprinids (Squalius cephalus, Leuciscus leuciscus, Chondrostoma nasus) and salmonids (Salmo trutta m. fario, Thymallus thymallus) (Borzęcka et al. 2012; Łaszewski 2016).

Fig. 1
figure 1

Study area with the location of measurement sites

Field data collection

Water temperature data was collected to encompass the thermal heterogeneity of lowland watercourses in two different spatial scales—main rivers (referred to as the MR scale) and tributaries scale (referred to as the TR scale). Measurement sites were selected on the basis of the spot measurements, conducted with the use of the handheld thermometer in April 2017 during a non-radiative, cloudy day. Such instantaneous measurements were carried out in a total of 30 sites in the length of the four small-sized lowland rivers—the Jeziorka, Rządza, Świder, and Utrata rivers, as well as in 30 first- or second-order streams in the Świder River catchment. The measurements indicated that across the entire MR scale sites water temperature contrasts were very small and differences resulted mainly from human impacts, such as reservoirs and sewage inflows. On the other hand, in the TR scale sites significant thermal heterogeneity was observed despite their spatial proximity. On the background of such spot measurements, as well as field observations and GIS analysis, 25 sites were subjectively selected, as in the case of similar hydrological studies related to water temperature, where purposive sampling scheme was applied (e.g. Brown et al. 2010; Broadmeadow et al. 2011; Imholt et al. 2013a, b). In that way, 10 sites were distributed across the MR; measurements were taken in the upstream, middle, and downstream parts of the Świder and Utrata Rivers, while in the case of the Jeziorka and Rządza Rivers—only in the upstream and downstream parts. Consequently, sites represent catchments with a diverse area (ranging from 35.6 to 854 km2), mean slope, land cover, as well as morphometric parameters of the channel and riparian shading. However, due to the urbanization of the catchments, some of the measurement sites (R3, R8, and R9) deliberately reflect different anthropogenic influences (Somorowska and Łaszewski 2017; Łaszewski 2018). Sites in the TR scale, with the catchment area ranging from 2.9 to 24.3 km2, were selected with a view to maximizing differences between catchment and channel properties, while their localization precluded direct anthropogenic modification of their thermal regime; moreover, an important criterion was whether these streams were permanently flowing during the summer period. Due to several random incidents during the measurement period (vandalism), one site from the MR scale and four sites from the TR scale had to be removed, so only 20 sites were taken into account in the further course of analysis; their detailed catchment and channel properties are presented in Table 1.

Table 1 Detailed catchment and channel properties of all of the measurement sites

Water temperature measurements were carried out from the 1st of May 2017 to the 31st of October 2017, covering the summer half of the Polish water year. The winter half-year was omitted due to the frequent occurrence of ice phenomena, which pose a direct threat to the temperature loggers, as well as the reliability of the measurements; on the other hand, spatial variability of the water temperature during the winter was usually considered as negligible (Dick et al. 2015). Temperature was recorded with a temporal resolution of 30 min by HOBO U22-001 and UA-001-08 data loggers (Onset Computer Corporation, accuracy of 0.2 °C and 0.44 °C, respectively). Each temperature recorder was installed in a perforated PVC pipe attached to a block of concrete or large wooden log and hidden on the streambed. The location of the loggers was chosen carefully to exclude places with cross-sectional water temperature; to this end, the devices were located in riffles at a depth of approximately 40–60 cm in a turbulent, mixing type of flow. Prior to deployment, the recorders were checked in an ice bath, a simple method described by Dunham et al. (2005). In the case of the U22-001 devices, differences of 0.07 °C were obtained, whereas in the case of UA-001-08 devices − 0.18 °C. Additionally, measurements carried out by temperature loggers were inspected in the field by a certified handheld electronic thermometer to confirm whether the water temperature recorded by the logger at a given point is uniform with the temperature measured in the cross section of the channel. Only small differences not exceeding ± 0.1 °C were noted, which makes the results obtained by the loggers comparable.

Calculation of the environmental metrics

Based on previous literature investigations (Hrachowitz et al. 2010; Grabowski et al. 2016; Żelazny et al. 2018) several catchment and channel metrics were chosen to assess their influence on spatial water temperature variability. These features, calculated for the area or watercourse section upstream from the measurement sites, have direct and indirect effect on heat transfer between water and the surrounding environment; however, some of them are considered as surrogate of other, more difficult to obtain parameters (Malcolm et al. 2004). Nearly all catchment and channel metrics were estimated with the use of raster and vector processing tools from Spatial Analyst extension in ArcMap 10.5 GIS software. In that way, catchment area (A) of the measurement site, as well as its distance from the source (DS) were calculated on the basis of vector layers of Polish hydrographic maps. Mean catchment elevation (E) and its mean slope (S) were computed with the use of raster, digital terrain model SRTM (spatial resolution of 30 m). Riparian shade (R) in a 1-kilom-long reach upstream from the measurement sites was assessed with the use of a modified method proposed by Johnson et al. (2014). Water surface covered by the line of riparian trees on the aerial photographs from public Geoportal of Poland was considered to be a shaded section, while the unshaded section was characterized by water surface completely visible from aerial level, as in this case trees did not form dense lines along the watercourses. The contribution of forested areas (F) was calculated on the basis of land use layer from Corine Land Cover 2012; both coniferous and deciduous forests were linked into one class. Width and depth of the channel, used further to compute the width:depth ratio (W:D), were measured during low flow period in August 2017 directly in the measurement site, as well as in transects of 50, 100 and 500 m upstream and were then averaged.

Statistical analysis

Initially, measurement data obtained from data loggers was inspected to find unnaturally high or low temperature values, related to the temporary sensor exposure to air temperature during data downloading (the time of such cases was precisely known). This values were omitted and interpolated with the use of the neighboring measurement values. Such corrected dataset was then used in the multivariate statistical analysis.

To asses similarities and differences of the water temperature in the main river and tributary scale sites, on the basis of verified, 30-min water temperature data, statistical distribution of magnitude and variability thermal parameters in daily time scale—mean, maximum, minimum, and range—was presented on the box and whiskers charts. Furthermore, hierarchical cluster analysis (CA) was used to support determining the similarity of temperature parameters between investigated sites more objectively. For this purpose, the Ward agglomeration method with Euclidean distance as a measure of similarity was applied. In the Ward method estimation of distance between clusters is based on the variance analysis approach (Ptak et al. 2018); clusters are combined to minimize the sum of total square deviations within them. The literature presents many different methods for determining the optimal number of clusters, however, Żelazny (2012) stated that the most important aspect in this process is whether the aggregated clusters can be subjected to a logical interpretation, based on the author's intuition and good knowledge of the research objects. In the current study, optimal number of clusters was indicated by the highest value of the qj index, used previously by Gutry-Korycka et al. (2009) and given by the formula:

$$q_{j} = \frac{{d_{pq}^{i} }}{{d_{pq}^{i - 1} }}$$

where dpq—distance between two neighbor clusters p and q, i—clustering step (i = 2, …, n).

CA for all measurement sites was conducted using mean, maximum, and minimum daily water temperatures, as well daily range values.

To visualize the relationships between thermal parameters and calculated catchment and channel metrics and thus provide environmental context for the spatial variability of water temperature across the investigated sites, principal component analysis (PCA) was applied in the study. The PCA transforms the dataset into a new, linear combination of original variables, called principal components, which allow to describe the direction and strength of the relationships between them (Maheu et al. 2016); to this end, it was broadly used in recent water quality studies (Parinet et al. 2004). Metrics, initially pre-selected for the PCA, included all catchment and channel metrics described above (A, DS, E, S, R, F, and W:D), while thermal parameters included magnitude parameters calculated for the whole investigated period—mean (T_mean), mean maximum (Tmax_mn), and mean minimum (Tmin_mn) water temperature—as well as variability parameters—standard deviation of water temperature (St_dev), maximum daily range (Max_rg), and mean daily range of temperature (Mean_rg). Additionally, there was estimated duration of water temperature (in % of the total time of study period) fell within 4–19 °C (T_opt), which is considered as the optimum temperature range for the growth and feeding of coldwater fish; this criterion was based on the thermal requirements of brown trout, as previously reported previously in literature (Elliott 1981). Before PCA, variables were tested to establish whether they meet the assumptions of factor analysis. The Shaphiro-Wilk goodness of fit test was used to check the normality of distribution of all of the variables; the results indicated that St_dev, A, R, and DS not represent a normal distribution (p < 0.05), so they were transformed with the use of a logarithmic function (Log10), similar to Biem Mori et al. (2015). After such transformation, only in the case of A and DS statistical distribution became satisfactory (p > 0.05). However, due to an unfavorable sample-to-variables ratio and simultaneously small informational value, DS, as well as R and St_dev variables were removed from the analysis. The same was done with Max_rg parameter, which is strongly correlated with the Mean_rg and its nature is definitely more accidental. Furthermore, Bartlett's test of sphericity was significant (p = 0.000), which indicates that the selected variables are statistically significantly related to each other and can be reduced to a smaller number of components, while the Kaiser–Meyer–Olkin test result was 0.712, which according to the Kaiser (1974) classification suggested that the sampling is adequate on the average (middling) level. Finally, both the temperature parameters and landscape metrics datasets were standardized with the use of the common z-score formula. Then, on the basis of a verified, 10-variable dataset, PCA was performed with a correlation matrix for component loadings (Cload) calculations, while a number of significant components was selected using an empirical Kaiser criterion. In consequence, all components with eigenvalues under 1.0 were removed from further interpretation (Ptak et al. 2018). Component loadings equal to or greater than 0.50 were considered as significant correlations, similar to Żelazny et al. (2018). Here, it must be emphasized that PCA conducted in the current study is based on the small sample size, which is not favorable, as sample sizes above 100 are generally recommended (Gorsuch 1983). However, Shaukat et al. (2016) suggested that in the case of experimental, field-based, and costly investigations, a smaller number of samples can be adopted in the PCA, especially when only the first two principal components are interpreted and visualized (Żelazny 2012).

In the extension to PCA, a linear regression models were performed on the basis of daily water and air temperature data. Such relationships, as suggested by O’Driscoll and DeWalle (2006), indicate the relative control of the thermal regime of watercourses (e.g. if water temperature is determined in larger degree by groundwater inflows or meteorological heat fluxes); however, it is worth noting that they are also commonly used to predict water temperature in various temporal scales (Webb et al. 2003; Lagergaard Pedersen and Sand-Jensen 2007) despite their limited modeling capability (Zhu et al. 2018). Regression models were established for all MR and TR sites using daily mean air temperature from the Warsaw Okęcie meteorological station (operated by the Institute of Meteorology and Water Management—National Research Institute), being the only station located in the central part of the study area (Fig. 1). Statistical analysis was conducted in the Statistica 12.0 and IBM SPSS statistical software.

Meteorological background for the study period was provided by daily mean air temperature and daily precipitation data from the Warsaw Okęcie meteorological station, while mean daily streamflow data was acquired from three available gauging stations located in the downstream parts of the investigated rivers (hydrological stations Piaseczno, Wólka Mlądzka, and Krubice).

Results and discussion

Hydrometeorological background

The summer of 2017 was characterized by a typical air temperature course, with maximum daily mean temperature occurring in August (27.4 °C) and minimum in October (3.1 °C). Mean air temperature in the investigated period was 16.0 °C, which was 0.5 °C higher than the mean value for multiyear period 1988–2017. The sum of precipitation almost reached 482 mm, which was approximately 135 mm higher than the multiyear period mean (1988–2017). Although the investigated period can be considered as relatively wet, nearly half of the precipitation sum was noted in September and October, while from June to August the precipitation sum was similar to the multiyear mean. As a result, high streamflow variability was observed, which makes the obtained results quite representative, as water temperature variability was captured in diverse meteorological and hydrological conditions (Fig. 2).

Fig. 2
figure 2

Mean daily air temperature and precipitation at the Warsaw Okęcie station (a) and mean daily streamflow rate in Piaseczno (Jeziorka River), Wólka Mlądzka (Świder River), and Krubice (Utrata River) (b) during the investigated period

Spatial variability of water temperature

During the summer half-year of 2017 the water temperature varied between sites and clear temperature contrast were found in the case of mean, maximum, minimum, and daily range parameters, as indicated by statistical distribution and cluster analysis (Figs. 3, 4). Generally, in the quasi-natural, undisturbed TR scale sites water temperature exhibited definitely greater spatial variability than in the MR scale sites. Mean water temperature in TR sites varied from 12.8 °C in site T7 to 15.4 °C in site T2, with a higher interquartile range of mean values observed in larger catchments (Fig. 3). In comparison, the mean temperature in MR scale sites was always higher than in TR sites, in the summer of 2017 ranging from 15.5 to 17.2 °C. Cluster analysis supports this findings and permitted to distinguish three main groups, which represent all MR scale sites, colder, upstream TR scale sites with lower interquartile range of mean values and warmer TR sites (Fig. 4). Maximum water temperature was uniform throughout MR sites and only a 1.5 °C difference was noted in both mean and absolute values. In the contrary, in the TR scale sites higher variability of maximum water temperature was observed. Therefore, the biggest cluster was combined from nearly all MR sites and two TR scale catchments (T2 and T3), characterized by wide interquartile range of maximum values, as well as high absolute and mean values of maximum temperature (even up to 26.2 °C and 17.9 °C, respectively). Most of the TR scale sites were linked into second cluster (with moderate mean and absolute values of the maximum), while the smallest group was combine from four headwater TR scale sites (with the lowest maximum values). In the case of the minimum values there was observed similar tendency to the mean temperature, as TR scale sites exhibited lower mean and absolute values of minimum water temperature parameter than MR scale sites (Fig. 3). Such a different statistical distribution of minimum values was clearly reflected by CA, which allowed to distinguish two groups, one linked from all MR sites (with the highest minimum values), while the second from all TR sites (with lower interquartile range of minimum values, as well as absolute minimum values). Finally, daily range parameter was the most varied across sites; in the upstream parts of MR (R1, R4, R7) it was similar to selected TR sites (e.g. T1, T5, T9), representing significant daily temperature fluctuations, while downstream parts of the MR (e.g. R2, R6, and R9) were characterized by small daily temperature fluctuations, similar to headwater TR sites, such as T7, T8, T11. This was confirmed by CA, which permitted to distinguish two groups (Fig. 4).

Fig. 3
figure 3

Distribution of mean (a), maximum (b), minimum (c), and daily range of water temperature (d) during the investigated period

Fig. 4
figure 4

Dendrograms showing the clustering results by the Ward method for daily mean (a), maximum (b), minimum (c), and range (d) of water temperature. Vertical dashed lines indicate optimal number of clusters

Observed thermal heterogeneity in lowland watercourses is quite consistent with findings from the previous studies, which reported higher variability of water temperature in the TR scale than in the MR scale in other geographical regions (Kanno et al. 2014; Dick et al. 2015). Thermal contrasts between tributaries are generally related to their different land use, morphology, and groundwater inflows, which determine vulnerability to atmospheric heating (Imholt et al. 2013b). In contrast, most of the MR sites illustrated relative thermal homogeneity, which was found not only in the river length (e.g. in site R5 and downstream R6), but also between catchments, as documented in sites R1 and R4. Such similarity of water temperature in the MR scale should be associated with the homogeneity of the study area in terms of climatic, hydrological, and hydrogeological factors, considered as overriding in terms of thermal regime control (Hannah and Garner 2015). Furthermore, in the most downstream MR sites there is expected an “averaging effect”—with the rise of the drainage area water with different thermal properties is mixed together, and, in consequence, thermal contrasts decrease, even between different catchments. An analogous tendency was suggested previously in terms of the flow regime—due to an increase in size, streams and rivers are expected to maintain a more homogeneous flow regime if they remain in one hydroclimatic region (Poff et al. 2006).

The results demonstrated that among the MR sites the largest water temperature deviations were related to direct anthropogenic influences. This was observed mainly in site R3, where inflow from upstream pond caused an increase of mean temperature and a decrease of its daily fluctuations; such an effect was reported previously by Dripps and Granger (2013) and Ignatius and Rassmusen (2016). Furthermore, in site R8, thermally stable sewage inflow (during investigated period its temperature was ranging from 18 to 24 °C, while maximum daily range was only 2 °C) dampened diurnal dynamics of water temperature and this effect was additionally strengthened by simultaneous impact of reservoir release. Therefore, in this case heating effect of sewage inflow could not be quantified, however, Kinouchi et al. (2007) stated previously that stream water temperature is not very sensitive to sewage inflows during summer, mainly because of that inflow temperature is similar to those of natural flow. Finally, in site R9 the water temperature was elevated as a result of channel regulation and the removal of riparian vegetation. It must be emphasized that similar effects of human impacts on stream water temperature were exhaustively documented in literature (Webb and Walling 1993; Bae et al. 2016).

Environmental drivers of stream water temperature

The results of the PCA analysis indicate that variability of water temperature during the summer half-year of 2017 for both the TR and MR scales was represented by two independent components (factors) with eigenvalues higher than 1.0, which explained in total 74.6% of the variance. The third (PC3) and fourth (PC4) principal components described 8.6% and 7.7% of the total variability, respectively, but their eigenvalues were less than 1.0, so consequently, they were not analyzed. The remaining percentage can be explained by other factors, which alter water temperature of the watercourses and simultaneously were not included in the analysis (Table 2).

Table 2 Component loadings (Cload), eigenvalues, and explained variance obtained in the PCA analysis

The first principal component (PC1) explained 56.0% of the variability in the data set. Generally, in the case of the investigated watercourses, as indicated by loadings of the PC1, higher catchment area (Cload = 0.94), as well as lower mean catchment slope (Cload = − 0.87) were reflected in higher mean (Cload = 0.98), mean maximum (Cload = 0.86), and mean minimum water temperature (Cload =  0.88), and lower duration of optimum temperature for coldwater fish guild (Cload = − 0.96) (Fig. 5). The second principal component (PC2) explained 18.5% of the variability in the data set and was associated with the thermal regime dynamics. According to this PC, higher mean daily range of water temperature (Cload = 086) was related with higher W:D ratio of the channel (Cload = − 0.61), and simultaneously, with higher mean elevation (Cload = 0.65).

Fig. 5
figure 5

Component loadings associated with all thermal and landscape metrics, used in the PCA analysis

Overall, the first two principal components provide a general view on the spatial variability of the water temperature, as its longitudinal changes across both the MR and the TR sites demonstrated some interesting tendencies, which were clearly related to landscape properties. As highlighted by the principal component PC1, an increase of mean, mean maximum and mean minimum water temperature was observed in the downstream direction, with a simultaneous increase of the catchment area (A) and W:D ratio (W:D), and decrease of the mean catchment slope (S). Here, it must be emphasized that the catchment area of lowland watercourses was not only considered as a proxy for stream flow rate (Hrachowitz et al. 2010), but its changes were a proxy for longer residence time of water and less sensitivity to riparian shade due to the increasing width of the channel. However, as water moves downstream and the streamflow rate greatly increases, mean and maximum temperature increase was limited due to increasing thermal capacity, as suggested by Webb and Walling (1986) and Imholt et al. (2013b); moreover, in some sites a reduction of maximum temperature was observed, which could be caused by inflows from colder tributaries (Poole and Berman 2001). A second tendency was demonstrated by the principal component PC2, where the highest maximum temperature and the highest daily fluctuations (indicated by Mean_rg and Max_rg) were characteristic for the largest catchments from the TR scale (e.g. T1, T2, T3), as well as for upstream sites of the MR. This tendency could be associated with a relatively low thermal capacity and sufficiently wide channels in such sites (high W:D ratio), which are susceptible to radiation heating and cooling (Malcolm et al. 2004). A clear distinction in the PCA bi-plot was noted between such catchments and headwater sites from the TR scale (e.g. T7, T10, T11) and the largest MR sites (e.g. R2, R5, R6), where daily fluctuations of water temperature were dampened due to groundwater inflows and significant thermal capacity, respectively (Fig. 6). For a similar reason, the minimum water temperature values were also the highest in those catchments, as they were negatively correlated with variability parameters. In consequence, the projection of certain watercourses on the plane defined by principal components PC1 and PC2 exhibited clear longitudinal zonation of water temperature (e.g. changes of magnitude and variability temperature parameters), related with the transition from small, upstream headwater TR sites, through downstream TR sites and upstream sites of the MR, to the biggest, downstream MR catchments (Fig. 6). The presented concept is quite consistent with the general concept of changes in stream temperature (Caissie 2006), however, it provides much more details about its spatial dynamics in the investigated lowland watercourses, which could be important in the ecological context. For example, upstream TR sites exhibited the best thermal conditions for growth and activity of coldwater fish, however, in the face of a rapid increase of the maximum temperature values above the upper incipient lethal threshold (25 °C) only a few kilometers downstream, great carefulness is needed in appropriate fish stocking in such catchments. Therefore, cold tributaries of the lowland rivers should be identified by fisheries managers, as they are considered to be important thermal refugia for fish and macroinvertebrates (Lessard and Hayes 2003). From this point of view, the presence of small reservoirs in a lowland landscape (such as upstream from site R3) is undesirable not only because of the interruption of the natural thermal continuum, but also due to the limited possibility of fish migration within main rivers and into the tributaries during stressful thermal conditions (Kanno et al. 2014).

Fig. 6
figure 6

Bi-plot presenting projection of certain cases (watercourses) on the factor plane in the PCA analysis

PCA showed that in the case of the investigated lowland watercourses no altitude zonation was found, as presented by Hrachowitz et al. (2010) and Żelazny et al. (2018) in the case of a mountainous landscape, which was evidenced mainly by the moderate component loading of mean elevation metric (E) in the PC2. Due to small differences in absolute values of mean elevation it could not be responsible for distinct air temperature variability (Isaak and Hubert 2001; Jackson et al. 2017), so obtained component loading was only the result of the accidental spatial distribution of catchments. Furthermore, although many studies demonstrated a moderating effect of the afforestation degree on the river water temperature in the summer period (Imholt et al. 2013a; Kristensen et al. 2015; Ptak 2017; Dugdale et al. 2018), contribution of forested area (F) metrics was not adopted as significant variable in the first two PC. Riparian shade degree (R), despite not being included in the PCA, was also not related with thermal parameters, as streams characterized by similar riparian shading (both with no shade or with high degree of shade) were definitely different in their thermal properties, such as T6 and T9, as well as T7 and T8, respectively. Nevertheless, in the investigated lowland watercourses, as documented by the field observations, most of the shade was associated with ephemeral shrubs and grasses, which could not be quantified from land-use/land cover maps. Moreover, dense communities of aquatic plants in selected stream and river channels were also observed, which were previously documented to affect water temperature (Willis et al. 2017). Thus, achieving a significant relationship between thermal parameters and the degree of channel shading, indicated by afforestation and riparian shade degree, as well as evaluation of the various types of riparian and aquatic vegetation impact on the seasonal modification of thermal characteristics, require more advanced and accurate methods for shade quantification. One of the possibilities is to assess shade associated with ephemeral shrubs and grasses with the use of unmanned aerial vehicles (UAVs), which are commonly used in various aspects of environmental research (Dufour et al. 2013; Zmarz 2014). Much attention should be also paid by the land managers to the preservation such natural riparian buffer zones of lowland watercourses, overgrown with trees, shrubs, and grasses, due to its significance in the context of counteracting global warming (Broadmeadow et al. 2011), expected to have a major impact on the distribution of cold-water fish due to an increase in water temperature (Almodovar et al. 2012; Isaak et al. 2018).

Finally, the results given by the PCA were complemented by the regression analysis. Linear regression models, established between mean daily air and water temperatures during the summer of 2017, exhibited varied parameters and explanatory power, and were statistically significant at p < 0.01. The explanatory power (R2) of regression models for all sites ranged from 0.73 to 0.88; however, greater variability of R2 was found in the TR scale (0.73–0.87), particularly in comparison to the MR scale (0.81–0.88) (Table 3). Such a tendency was also observed in the case of regression parameters. The values of the slope and intercept of the MR sites, excluding those with clear anthropogenic alternations (R8, R9), were generally similar, while in the case of the TR they were definitely much more varied (Table 3). Overall, for both the TR and the MR measurement sites, an increase of the regression slope resulted in a decrease of the intercept (Fig. 7). Such a shift of regression parameters was clearly associated with the increase of the catchment area and the decrease of its slope, resulting in a longer water residence time, and in consequence, more vulnerability for atmospheric heating. Thus, the smallest and the coldest TR catchments (T6, T7, T8, T10, T11), which have lower slopes and higher intercepts, were less meteorologically controlled and in turn their thermal regime reflected more significant groundwater impact, as it was highlighted by Erickson and Stefan (2000) and Krider et al. (2013). In contrast, the thermal regime of the largest catchments (with the highest mean and mean maximum temperature as indicated by PCA) was driven mostly by short-wave and long-wave atmospheric heat fluxes (Benyahya et al. 2012), as they were characterized by lower regression slopes and higher intercepts. In addition, in such meteorologically controlled sites the explanatory power of the regression measured by R2 was generally better, which indicates a closer relationship between air and water temperatures (Erickson and Stefan 2000) (Table 3).

Table 3 Summary of linear regression models between daily mean air and water temperatures for both MR and TR sites
Fig. 7
figure 7

Relationship between slope and the intercept of the linear regression performed for all of the investigated sites

Conclusions

The present study provided a new insight into spatial water temperature variability across lowland temperate watercourses and examined the effects of environmental properties on their thermal heterogeneity. The combination of three relatively simple statistical methods, namely cluster analysis, principal component analysis, and regression analysis, proved to be helpful in the identification and visualization of general spatial water temperature behavior of lowland watercourses. They revealed clear longitudinal zonation of their thermal properties, driven mainly by changes of landscape properties, such as catchment area, its mean slope, and W:D ratio of the channel. Moreover, the results suggested that the relative importance of landscape drivers in shaping thermal heterogeneity of lowland watercourses seems to be the highest in the source zone, while as water moves downstream, meteorological conditions and streamflow rate have dominant control over the thermal regime, which was indicated by the parameters of the linear regression, established between daily mean air and water temperatures. In consequence, great variability of the water temperature parameters was found in the tributary scale catchments, covering up to 30–40 km2, especially in terms of maximum and daily range water temperature values, while further, in the downstream direction, main river sites were characterized by similar thermal properties, even between separate, independent catchments. Among tributaries two groups of catchments were clearly identified—headwater sites with the lowest mean and maximum temperatures and the largest tributaries, more similar in their thermal parameters to the main rivers. In the main river sites, the highest thermal contrast were related to anthropogenic impacts, such as reservoir releases and sewage inflows. Due to the fact that the spatial water temperature behavior explored in current study is expected to appear in similar Polish lowland catchments, the obtained results provide a unique basis in the context of fisheries management and land management practices, as well as for further investigations on the thermal regime of lowland temperate rivers.