Compound changes in temperature and snow depth lead to asymmetric and nonlinear responses in landscape freeze–thaw

Cycles of freeze–thaw (FT) are among the key landscape processes in cold regions. Under current global warming, understanding the alterations in FT characteristics is of a great importance for advising land management strategies in northern latitudes. Using a generic statistical approach, we address the impacts of compound changes in air temperature and snow depth on FT responses across Québec, a Canadian province ~ 2.5 times larger than France. Our findings show significant and complex responses of landscape FT to compound changes in temperature and snow depth. We note a vivid spatial divide between northern and southern regions and point to the asymmetric and nonlinear nature of the FT response. In general, the response of FT characteristics is amplified under compound warming compared to cooling conditions. In addition, FT responses include nonlinearity, meaning that compounding changes in temperature and snow depth have more severe impacts compared to the cumulative response of each individually. These asymmetric and nonlinear responses have important implications for the future environment and socio-economic management in a thawing Québec and highlight the complexity of landscape responses to climatic changes in cold regions.

www.nature.com/scientificreports/ In each ecozone, the annual number of Frozen Days (FD year ) and the number of transient days in the winter season (December, January, and February; FTD DJF ) at grid scale are extracted. FD year is defined as the total number of days with frozen soil in both AM and PM of a given day, spanned from September 1st to August 31st. FTD DJF is also defined as the total number of days with the transient state over a typical winter season (months December, January and February). Transient days are defined above as days with frozen soil in AM and thawed soil in PM, or vice versa (thawed in AM, frozen in PM). In line with earlier studies 41,42 , our empirical results also show that these characteristics are the two most sensitive regional FT characteristics to changes in climate conditions, and demonstrate the strongest interdependence with temperature and snow depth at common scales and across the eight ecozones of Québec. These FT characteristics have also direct implications for land management. The changing annual number of frozen days can affect the length of phenological and agricultural activities 43 , accessibility to natural resources 44 , and the emission of greenhouse gasses from thawing permafrost 45 . In addition, the increasing number of transient days can be a proxy for land subsidence and erosion 46 as well as the deterioration rate of civil infrastructures such as buildings, roads, and pipelines 47 .
By applying the proposed framework (see "Methods" below), FT characteristics at each ecozone can be conditioned to air temperature and snow depth using developed copula models 48,49 . Before applying these models for impact assessments, we evaluate their performance in representing the marginal and joint characteristics of the two considered FT characteristics, observed across the eight ecozones. Figure 2 summarizes the results, in which the top row depicts the results for the annual number of frozen days (FD year ) and the bottom row is related to the number of transient days during the winter season (FTD DJF ). In each row, panels from left to right show the expected estimations (bars) and observed values (thresholds) for the mean, standard deviation, and skewness of the considered FT characteristics as well as their dependence with temperature and snow depth across the eight ecozones. Figure 2 clearly shows that the parametrized copulas can effectively represent the first three moments of the empirical probability distributions of both FT characteristics considered. For FD year , the overall expected relative error of 0.2%, 2.2%, and 13.5% is observed for mean, standard deviation, and skewness, respectively. For FTD DJF , the overall expected relative error for the first three moments are 0.4%, 1.1%, and 3.1%. The interdependencies between the considered FT characteristics with temperature and snow depth are also preserved with average relative errors of ~ 1% or less in all cases. Such performances in reconstructing the historical FT characteristics justify the application of the model for assessing the compounding impacts of changing temperature and snow depth.
Apart from benchmarking the performance of copula models in reconstructing the historical dependencies between FT characteristics and climatic conditions, Fig. 2 shows a clear increase in the magnitude of dependencies between FD year and both temperature and snow depth by moving from north to south. Similarly, the magnitude of dependence between FTD DJF and the considered climatic drivers generally increases in southern ecozones. Using the one-way ANOVA test, we assess the differences between the estimated dependencies across different ecozones. The results highlight the uniqueness of interdependencies between FD year , temperature, and snow depth across all ecozones. In addition, dependencies between FTD DJF and climate drivers are almost unique. The only statistically significant similarity is in the case of FTD DJF and temperature in EZ4 and EZ6. Significant variations in the dependencies between FT characteristics and hydroclimatic drivers can reveal the effect of ecosystem conditions on regulating the impacts of changing conditions on FT characteristics. To showcase this  www.nature.com/scientificreports/ empirically, we consider two compound scenarios related to opposing warming and cooling conditions, consisting of simultaneous changes of 2 °C warming in long-term mean temperature along with 10 cm thinner snow on the ground (compound warming), as well as 2 °C cooling in long-term mean temperature along with 10 cm thicker snow on the ground (compound cooling). Figure 3 summarizes the results of this impact assessment and compares them with historical conditions. Whiskers span the range of expected FT characteristics obtained by 1,000 resampling trials-see "Methods" for the details of the proposed impact assessment framework. Historical baseline (no change condition) as well as compound cooling and warming scenarios are shown with gray, blue, and red colors, respectively. Using Fig. 3, a couple of key findings are made. First, there is a clear geographic departure between FT responses to unique compound changes in temperature and snow depth. While the response of FD year to compound cooling and warming is more vivid in the north, the impacts on FTD DJF are more pronounced in the south.
Considering FD year and under compound warming conditions, on average between 48 (EZ1) to 22 (EZ8) fewer days with the frozen condition are expected in a typical year. The north-south decline in the response of FD year is less obvious under compound cooling conditions, ranging from roughly 10 (EZ5) to 34 (EZ2) days extension in frozen conditions per year. Considering FTD DJF , in contrast, the impact of compound warming increases by moving from north to south. Under this scenario and on average, 10 more days during a typical winter with the transient condition is expected in EZ8, where the majority of Québec's population is concentrated, as opposed to only 5 days in EZ1. Under cooling conditions, this change ranges from roughly 2 (EZ3) to 6 (EZ8) fewer days with transient conditions in a typical winter. These spatially heterogeneous impacts on FT characteristics under unique compound changes in temperature and snow depth can reveal how changes in ecozonal features (e.g. vegetation and soil types, land cover, climate, etc.) regulate FT responses to compound changes in climate conditions. Second, Fig. 3 clearly shows that in those ecozones in which FT characteristics are more sensitive to compound changes, responses to opposite warming and cooling scenarios become asymmetric and are more vivid under the compound warming scenario. For instance, looking at FD year , the absolute shift under compound warming scenario is 30 days more in the Northern Arctic (EZ1), compared to the shift under cooling scenario. The higher impact of warming on FD year in northern ecozones is of great importance in the context of the landatmospheric Carbon emissions due to permafrost degradation 50 . Similarly, in the case of FTD DJF , the magnitude of positive shift due to the compound warming scenario is about 4 days more in the Mixed Wood Plains (EZ8) compared to the negative shift caused by the compound cooling scenario. More sensitivity to the compound warmings has important implications regarding soil stability 51 and deterioration of critical infrastructure 52 that are rapidly aging 53 .
To better understand climate controls on FT characteristics, we consider additional compound and individual changes in temperature and snow depth by mixing-and-matching long-term shifts in mean temperature, ranging from −2 °C to +2 °C (sampled every 0.5 °C), and mean snow depth, ranging from − 10 cm to 10 cm (sampled every 2.5 cm). This configuration results in 81 different scenarios, with which FD year in the Northern Arctic (EZ1) and FTD DJF in the Mixed Wood Plains (EZ8) are conditioned. We choose these two ecozones due to their largest www.nature.com/scientificreports/ sensitivity in changing hydroclimatic conditions with respect to FD year and FTD DJF , respectively (see Fig. 3); and accordingly, they can better reveal how marginal and joint variations in temperature and snow depth can result in alterations of FT characteristics. Figure 4 presents the results of this analysis, in which the top (panels a to c) and the bottom rows (panels d to f) depict the results related to FD year and FTD DJF , respectively. Panels a and d show response surfaces for changes in FT characteristics with respect to changes in temperature and snow depth. These response surfaces are reconstructed using the proposed C-vine copula and 1000 resampling-see "Methods"below. Panels b and e, as well as c and f, show the partial derivatives of these response surfaces, derived at known temperature and snow depth conditions, respectively. Looking at the results obtained for FD year across EZ1 (Fig. 4, first row), several interesting observations can be made. First and foremost, the non-symmetric response of FT characteristics to warming vs. cooling again reveals itself and is even more vivid in this round of analysis. While cooler temperature and increased snow depth increase the length of the frozen period in a typical year by around 18 days, warmer temperature and decreased snow depth can decrease the length of the annual frozen period by around 48 days. Projections of response surfaces, shown in panels (b) and (c), clearly present higher sensitivity of FD year to warmer temperatures and thinner snow depth. Second, despite having different signs, the magnitude of interdependencies with temperature and snow depth is more or less the same (− 0.27 and 0.22, respectively; see Fig. 2); and therefore, the sensitivity of FD year to changes in temperature and snow depth are rather similar. This is empirically the case. The expected change in FD year under constant temperature conditions when snow depth is varying from − 10 cm to 10 cm of its long-term average value is about 25 days (Fig. 4b). Similarly, the expected change in FD year under constant snow depth when the temperature is shifted from − 2 °C to + 2 °C of the annual long-term average is about 30 days (Fig. 4c). Having said that, it should be mentioned that the impacts of both snow depth and temperature variations are intensified as the temperature warms and snow depth declines. This is another line of evidence for asymmetric marginal impacts that are further controlled by the state of other driving variables.
A similar rationale can be used to interpret the results obtained for FTD DJF across the Mixed Wood Plains (EZ8; the bottom row of Fig. 4). The changes in the expected FTD DJF range from 9 more transient days under warmer climate with thinner snow depth to 7 fewer days with transient days under cooler climate and deeper snow depth. This observation again points at the asymmetric response of the FT characteristics, although it is less obvious compared to FD year across EZ1. In this case, however, the dependence between FTD DJF and temperature is stronger than the dependence between FTD DJF and snow depth; and therefore, this FT characteristic is more sensitive to temperature changes, compared to changes in the snow depth. As an example, the expected

