Variability of Mean Annual Flows in Southern Quebec (Canada)

: Snow is the main source of streamﬂow in temperate regions characterized by very cold and snowy winters. Due to global warming, these regions are experiencing a signiﬁcant decrease in snowfall. The main objective of this study is to analyze the impacts of snowfall on the spatio-temporal variability of mean annual ﬂows (MAFs) of 17 rivers, grouped into three hydroclimatic regions, from 1930 to 2019 in southern Quebec. In terms of spatial variability, snowfall is the variable most correlated with MAFs (positive correlation), followed by drainage density (positive correlation) and wetland surface areas (negative correlation). Due to the inﬂuence of these three factors, MAF values are generally higher in the most agricultural watersheds of the southeastern hydroclimatic region on the south shore than in the less agricultural watersheds of the southwestern hydroclimatic region on the north shore of the St. Lawrence River. As for temporal variability, the four statistical tests applied to the hydrological series detect no signiﬁcant downward trend in MAFs, despite having reduced snowfall. Instead, they suggest an evolution toward an increase in mean annual ﬂows, as a result of increased rainfall due to the increase in temperature. This evolution is more pronounced on the north shore than on the south shore, likely due to the presence of wetlands and others water bodies, whose runoff water storage capacity does not change over time to be able to store the surplus of the quantity of water brought by the increase in rain.


Introduction
In cold temperate regions characterized by very cold and snowy winters, one of the major impacts of global warming is a significant decrease in snowfall. This decrease has already been observed in many parts of the world [1][2][3][4][5]. From a hydrological perspective, snow is the main source of streamflow and aquifer recharge in these regions. The consequences of this decrease in snow cover have mainly been analyzed in relation to the temporal variability of extreme flows (floods and low flows). For example, Blöschl et al. [6] clearly demonstrated that the decrease in snow cover resulted in a significant decrease in the magnitude of spring freshets caused by snowmelt in eastern Europe. Similar impacts have been observed in other regions [7] and predicted by all climate models (e.g., [8]).
However, none of these studies analyzed the impacts of this decrease in the magnitude of floods on mean annual flows. In general, while studies exist in the scientific literature on extreme flows, very few have focused on the impacts of climate change on mean annual flows (e.g., [9][10][11]). However, in the current context of global warming, mean annual flows can be used to detect the impacts of changes in temperature and precipitation regimes on river flows as a result of this warming [12,13]. In this context, their analysis is fully justified along with that of extreme flows.
Nevertheless, very few studies have focused on temporal variability in general and the impacts of global warming in particular on mean annual flows in Canada [28,29]. Others have analyzed this temporal variability to determine the climatic indices that influence it (e.g., [30][31][32][33][34]). Furthermore, there have been no studies on the influence of snowfall on the spatio-temporal variability of these flows. Moreover, on the methodological level, none of these studies analyzed the effects of persistence, either in the short (STP) or the long term (LTP), on the stationarity of these average annual flows.
In light of these considerations, this study aims to fill these gaps by pursuing the following objectives: 1.
Identify the physio-climatic factors and the factors related to land use and land cover that influence the spatial variability of mean annual flows in southern Quebec.

2.
Analyze the temporal variability of mean annual flows to determine the influence of short-term persistence (STP) and long-term persistence (LTP) effects on their stationarity in the current context of global warming. This objective aims to test the hypothesis that the significant decrease in snow cover led to the decrease in mean annual flows in southern Quebec. This is the primary objective of this study.

3.
Determine whether the factors influencing spatial variability can amplify or mitigate the effects of global warming on the temporal variability of these flows. 4.
Determine the influence of spring flows resulting primarily from snowmelt on the spatio-temporal variability of mean annual flows in the various hydroclimatic regions of southern Quebec.
This study is part of the research program aimed at determining the impacts of climate change on the spatio-temporal variability of flows at different time scales in relation to land cover and land use in southern Quebec.