Discussion
The results provided above illustrate two important findings. First, the impacts of changing temperature and snow depth conditions on FT characteristics differ spatially, with distinct manifestations and natures of response in northern and southern regions. Second, asymmetric responses to compound warming and cooling conditions are seen across different FT characteristics and/or regions. In general, when there is a significant sensitivity to changing climate, the response is more intense under warmer and/or decreased snow depth conditions. In this section, we look at another key feature of FT responses to compound changes in temperature and snow depth, which is rather overlooked in current literature. One key feature of the response to compound events is the nonlinear nature of the response, meaning that the impact of compound events can be more intense compared to the cumulative impacts of individual events when occurring independently 54,55 . Here we formally address the nonlinearity in FT responses at the grid resolution, by inspecting the deviation from the superposition principle of linear systems. This is through simulating FT responses to historical baseline and six opposing individual and compound scenarios of change in temperature and snow depth. Table 1 Figure 5 summarizes the results of this analysis. The maps in the top row demonstrate the results related to deviation from the superposition principle at the grid scale under compound cooling (left column) and warming (right column). Under compound cooling, shades of blue and red show amplifying and dumping effects. For compound warming, blue and red colors represent dumping and amplifying effects, respectively. The white color highlights grids in which a linear response is observed. Black dots show the grids where the autocorrelation is present in temperature data and accordingly these grids are excluded from our assessment. The bar charts in the bottom row show the percentage of grids at each ecozone, where amplification in compound response is observed.
Considering this figure, although the amplified response of FD year to compound changes is observed across the province, it is mainly concentrated in the northern ecozones for both cooling and warming scenarios. Looking at both compound cooling and warming, more than 70% of the area in the three northern ecozones, i.e. Northern Arctic (EZ1), Southern Arctic (EZ2), Arctic Cordillera (EZ3) demonstrate an amplified response. While the expected values for amplifying effect, i.e. �FD year (−2, 10) − �FD year (−2, 0) − �FD year (0, 10) in EZ1, EZ2, and EZ3 are 5, 8, and 7 days, respectively, the amplified cooling effect can reach up to 16 days in some grids in northern ecozones. In contrast, the percentage of grids with amplified effect reduces to less than 10% in the EZ7-see Fig. 5c. Similar findings are obtained with regard to the compound warming scenario, although the amplification can also dominate the response of FD year in the south, e.g. in EZ8 (see Fig. 5d). Having said that, still more amplifications are observed in the three northern ecozones-see Fig. 5d. Although the expected values for amplified warming, i.e. �FD year (2, −10) − �FD year (2, 0) − �FD year (0, −10) are − 4, − 6, and − 6 days in EZ1, EZ2, and EZ3, the amplified response can get to − 11 days in a grid located in Southern Arctic. This finding has some important implications for the thawing landscape in the northern regions and a www.nature.com/scientificreports/ wide suite of environmental changes that can be initiated by the amplified response of FT to compound warming. On the one hand, an amplified response can facilitate the access to untapped natural resources of the north and may unleash an opportunity for northern agriculture. On the other hand, however, it intensifies permafrost degradation, which results in emissions of greenhouse gasses 19,20,45 . In contrast, the dumping effect dominates the FD year response in the southern ecozones, with the exception of EZ8 under compound warming. For both cooling and warming scenarios, the most severe dumping effect is observed in the Atlantic Maritime (EZ7), with the expected and extreme dumping of − 5 and − 9 days under compound cooling as well as 6 and 10 days under compound warming. The moderate amplification in FD year response in southern regions can be justified by the lower sensitivity of FD year to compound changes in temperature and snow depth shown in Fig. 3. Accordingly, we look at the grid-  . Lower sensitivity to snow depth condition in EZ8 was also shown in Fig. 4(f). Under compound changes, the shift in the mean of FTD DJF is approximately twice higher under (2, − 10) compared to (− 2, 10). The asymmetric response of FTD DJF to compound events is also manifested in the exceedance probabilities. For instance, in Montréal the exceedance probability of FTD DJF changes by 42% under (2, − 10), while the corresponding shift in the exceedance probability under (− 2, 10) is only − 31%. This is another line of evidence for complex responses of FT characteristics to compound cooling and warming. Comparing the expected compound responses with hypothetical linear responses under cooling and warming scenarios, an amplified response to compound cooling is observed in all three cities with expected values of − 2, − 3, and − 1 days in Montréal, Québec City, and Gatineau, respectively. Considering the compound warming scenario, the amplifying effect in Montréal and Gatineau are 2 and 1 days, respectively. It should be noted however that a dumping effect of − 1 day is observed in Québec City under the considered compound warming.
It should be noted that the nonlinearity in the response transcends to exceedance probabilities as well. In Montréal, for instance, the amplifying effects in exceedance probabilities under compound cooling and warming conditions are − 12% and + 3%, respectively. The amplifying response of FTD DJF to compound warming in a place like Montréal has an important relevance to the regional climate change, manifested by warmer temperatures and less snowfall 56 . Montréal marks one of the most populated regions in Canada with a high concentration of aging infrastructures, vulnerable to increasing FTD DJF . In Montréal, for instance, it is shown that under the historical condition, 216 km length of the city's watermains need to be replaced and another 2400 km are required www.nature.com/scientificreports/ to be repaired during the next 20 years 57 . This rate will significantly change under increasing FTD DJF . Current Findings in other cold regions indicate that more swings in FT cycles can translate to multiple million dollars for infrastructure maintenance and repair 47 . In addition, Marshall test of asphalt concrete shows reduced stability by up to 57% when the materials are exposed to only 6 more transient days 58 . This poses an addition vulnerability to aging infrastructures and the built environment 59 .

Summary and concluding remarks
Landscape FT is arguably the most important land-surface feature in cold regions, controlling physical, biological, and socio-economic processes along with their interactions with other elements of the environment. Historical characteristics of FT, however, are under significant alterations due to climate change, which has important implications for land and resource management in the north. Despite existing understandings of the impact of changing climate on FT cycles, current assessment frameworks are rather limited. On the one hand, in-situ data networks, often used as the main data support for understanding FT responses to changing climate conditions, are rather sparse and therefore fail to provide a synoptic view on different modes of landscape responses. On the other hand, the current projection paradigm based on the use of physically-based land-surface models is incomplete due to several assumptions and/or simplifications in the conceptualization, representation, and parametrization of the interacting processes that determine the state of FT at a given time and space. Our points of innovation are in the use of satellite remote-sensing data in conjunction with a formal statistical technique to overcome some of the above-mentioned limitations, at least at coarser spatial and temporal scales. First, we propose pairing gridded FT characteristics obtained from satellite retrievals with corresponding gridded data of temperature and snow depth, the two most influential hydroclimatic drivers of FT, to study the linkage between FT and climate characteristics. Second, we suggest C-vine copulas to develop a conditional model, with which the impacts of individual and compound changes in hydroclimatic drivers on FT characteristics can be quantified probabilistically.
We showcase the application of this framework in the province of Québec, Canada, a sub-country jurisdiction with an area comparable to Mongolia, the 18th largest country in the world. We show that the simulations obtained by parametrized C-vine copulas can capture the empirical moments of observed FT characteristics as well as the interdependencies between FT, temperature, and snow depth characteristics from 1998 to 2016. The application of the proposed framework for assessing the impact of individual and compound changes in temperature and snow depth reveals some important features of the FT response to climatic changes that have remained rather obscured. First, we highlight an ecozonal differences in the interdependencies between FT characteristics, temperature, and snow depth, pointing at the role of ecosystem conditions and/or latitudinal gradients in regulating the impact of changing climate conditions on FT characteristics. Moreover, a strong north-south divide in the FT response to changing climate conditions is observed. While the impacts of changing climate conditions are manifested in the annual extent of the frozen period in northern regions, they are revealed through alteration in the extent of the transient period during a typical winter season in the south. Through sampling the response surfaces related to FD year in Northern Arctic and FTD DJF in Mixed Wood Plains, we also demonstrate a different nature of the response to changing temperature and snow depth conditions between the north and the south. Our results show an asymmetric response to changing temperature and snow depth conditions despite differences in FT variables, regions, and/or spatial scales. In general, where there is a considerable sensitivity in FT response, the alteration due to warming temperature and/or thinning snow depth is more intense compared to corresponding scenarios with cooling temperature and/or thickening snow depth. Closer examination of the FT responses at the grid scale also revealed nonlinear, mainly amplifying, responses of FT characteristics to compound changes in climate conditions. This means that the response of FT to compound changes in temperature and snow depth is often more severe than the cumulative response of FT to changes in temperature and snow depth individually. These amplifying impacts can result in up to two weeks of alterations in the frozen period during a typical year in the north, and up to one week change in the transient period during a typical winter season in the three most populated regions of Québec, i.e., Montréal, Québec City, and Gatineau. While the amplified shortening in the annual frozen period in the north can result in unprecedented changes in the northern environment under compound warming, the similar climatic condition can result in additional stress to already aging infrastructures in the south. We show that the amplified response is not only observed in the expected values of the transient period in a typical winter, but also in the likelihood of exceeding the long-term historical values.
Our proposed methodology is generic and can be applied globally. We encourage inspecting the asymmetry, nonlinearity, and spatial variability of FT response to changing climate conditions in other Canadian regions and globally. In the context of Québec, landscape responses to freeze and thaw can have a wide range of implications from endangering already aging and vulnerable community infrastructures to substantial environmental changes due to amplified rates of permafrost degradation. The latter can be even globally relevant due to massive land-induced carbon emissions. Facing these challenges requires integrated and inclusive approaches that are supported by relevant and evidence-based scientific information. In that line, we will soon report the application of the proposed framework for estimating future FT characteristics in Québec using available downscaled climate projections: donc à bientôt! www.nature.com/scientificreports/ PM thawed), and inverse-transitional (AM thawed and PM frozen) 60 . We categorize transitional and inversetransitional states into one combined transient state that shows whether landscape switches between the frozen and thawed conditions in a diurnal cycle or not. Knowing the gridded daily states of FT, two critical FT characteristics, namely FD year and FTD DJF , can be extracted at the grid scale or each ecozone, and accordingly paired with corresponding gridded temperature and snow depth data. For air temperature data, we use the Global Meteorological Forcing Dataset (GMFD) provided by Princeton University available at https:// hydro logy. princ eton. edu/ data. pgf. php. GMFD dataset is constructed by blending the reanalysis data from the National Centers for Environmental Prediction, National Center for Atmospheric Research with a group of recent global observation-based data 37 . GMFD provides daily maximum and minimum air temperature at 0.25° × 0.25° for the period of 1948-2016. The daily mean temperature is calculated by averaging the daily maximum and minimum temperature at each grid. Monthly snow depth data are obtained from the Canadian Meteorological Center (CMC; https:// doi. org/ 10. 5067/ W9FOY WH0EQ Z3). CMC dataset is constructed by combining the information from in-situ snow depth measurements with optimal interpolation results of a simple physical snow accumulation and melt model 38 . The data is provided for the period of 1998-2020 at 24 × 24 km 2 across the northern hemisphere. As the grid size and centroid locations of the three considered data sources are not the same, we implement k-nearest neighbor interpolation 61 to re-grid the three data sets into a unique spatial scale and over the common period of 1998-2016. In brief, the k-nearest neighbor is a non-parametric approach to estimate a variable in a given point in time and space based on its neighboring values. The nearest neighbors are identified as those with the smallest Euclidian distance to the center of a reference grid, here the grid of FT, to which climate data are regridded. After finding the optimal nearest neighbors, a weighted averaging is applied to rescale the variables. The weight function gains its maximum value where the distance from the interpolated point is zero and decreases as the distance increases 62 . By implementing some numerical experiments, we find k = 4 as the optimal number of nearest neighbors to achieve the highest accuracy in modeling the mean and standard deviation of temperature and snow depth over different ecozones of Québec. The re-gridded monthly mean temperature and monthly mean snow depth are then matched with the corresponding FD year and FTD DJF at the same temporal and spatial scales. Before developing dependence models, we investigate the existence of autocorrelation in temperature, snow depth, and FT characteristics and exclude those grids in which the autocorrelation is significant. We find that this is the case in less than 17% of grids, and only for annual temperature. These grids are excluded from our analysis.