Description of the Watersheds
This study is based on an analysis of mean annual flows measured continuously from 1930 to 2019 in 17 watersheds ( Figure 1 and Table 1). These flows were minimally disturbed by human activity. These watersheds were grouped into three homogeneous hydroclimatic regions [22] on both shores of the St. Lawrence River in southern Quebec. The first hydroclimatic region (southwest hydroclimatic region) is located on the north shore, and the rivers in its watersheds cut through the Canadian Shield, which mainly consists of metamorphic and igneous rocks covered by glaciofluvial deposits. The climate is continental temperate. The forest consists mainly of sugar maple trees (sugar mapleyellow birch forest and basswood-sugar maple forest). The second hydroclimatic region (southeastern hydroclimatic region) is located on the south shore, south of 47 • N. The rivers of this region drain the geological formations of the Appalachians and the St. Lawrence Lowlands, which consist mainly of sedimentary rocks covered in sediment of different origins (glacial, lacustrine, marine, etc.). Forest formations consist of sugar maple trees (sugar maple-hickory stands, basswood-sugar maple forests, and sugar maple-yellow birch forests). The climate is mixed temperate with maritime and continental features. North of this parallel on the same shore is the third hydroclimatic region (east hydroclimatic region), which has watersheds contained in the same two geological formations. However, the sugar maple forest is gradually replaced by the fir forest (balsam fir-white birch stands and balsam fir-yellow birch stands). The climate is becoming increasingly maritime at higher latitudes.   20). Finally, physiographic data and data related to types of land use and land cover were drawn from the Université du Québec à Trois-Rivières GEOLAB laboratory [26,35]. These data were described by [24]. It is important to mention that measurements of the daily flows of the Du Nord, Eaton, Matawin and Blanche rivers began in 1931, 1932, 1933 and 1934, respectively. For the Matawin River,  Finally, physiographic data and data related to types of land use and land cover were drawn from the Université du Québec à Trois-Rivières GEOLAB laboratory [26,35]. These data were described by [24]. It is important to mention that measurements of the daily flows of the Du Nord, Eaton, Matawin and Blanche rivers began in 1931, 1932, 1933 and 1934, respectively. For the Matawin River, data from 1940, 1941 and 1943 were incomplete. This is also the case for 1940 on the Matane River.

Analysis of Spatial Variability of Data
This analysis consisted of three main steps: In the first step, annual means for annual flows (means calculated over the 365 days of each year) for each river from 1930 to 2019 were calculated using daily flows. These means were calculated in two ways: from October to September and from January to December. The approaches led to similar results. These means were then compared using parametric (ANOVA) and non-parametric (Kruskal-Wallis) tests. It is important to note that at this stage, flows measured in m 3 /s were converted into specific flows expressed in L/s/km 2 to eliminate the influence of watershed size on the results.
In the second step, these means were correlated with physiographic and climatic variables and types of land use and land cover in watersheds. Linear correlation (parametric) and Spearman rank (non-parametric) methods were applied. As both methods led to the same results, only those obtained from the first method are presented.
In the third step, mean annual flows were compared to mean seasonal spring flows. Note that in Quebec, due to its cold temperate climate, mean annual flows are strongly influenced by snowmelt in the spring. Snowmelt is the main source of streamflow [36] and aquifer recharge [37,38]. For this comparison, two ratios were calculated: the ratio of mean annual flows (MAFs) to spring mean daily flows (SMDFs) and the ratio of mean annual flows to spring maximum daily flows (SMAXDFs). These two ratios, expressed as percentages, were labelled as R-Mean and R-Max, respectively. They measure the contribution of the spring daily mean flows and the annual daily maximum flows generated mainly by snowmelt in Quebec on the mean annual runoff (MAF). If the values of these two ratios are high, the contribution of spring flows to annual runoff is relatively high. Thus, the annual flow tends to concentrate during the spring season when the snow melts. This thus becomes the main source of this annual flow. To complete this comparison, the correlation between mean annual flows and the two spring flow series measured from April to June was calculated. The purpose of this correlation calculation was to identify the type of spring flows that could potentially be used to estimate mean annual flows based only on those measured in spring.