Methods
Proposed copula-based impact assessment framework. To model the joint dependencies between FT characteristics, air temperature, and snow depth, trivariate copulas can be used. Based on the Sklar's Theorem, the joint dependencies between FT characteristics ( FT ), mean temperature ( T ), and snow depth ( SD ) can be described as: where F 1 (FT) , F 2 (T) and F 3 (SD) are the Cumulative Distribution Functions (CDFs) for FT , T and SD , respectively and C is the trivariate copula function 63 . Among different alternative multivariate copulas, we employ canonical vine (C-vine) to represent joint distribution between the three above-mentioned variables in Eq. (1) 64 .
In brief, C-vine copulas decompose a high-dimensional joint distribution into a d(d − 1)/2 bivariate pairs of copulas arranged into (d − 1) trees; and accordingly, the joint distribution between FT , T and SD can be described as: where f (.) is the marginal PDFs and the conditional distributions that can be estimated as 65 : and where c(·) is the 3-dimensional copula density. As a result, the three dimensional joint density can be represented in terms of bivariate copulas as the following 66 : where c 2,1 (F 2 T , F 1 (FT)) and c 2,3 (F 2 T , F 3 (SD)) are simply written as c 2,1 and c 2,3 ; and the conditional pairwise copulas between F 1 (FT) and F 3 (SD) conditional to F 2 T , i.e. c 1,3|2 F 1 (FT), F 3 (SD)|F 2 T is shown by c 1,3|2 . In addition, c 2,1 , c 2,3 and c 1,3|2 are the densities of bivariate pairs. Having the C-vine copulas, the probability distribution of FT characteristics due to different quantitative changes in temperature and snow depth can be obtained through conditional modeling as: (1) F FT, T, SD = C F 1 (FT), F 2 T , F 3 (SD) (2) f FT, T, SD = f 2 T f 3|2 (SD|T)f 1|2,3 (FT|T, SD)  www.nature.com/scientificreports/ where F FT|T, SD is the conditional distribution function. Moreover, and Using Eqs. (7) and (8), Eq. (6) can be rewritten as: The estimated CDF of characteristics can be back transformed to the original quantile space using the inverse CDF function, assuming empirical distributions for FT , T and SD at each ecozone. The inverse form of h-function given in Eq. (9) is applied for this purpose. To extract the probability distribution of FT characteristics, a Monte Carlo-based simulation is adopted by generating 1000 random set of FT characteristics under known values of temperature and snow depth 67 . Given random uniform numbers of ε , given FT characteristics can be sampled as: Computer models for conducting these simulations are developed in CRAN R with the use of packages of VineCopula, CDVine, and copula [68][69][70] . Tree structures for C-vine copulas are selected based on the maximum spanning tree algorithm, in which copula parameters are chosen with respect to the interdependencies between pairwise variables 71 . A set of well-known parametric copula families (i.e. Frank, Gaussian, Student t, Clayton, Gumbel, and Joe) are used to develop, falsify and select pairwise bivariate copulas. The formulations of these copulas are provided in detail in other sources 72 . These copulas are parameterized using the Maximum log-Likelihood Method and considering empirical margins along with the Bayesian information criteria as the Goodness-of-Fit measure 73 . The pool of developed structures at each ecozone are then compared and evaluated based on their capability in representing the marginal FT characteristics and preserving empirical dependencies between a given FT characteristic and T or SD. Dependencies are quantified using the Kendall's tau nonparametric dependence measure and the associated hypothetical test 74 . The best C-vine copula structure is then used for conditioning the control of compounding changes in temperature and snow depth on FT characteristics. The one-way ANalysis Of VAriance (ANOVA) with Bonferroni correction is used to formally examine any change in the estimated dependencies across different spatial regions 75 .

Data availability
All data used in this study along with additional data related to elevation, land-use, and land-cover, as well as existing in-situ climatic and hydrometric networks are available through the Cold Region Data Accessibility Portal for Québec (CRDAP -QC; http:// wscc. encs. conco rdia. ca/ home. html).