Analysis of the Temporal Variability of the Data
To analyze the long-term trend for mean annual flows, five different statistical trend analysis tests were carried out.
First, the trend-free prewhitening (TFPWMK) method was applied to eliminate autocorrelation effects on the long-term trend. It was described in detail by [39]. The second test, the modified Mann-Kendall test (MMKH), was applied to eliminate all autocorrelations; however, it was completed by correcting the variance of the hydrological series [40]. The third test, long-term persistence (LTP), eliminated the effects of long-term persistence (Hurst effect [41], while the first two tests only eliminated the effects of short-term persistence (STP). The purpose of the fourth test was to detect the dates of the shift by means of the hydrological series ( [42]. The rationale for selecting this test is that it can detect abrupt or gradual breaks in means, unlike all other statistical tests used in hydrology. These four tests are widely described in the scientific literature and are commonly used to analyze the long-term trend and breaks by means of the hydrological series (e.g., [43][44][45][46]). Note that other tests based on Bayesian statistics (e.g., [19,21]) lead to the same results as the most widely used tests described previously. The fifth and final test was applied to verify the spatial autocorrelation of flows. It has been described in detail by [19,47]. This test (no spatial autocorrelation) confirmed the results of the previous four tests; the results of the fifth test will therefore not be presented here.

Spatial Variability of Mean Annual Flows
The mean annual flow values are presented in Figure 2. They vary spatially in southern Quebec even if the differences between the watersheds are not very important. They ranged from 14 L/s/km 2 (Châteauguay River) to 25.5 L/s/km 2 (Blanche River). Generally speaking, in the two hydroclimatic regions in the southwest and east, these values were generally less than 20 L/s/km 2 , except for a few rivers. However, they exceeded this threshold in the southeastern hydroclimatic region-the most agricultural region in Quebec. The application of ANOVA and Kruskal-Wallis tests revealed that the means of the flows of 17 rivers are globally significantly different at the 5% threshold.
Water 2022, 14, x FOR PEER REVIEW 5 of 13

Spatial Variability of Mean Annual Flows
The mean annual flow values are presented in Figure 2. They vary spatially in southern Quebec even if the differences between the watersheds are not very important. They ranged from 14 l/s/km 2 (Châteauguay River) to 25.5 l/s/km 2 (Blanche River). Generally speaking, in the two hydroclimatic regions in the southwest and east, these values were generally less than 20 l/s/km 2 , except for a few rivers. However, they exceeded this threshold in the southeastern hydroclimatic region-the most agricultural region in Quebec. The application of ANOVA and Kruskal-Wallis tests revealed that the means of the flows of 17 rivers are globally significantly different at the 5% threshold. The correlation analysis between mean annual flow values and physiographic variables showed that flows were positively correlated with drainage density but negatively correlated with the surface areas of wetlands in watersheds ( Table 2). For climate variables, flows were only positively correlated with snowfall. The highest correlation was observed between flows and snowfall in the winter/spring season.  The correlation analysis between mean annual flow values and physiographic variables showed that flows were positively correlated with drainage density but negatively correlated with the surface areas of wetlands in watersheds (Table 2). For climate variables, flows were only positively correlated with snowfall. The highest correlation was observed between flows and snowfall in the winter/spring season.

Temporal Variability of Mean Annual Flows in Quebec
The results of the four statistical tests are presented in Table 3. The three Mann-Kendall tests produced similar results and were supported by the Lombard test. Two aspects must be highlighted: the signs of the Z-score values of the tests and their significance. For the southwest hydroclimatic region on the north shore, all Z-scores were positive. This reflected an evolution toward an increase in mean annual flows in this region ( Figure 3). Nevertheless, this trend was statistically significant for only three rivers (Figure 2). The same evolution was also observed for almost all rivers in the southeast hydroclimatic region with Z-scores that were also positive. However, only the Z-score for the Châteauguay River was significant at the 5% threshold ( Figure 2); the Nicolet Sud-Ouest River and the Beaurivage River were only significant at the 10% threshold. Moreover, the LTP test did not detect any significant trends for the two last rivers, even at this threshold. In this hydroclimatic region, only the Rivière Du Sud had a negative Z-score. A negative Z-score was observed for over half of the rivers in the east hydroclimatic region. This result suggests an evolution toward a decrease in mean annual flows in this hydroclimatic region, in contrast to the other two hydroclimatic regions. Moreover, no trends were statistically significant, even at the 10% threshold. The Lombard test detected breaks in means for three rivers in the southwestern hydroclimatic region and for only one river (Châteaugay river) in the southeastern hydroclimatic region. Apart from the Matawin River, breaks in means of mean annual flows occurred during the second half of the 1960s in both hydroclimatic regions.

Relationship between Mean Annual Flows and Spring Flows
Note that in cold temperate regions characterized by very snowy winters, the annual runoff depended heavily on the melting of snow accumulated during the cold season. However, many factors were at play. To determine the spatial variability of this influence, the relationship between mean annual flows and spring flows was analyzed. Two variables were considered for these types of flows: spring mean daily flows (SMDFs) and spring maximum daily flows (SMAXDFs). The latter variable frequently corresponds with the magni-

Relationship between Mean Annual Flows and Spring Flows
Note that in cold temperate regions characterized by very snowy winters, the annual runoff depended heavily on the melting of snow accumulated during the cold season. However, many factors were at play. To determine the spatial variability of this influence, the relationship between mean annual flows and spring flows was analyzed. Two variables were considered for these types of flows: spring mean daily flows (SMDFs) and spring maximum daily flows (SMAXDFs). The latter variable frequently corresponds with the magnitude of the annual freshet usually generated by snowmelt.
The correlation coefficients calculated between the mean annual flows and the two spring flow variables are presented in Table 3. For spring mean daily flows, correlation values were relatively lower in the southeast hydroclimatic region-the region that is the most agricultural-than in the other two hydroclimatic regions. In fact, the value was less than 0.700. This threshold was exceeded in the southwest hydroclimatic region on the north shore. For the last hydroclimatic region in the east, the correlation coefficient values varied much more than those in the other two hydroclimatic regions. They ranged from 0.65 (Rivière des Trois Pistoles) to 0.97 (Blanche River). For the second variable, the same trend was observed. Correlation values were lower in the southeast hydroclimatic region than in the other two regions. It is important to note that correlation coefficient values were much lower than those calculated with spring mean flows. Thus, mean annual flows were very weakly correlated with spring maximum daily flows in southern Quebec.
To compare mean annual flows with spring flows, ratios between the mean annual flows and the two spring flow variables were calculated. These ratios are presented in Table 4. For spring mean flows, the ratios varied between 39.3% (Blanche River) and 65.7% (Châteauguay River). Table 4 clearly demonstrates that the values of these ratios were generally greater than 50% in the southeast hydroclimatic region, 50% on average in the southwest hydroclimatic region on the north shore, and less than 45% in the east hydroclimatic region. As for spring maximum daily flows, all ratios were less than 15%, except for those of the Vermillon River (19%) on the north shore and the Etchemin River (27%) on the south shore. No spatial disparity was observed in relation to these ratios. To determine the factors that influenced the spatial variability of these ratios, correlation coefficients between these ratios and physio-climatic variables were calculated ( Table 2). The ratios between mean annual flows and spring mean daily flows were negatively correlated with forest surface areas and winter-spring snowfall but positively correlated with rainfall and daily maximum temperatures. These variables were very strongly correlated with the ratios. As for the ratios between mean annual flows and spring maximum daily flows, no correlation coefficients were statistically significant at the 5% threshold.

Analysis of Spatial Variability of Mean Annual Flows
The main objective of this study was to analyze the influence of snow cover on the spatio-temporal variability of mean annual flows in southern Quebec from 1930 to 2019. In terms of spatial variability, the correlation analysis between these flows and many physiographic and climatic variables revealed that snowfall was the variable most correlated with flows, followed in succession by watershed drainage density and wetland surface area.
The positive correlation observed between flows and snow cover is entirely normal because, due to the cold temperate climate in Quebec, annual runoff depends mainly on snowfall that melts in the spring. Thus, the more snowfall in a watershed, the higher the annual runoff. In Quebec, snowfall was generally higher on the south shore than on the north shore; thus, it was the same for mean annual flows. In addition to snowfall, the second factor influencing spatial variability of flows was watershed drainage density. The influence of this factor was reflected in the faster transfer of runoff to river channels, resulting in increased flows. Since the modernization of agricultural practices, which began in 1950 [48] in Quebec, many canals have been dug to drain wetlands with surface areas that have decreased significantly (<5%). In addition, river channel adjustments (avulsions) were made in agricultural watersheds, most of which were located on the south shore, south of 47 • N in particular. In this region, the most agricultural region of Quebec (agricultural surface areas exceeding 10%), mean annual flow values were the highest due to the increase in this substantial amount of drainage. However, it is important to note that mean annual flows were not significantly correlated with agricultural surface areas as would normally be expected due to soil sealing. This lack of correlation between the two variables could be explained in part by the fact that the modernization of agriculture led to a significant reduction in cultivated land and an increase in uncropped land and reforested areas. These land use changes have had the effect of promoting water infiltration at the expense of runoff. This resulted in a significant increase in minimum flows at the expense of maximum flows in agricultural watersheds (e.g., [21]). However, despite this increase in infiltration, runoff remained consistently higher in this region-the most agricultural in nature-than in the other two hydroclimatic regions.
In contrast, intensive drainage into agricultural watersheds on the south shore resulted in a significant decrease in wetlands and others water bodies (lakes) with surface areas that now occupy less than 4% of the watersheds, while they occupied more than 8% of the watersheds on the north shore. The presence of larger wetland and other water bodies surface areas on the north shore resulted in significantly less runoff during snowmelt because the water was stored in large quantities in these wetlands and water bodies (e.g., [49][50][51][52][53][54][55][56]). This explains the negative correlation observed between wetland surface areas and mean annual flows. Thus, mean annual flows were generally lower on the north shore than on the south shore, partly due to the presence of these wetlands and other water bodies.
Furthermore, snowfall in the winter/spring season was also significantly correlated with the ratios between mean annual flows and spring mean daily flows. This correlation is negative because these ratios were higher south of 47 • N (>50%) than north of this parallel (<45%) on the south shore, whereas the winter-spring snowfall was higher north of this parallel than south of it. However, these ratios were correlated with the amount of rainfall during the same season because rainfall was greater south of this parallel than north of it. This spatial variability in rainfall was due to that of maximum daily temperatures. They were also higher south of 47 • N than north of it on the south shore due to the influence of latitude. They were strongly correlated with flow ratios. Finally, the negative correlation between these ratios and forest surface areas can be explained by relatively smaller forest surface areas south of 47 • N than north of it due to agriculture, which promoted surface runoff. From a hydrological perspective, these ratios indirectly reflected the contribution of spring flows generated by snowmelt to total annual runoff. The significant influence of these spring flows resulted in a low ratio between mean annual flows and spring mean daily flows. This was the case in the eastern hydroclimatic region, where the influence of spring flows on total annual runoff exceeded 60%, whereas in the southeastern hydroclimatic region located to the south, it was less than 50% due to the increasing influence of summerfall rainfall to annual runoff in the region.

Analysis of Temporal Variability of Mean Annual Flows
It is important to remember that all previous studies have shown a significant decrease in the amount of snow in Quebec ( [14][15][16], like many cold temperate regions. This decrease was confirmed within the framework of this study by analyzing the data available in a few watersheds (these results are not presented here). In these cold temperate regions where the annual runoff strongly depends on the amount of snow, it has been shown that the impacts of the decrease in this amount of snow have caused a significant drop in these flows, in particular the magnitude of the floods (see [6]).
However, analysis of the stationarity of the hydrological series using three statistical tests from 1930 to 2019 revealed that the interannual variability of mean annual flows during the period was characterized by an upward trend in flows (positive Z-scores) for 12 of the rivers analyzed (71%) but a downward trend (negative Z-scores) in only five of these rivers. An evolution toward an increase in mean annual flows was observed in the two hydroclimatic regions south of 47 • N on both shores of the St. Lawrence River, while all rivers characterized by an evolution toward a decrease in mean annual flows were located north of this parallel on the south shore. Nevertheless, with regard to the statistical significance of the tests, there was a significant positive trend for only four rivers, three of which were located on the north shore and only one (Châteauguay River) of which was located on the south shore, south of 47 • N. It is important to note that this river's watershed was the most urbanized (urbanized area of >40%) of all the 17 watersheds analyzed. These trend analysis results clearly demonstrate that, despite the significant decrease in the amount of snow, the average annual flows do not decrease significantly over time in southern Quebec. On the contrary, they increase or remain constant over time. The absence of an overall decrease in these flows, despite a decrease in snowfall, can only be explained by the increase in rainfall observed on both shores of the St. Lawrence River. During the hot season (summer and fall), this increase in the amount of rainfall has been well documented throughout the northeastern part of North America [57,58], including Quebec [59,60]. This increase is also observed during the cold season (winter and spring) and its impacts translate into an almost general increase in the magnitude of winter floods and, to a lesser extent, that of spring floods. This increase in the amount of rainfall during the two seasons was confirmed by the analysis of data in a few stations as part of this study. It compensates for the deficit in the quantity of snow, thus making it possible to maintain or increase the annual flow. The increase in the amount of precipitation in the form of rain results from the almost general increase in temperature observed in Quebec and Canada as a whole.
In this context of increased rainfall following warming, the storage capacity of runoff water by wetlands and other water bodies does not change over time. Thus, the excess rainfall runs off, causing an upward trend in mean annual flows over time, which is more pronounced on the north shore (southwestern hydrologic region) than on the south shore [61]. On the latter, the significant decrease in agricultural areas since 1950 has reduced this upward trend in flows due to greater water infiltration, which easily absorbs excess rainfall.
This increase in water infiltration in the most agricultural watersheds on the south shore translates into a significant increase in minimum flows over time, unlike the watersheds on the north shore [24]. It follows that the presence of wetlands and other water bodies would tend to favor the increase in means annual flows over time, while the change in land use which results in a significant reduction in agricultural areas tends to attenuate this increase.

Conclusions
Like many cold temperate regions characterized by very snowy and cold winters, Quebec is experiencing a significant decrease in snowfall and an increase in winter temperatures due to global warming. The main objective of this study was to determine the impacts of these changes in temperature and precipitation regimes on the spatio-temporal variability of mean annual flows from 1930 to 2019.
With regard to spatial variability, the correlation analysis revealed that snowfall was the climate variable that was most positively correlated with mean annual flows. They were also positively correlated with watershed drainage density but negatively correlated with wetland/water bodies surface areas. As a result of these three factors, mean annual flows are generally higher in the most agricultural watersheds of the southeastern hydroclimatic region on the south shore than in those less agricultural of the southwestern hydroclimatic region on the north shore of the St. Lawrence River. On the north shore, the presence of larger wetland and other water body surface areas supported runoff storage from snowmelt, which limited its transfer to river channels. On the south shore, intensive drainage of these wetlands into agricultural watersheds, since the modernization of agricultural practices began in 1950, facilitated the rapid transfer of runoff to river channels.
As for temporal variability, the most significant finding of this study is the upward trend in mean annual flows for most of the rivers analyzed despite the decrease in snowfall. This upward trend is explained by the increase in temperature and rainfall, which compensated for the decrease in snowfall. This increase in temperature and rainfall appears to have had a greater impact on the rivers on the north shore than on those on the south shore, likely due to the presence of wetlands, which are an important source of the moisture needed to maintain rainfall intensities.
This study found that decreased snowfall did not lead to a widespread decrease in mean annual flows in Quebec as expected due to the rise in temperatures and the resulting increase in rainfall. However, no hydrologic and climate models have taken this increase in temperatures and rainfall into account to predict the hydrological impact of decreased snowfall in Quebec. In addition, in their predictions, these hydroclimatic models do not take into account the amplification effects of wetlands and water bodies on the north shore on the temporal variability of means annual flows on the one hand and the attenuation effects exerted by the reduction in agricultural areas on this variability of flows on the south shore on the other hand. This aspect must be taken into account to predict, with fewer errors, the temporal variability of mean annual flows in future decades in southern Quebec.