Evaluating the Dependence between Temperature and Precipitation to Better Estimate the Risks of Concurrent Extreme Weather Events

School of Geography and Earth Sciences and McMaster Centre for Climate Change, McMaster University, 1280 Main St. West, Hamilton, ON L8S4K8, Canada Department of Earth and Atmospheric Sciences, University of Québec at Montréal (UQAM), Québec, Canada ESCER (Étude et Simulation du Climat à l’Échelle Régionale) Centre, University of Québec at Montréal (UQAM), Québec, Canada Department of Civil Engineering, McMaster University, 1280 Main St. West, Hamilton, ON L8S4K8, Canada


Introduction
Extreme weather events can have serious physical and economic impacts on urban and rural communities [1][2][3]. According to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change [4], warm temperature extremes are expected to increase and cold temperature extremes are projected to decline over the next several decades [5]. Additionally, significant increases in extreme precipitation events have been projected to occur in many areas over the 21st century but with strong variability across the years [6]. Estrella and Menzel [7] found that interdependence of climate variables may have a more severe influence on spatial climate change rather than the influence of any single climate variable, for example, the combined effects on both temperature and precipitation changes on drought occurrence and severity over the Canadian Prairies [8]. Other studies have also assessed the dependence between climate variables and extreme events. AghaKouchak et al. [9] found that the impacts of drought and heat events in the United States significantly increased when they occurred concurrently. Little et al. [10] found that sea-level rise and changes in the frequency and intensity of tropical cyclones will increase flooding risk in the future in the East Coast of the United States.
Most climate change impact studies consider variations in temperature and precipitation independently. However, these two variables are physically dependent through several mechanisms. For example, rainfall influences soil moisture, which in turn may have an impact on the surface and lowlevel air temperature through the effects of diabatic fluxes and partition between sensible and latent heat fluxes or lower/higher Bowen ratio [11,12]. e interaction of extreme temperature and precipitation may lead to high impact weather events and associated natural disasters with significant consequence on agriculture and other sectors of the economy [13]. For example, a drought accompanying heat wave can affect water availability for food production [14,15] as well as drinking water resources. In addition, temperature-precipitation interdependence in models may influence snow cover distribution and duration simulations [16,17], as well as flood occurrences and duration in spring over many regions particularly in Canada [18]. As a result, attention should be directed towards the combined effects of temperature and precipitation changes and their confounding impacts. Recent analyses of the dependence between temperature and precipitation are becoming a focus of meteorology and disaster prevention reduction research [19,20].
Many studies have shown that links between temperature and precipitation vary spatially and seasonally [21][22][23][24]. Johns et al. [25] used scatter diagrams of annual mean precipitation and temperature anomalies to show a global linear correlation which was simulated by most climate models for the period of [1980][1981][1982][1983][1984][1985][1986][1987][1988][1989][1990]. Based on the HadGEM2 model, Caesar and Lowe [26] analyzed the correlation between average annual temperature and extreme precipitation. ey found a proportionally high correlation between the two. Dai et al. [27] found a strong negative correlation between precipitation, maximum temperature (Tmax), and diurnal temperature range at short timescales during the warm season globally. A negative correlation between summer precipitation and temperature was found for most of the continental United States, which indicates that warm summers tend to be dryer and colder summers tend to be wetter [24]. In Canada, the temperature-precipitation interdependence tends to increase with latitude and is particularly strong in the Northwest Territories and relatively weak in the Prairies [21], but these correlation patterns vary between winter and summer months [28]. Southeastern Ontario and Quebec show almost no dependence of precipitation upon temperature during the summer (e.g., July), but there is a general negative correlation between monthly mean anomalies of temperature and precipitation for groups of months occurring over the period of May/June/July/ August/September [28]. In general, less/more precipitation falls during warm summer/winters months. In the eastern Rockies region, Isaac and Stuart [21] found that more precipitation occurs when the temperature is generally colder, regardless of season, as shown also in the study of Trenbeth and Shea [28]. is can be due to regular low-pressure systems (i.e., presence of the Alberta Clippers) moving from the west towards the Alberta area.
Most studies that analyze the dependence between temperature and precipitation have assumed a linear relationship and a normal distribution of both variables (e.g., [21,24,28]). erefore, the related results are likely to be inaccurate if one or both variables fail to satisfy the normal distribution condition, especially as precipitation is not normally distributed as well as low-level air temperatures for winter months due to nonlinear feedbacks induced by the presence of snow on the ground. Additionally, the linear correlation ignores high peak fluctuation and dependence structure [29]. In this study, we apply a copula-based approach to analyze the significance of the interdependence relationship between precipitation and temperature to overcome the issue of normality undistributed data. e interdependence of precipitation and temperature using copula has not been extensively investigated in the literature, especially for southern Ontario [11,30]. Cong and Brady [11] used copulas to study the dependence between temperature and rainfall for Scania (Sweden's southernmost province). ey found a significant negative correlation from April to September. In terms of copula, Student was selected as more appropriate for this region. ey stated that their results are strongly related to the studied region, and those results vary spatially (from region to region) and temporally. In this study, a copula-based approach was used to determine the inherent relationship between daily mean temperature and total precipitation (precipitation >0.1 mm). Five copula functions belonging to three different families were constructed to identify the joint distribution or interdependence of precipitation and temperature (see the following section). e proposed approach was applied to gridded (0.0833 degrees or ∼ 10 km × 10 km) climate data available over the whole Canada but was used for a specific region in southern Ontario, Canada. is study is the first that focuses on the examination and quantification of the nonlinear dependence between precipitation and temperature using copulas in southern Ontario.

Study Region and Datasets.
e study region is located in southern Ontario, Canada ( Figure 1). e study area is located between 42°18′N 83°01′W and 45°31′N 74°06′W and bounded on the south by the international border with the United States. e region encloses approximately 100,000 km 2 . is region is home to about 11 million people (2016 Census of Canada) and represents almost one-third of Canada's population and approximately 75% of the province of Ontario's population. e climate in eastern Canada, including the study area, is partly influenced by topographical effects, the Niagara Escarpment, and large lakes [31]. e Niagara Escarpment crosses the study area from Georgian Bay in the northwest to Lake Erie in the east and influences rainfall patterns in the region. In addition, the Laurentian Great Lakes have major effects on the climate in the region including the following: (1) the moderation of extreme temperatures in all seasons, (2) an increase in cloud cover and precipitation during winter (snow-lake effect), and (3) decreases in summertime convective clouds and rainfall events. ese effects are due to differences in the heat capacities of water and land surfaces and the rate of the Bowen ratio among the two surfaces and due to the large differences in moisture source especially from the Great Lakes [32,33]. e observed daily precipitation and maximum and minimum 2 m air temperature datasets used for this study were extracted from gridded historical weather station data (CanGrid) produced by Natural Resources Canada using interpolated Environment Canada's observed station across Canada. is dataset [34] covers Canada (south to 60°N) and contains daily variables from 1951 to 2013 on a Lambert conformal conic projection with 5′ arc minute spacing (equivalent to a resolution of roughly 10 km). ese climate data were generated using the "ANUSPLIN" software [35].
is dataset has been used in many climate change studies in the recent past (e.g., [36][37][38] and [8,39,40]). Canadian gridded data from 1699 CanGrid grid points located in southern Ontario corresponding to the boundaries of the study site over the period of 1951-2013 were used in this study ( Figure 1). Based on the grid-points information, the total annual precipitation for this period varies from the lowest recorded value of 725 mm to 1162 mm.

Copula Concepts.
Copula is a statistical notion used to describe the nonlinear dependence between random variables and to establish joint distribution of these variables using their marginal functions. It is also described as a function that connects univariate distributions to a multivariate distribution describing the dependence among correlated variables. e main advantages of using copula approach are as follows: (1) flexibility in choosing arbitrary marginal and structure of dependence, (2) extension to more than two variables, and (3) separated analysis of marginal distributions and dependence structure [41,42]. e joint distribution is fitted by estimating the marginal distributions of variables and their correspondence separately, which is not restricted to any parametric distribution (e.g., Gaussian distribution). Copula method has been widely used to examine interactions and associations between hydrological and/or climatological variables. De Michele and Salvadori [43] identified the relation between intensity and duration of storm rainfall in Italy by using Frank copula. Hao and AghaKouchak [44] used a set of parametric copulas to derive the Multivariate Standardized Drought Index from vectors of precipitation and soil moisture. Lee et al. [45] studied the influence of the tail shape of various copula functions (i.e., Gumbel, Frank, Clayton, and Gaussian) on drought bivariate frequency analysis. Renard and Lang [46] suggested applications of Gaussian copula on flood mitigation in France. Grimaldi and Serinaldi [47] have proved the adequacy of two copulas (Frank and Gumbel) on the flood characteristics analyzed for Kanawha River in West Virginia (United States). Chebana and Ouarda [48,49] presented regional multivariate flood analysis using copula and multivariate L-moments, as also used in the study of Saad et al. [18] who have developed trivariate copula for flood analysis over the Richelieu River located in southern Quebec. A Gumbel copula was used by Leonard, Metcalfe, and Lambert [50] to couple the seasonal rainfall maxima marginal distributions on the Murray-Darling Basin, Australia. Adlouni and Ouarda [51] proposed the application of copula to analyze the dependence of the water level of Saint-Louis Lake on the maxima flow on the Chateauguay River in Quebec (Canada). Rosa and Leite [52] presented that Frank and Clayton copulas fit well in studying a relationship between maximal flow and volume in Portugal. To establish the relation between the different flood characteristics, Salarpour et al. [53] applied the t-copula on the Johor River in Malaysia. A Gumbel copula was selected as the most appropriate model for trivariate frequency analyses of peak discharges, hydrograph volumes, and suspended sediment concentrations in Bezak et al.'s work [54]. Also copula was used to describe flood peak and volume, flood peak and duration, and flood volume and duration [55,56]; drought severity and duration [45]; and heat waves and drought [57]. erefore, copula function has been proved to a be very useful and effective tool for multivariate hydrological and climatological analysis and simulation. Mathematically, a copula is a multivariate probability distribution linking standard uniformly distributed marginals. Assuming that X and Y are pairs of random variables with cumulative distribution functions (CDF), . By the Sklar [58] theorem, the joint two-dimensional distribution function of X and Y, symbolized as H(x, y) with the cumulative joint probability p, can then be generated as follows: Here, C(u, v) is an arbitrary two-dimensional copula function.
e function C has the following elementary properties [59]: A variety of copula families have been described in the literature (e.g., [59,60]. Copula families differ in their parameter numbers and in their dependence structure, which have bearing on their complexity. e parameters of copulas control the strength of dependence. ese parameters are generally estimated using local optimization algorithms (e.g., [61][62][63], Bayesian (e.g., [64,65]), L-moments approach (e.g., [66]), and exact maximum likelihood (EML) methods [67]). In this study, five copula functions were used to describe interdependence between precipitation and temperature variables: Gaussian, Student, Frank, Clayton, and Gumbel (Table 1).
ese copulas belong to elliptical (Gaussian and t-copula) [60] and Archimedean (Frank, Clayton, and Gumbel copulas) families [59]. Note that elliptical copulas are often employed for simple dependency structure [68,69]. erefore, Archimedean family of copulas is more desirable for hydrological analyses because it can be easily constructed and can be applied whether the correlation between the variables is positive or negative [56]. ese copulas were selected for analysis because they have been commonly used to evaluate climate variable interdependence in hydrological and climatological studies.

Analysis Procedure.
Available climate data were first treated with the following criteria: days with no precipitation (<0.10 mm) were removed from analysis and all other days were grouped into seasons (December-January-February (DJF), March-April-May (MAM), Jun-July-August (JJA), and September-October-November (SON)). To fit the joint distribution between temperature and precipitation using copula approach, the following procedure was used: (1) During the first step of analysis, adequate marginal distribution (probability distribution function, PDF) is chosen for each of the studied variables. In fact, the marginal distribution of daily average temperature and daily total precipitation in the study region was identified. Identifying the fitted probability distribution allows predicting the probability of exceedance for a specified magnitude (quantile) or the magnitude associated with a specific exceedance probability. ere is no theoretical basis for the choice of probability distribution and the parameter estimation method. In this study, twelve probability distributions were fitted for each grid and season. e fitted PDFs are Beta, Exponential, Extreme value, Gamma, Generalized extreme value, Generalized Pareto, Inverse Gaussian, Logistic, Log-logistic, Log-normal, Normal, and Weibull (see Table 2). Several methods to estimate the parameters associated with these distributions are presented in the literature, that is, the maximum likelihood method (e.g., [70,71]), method of moments (e.g., [72]), and L-moment method (e.g., [73]). In this study, the maximum likelihood method was used. e selection of the best fit distribution for each variable (i.e., best marginal distribution fit of the grid data) is based on the Bayesian information criterion (BIC) proposed by Schwartz (1978). e smallest BIC values identify the best fit distribution [74].
(2) In the second step of the analysis, the cumulative distribution function (CDF) of these distributions is used to compute the marginal cumulative probabilities . is computation assumes that F i and G i are the selected PDFs for total precipitation and daily average temperature for the grid i, respectively. Note that u i and v i are strictly increasing functions and range within the interval [0, 1]. (3) In the third step of the analysis, five copula functions (as defined in Table 1) are fitted to the marginal cumulative probabilities of daily average temperature and total precipitation. To select the suitability (best fit) copula function, a method based on the comparison of theoretical and empirical copula functions was used. In detail, for each copula function, the root mean square error defined as the average square distance between theoretical and empirical copulas was calculated. e appropriate copula of each grid is the one that has the lowest root mean square error. Note that the parameters of the copula were estimated using the exact maximum likelihood (EML) method [67].
Once the joint distribution of each grid is established (including the margins, copula, and their parameters), it can be used, for example, to calculate the joint risk (probability). Joint probability of precipitation and mean temperature is very important for the management and assessment of the risks imposed by extreme meteorological and hydrological events. It helps in the management of resources such as agriculture sector productivity. Generally, in the bivariate case, four simultaneous events can be of interest: (1) simultaneous exceedance {P ≥ p, T ≥ t}, (2) exceedance-nonexceedance {P ≥ p, T≤ t}, (3) non-exceedance-exceedance {P ≤ p, T ≥ t}, and (4) simultaneous nonexceedance {P ≤ p, T ≤ t}. In this study, we mainly focused on computing the probability of concurrent extreme temperature and precipitation events including wet/cool (denoted by P w/c ) in the winter (DJF) and dry/hot (denoted by P d/h ) in the summer (JJA). Following Zhou and Liu [75], in order to capture a large number of events, we used the 25 th and 75 th quantile thresholds of T and P to define these two probabilities: where p 25 , p 75 , t 25 and t 75 are, respectively, the 25 th and 75 th quantile thresholds of precipitation and temperature. To quantify the usefulness of introducing the dependence between precipitation and temperature, these two probabilities were computed for each grid by assuming that P and T are independent (denoted as (P w/c ) I and (P d/h ) I ) and dependent (denoted as (P w/c ) D and (P d/h ) D ) variables. erefore, using some probability manipulation, they can be defined as

Name Function
Beta Exponential where µ is the mean.
Extreme value where µ and σ are, respectively, the location and scale parameters. Gamma where a and b are the shape and scale parameters and Γ(·) is the Gamma function.
Generalized extreme value where µ, σ, and k are, respectively, the location, scale, and shape parameters.
Generalized Pareto where k, σ, and θ are, respectively, the shape, scale, and threshold parameters.
, for x > 0, where µ and λ are, respectively, the mean and shape parameters.
x ≥ 0 where µ and σ are, respectively, the mean and scale parameters.
, for x > 0 where µ and σ are, respectively, the mean and scale parameters.

Normal
If X follows the log-normal distribution with parameters µ and σ, then log(X) follows the normal distribution with mean µ and standard deviation σ. Weibull where σ and k are, respectively, the scale and shape parameters.
x, y ∈ R, ρ is the linear correlation, and Θ is the normal standardized function.
dxdy; x, y ∈ R, ρ is the linear correlation coefficient, and t − 1 is the inverse cumulative distribution function (CDF) for Student's t-distribution with κ degree of freedom.
where θ measures the dependency between u and v. where F and G are the selected PDFs for daily total precipitation and daily average temperature and C represents the best fitted copulas for the selected grid. After computing these probabilities for each grid of the study region, two  these concurrent extremes by ignoring the dependence between temperature and precipitation. More details are presented in the Results section.

Results
e spatial variability of monthly mean temperature that was averaged over the study period for all four seasons is shown in Figure 2(a). A high variability in temperatures between seasons demonstrates the necessity of seasonal-scale copula analyses for this study to capture variability and extreme events. Monthly average temperature varied from − 20 to 3°C in the winter (DJF), from − 10 to 19°C in the spring (MAM), from 13 to 26°C in the summer (JJA), and from − 5 to 20°C in the fall (SON). During all seasons, the highest temperatures were found in southwestern area (urban centers) close to the Lake Ontario and Lake Erie in the western portion of the study area. e lowest temperature was found in the northern (rural) areas and at high elevations in the northwest portion of the study area. Overall, as spatial variability, we observed approximately +5°C gradient from the southwest to the north across the study region. e spatial variability of total precipitation by season and temporally averaged over the studied period is shown in Figure 2(b). Total monthly precipitation varied from 178 to 285 mm for DJF, from 185 to 240 mm for MAM, from 195 to 265 for JJA, and from 208 to 315 mm for SON. Variability in seasonal total precipitation was spatially more complex than average temperature, reflecting combined effects in surface condition and topography. e lowest precipitation occurred in the eastern and southern parts of the study region close to the Lake Ontario and Erie Lakes but with a maximum in precipitation in the western part close to the Huron Lake and Georgian Bay from fall to spring months (i.e., snow-lake effect). Day-and-night heating/cooling and wind patterns adjacent to the lakes were more variable compared to areas located further inland. Increased variability in these on-shore winds may have contributed to lower/higher precipitation adjacent to the lakes in southernmost/westernmost Ontario. e area of highest monthly precipitation was observed in the western and northern parts of the study region east of Georgian Bay and Lake Huron is also a high elevation sector relative to the southern part of the study site and higher forest cover. Weather patterns in this area were also more influenced by precipitation or known snow bands extending eastwards of Lake Huron and Georgian Bay in the west.
ese regional influences may have contributed to higher precipitation in this part of the study area. In general, a high seasonal variability (standard deviation not shown) was found during the fall and winter seasons when compared to variability during spring and summer. Figure 3(a) shows the spatial variability of Spearman rank correlation between daily data with total precipitation and daily average temperature for all four seasons.
e p values of corresponding significance test are shown in Figure 3(b). e two main advantages of the Spearman rank correlation test as opposed to simple linear correlation test are the following: (1) it can be computed without any assumption about the normal distribution of data and (2) it has low sensitivity for inhomogeneous time series. is correlation has been used for trend analyses in climatologic and hydrologic time series [76]. A positive rank correlation was found between daily total precipitation and average temperature across the studied area and for all four seasons. Spatially, correlation was in most significant grid point within the entire study area. e highest correlation was observed in western parts in the spring and autumn. Also a high correlation was observed in the central parts of studied area adjacent to the Erie Lakes and Huron in the summer but with low correlation values outside of this area.
is central and southern region was proportionally warmer and fog occurrence was higher as compared to elsewhere in the study area. Isaac and Stuart [21] show that, on 90% of days during which fog occurred, precipitation was also reported. erefore, the portion of the study site with common fog days was expected to have a large number of days with precipitation. In general, as suggested by Isaac and Stuart [21] from a station-scale study, more precipitation falls during warm weather conditions during the winter over southern Ontario and Quebec. Figure 4(a) shows the best fitted distribution (PDF) for seasonal daily average temperature from 1951 to 2013 for the study area. To show the spatial and temporal change of PDF that may exist at the extreme ends of the study region, two grids were selected: one in the furthest north east (grid 1) and the other in the furthest southwest (grid 1699) (Figure 4(b)). Selected distribution was found to vary spatially and seasonally between the two ends of the study area. During the winter season, based on average temperature, the studied area was grouped into two different regions based on temporal data distribution: (1) the southwestern part that followed a normal distribution curve and (2) the northeastern part that followed a GEV distribution curve with a high negative shape parameter (GEV type III). e shape parameter belonged to a narrow range, approximately from − 8 to − 10. e extreme cold temperature in the northeastern part of the study region during winter can explain the use of GEV III (Figure 4(b), red curve). GEV with a positive shape parameter (GEV type II) was found to have the best GEV distribution for spring and fall seasons. e shape parameter varied spatially between 2 and 5 during spring and between 5.5 and 9.5 during fall. e observed variation in shape parameters can be explained by the large temporal variability of daily average temperature that occurred during these two seasons, with nonlinear feedbacks according to the snow cover and frozen conditions of the ground, which affects diabatic fluxes and underlying heating/cooling process (Figure 4(b)). e positive shape parameter value implies a heavy tail distribution. e spring and fall data distributions show that extreme values occurred more frequently relative to other values. A two-parameter Weibull distribution with large positive shape was observed for summer data (Figure 4(b)). Figure 5 shows the best fitted distribution for seasonal total precipitation from 1951 to 2013 for the study area (5A) and the PDFs of daily total precipitation at grids 1 and 1699 for different seasons (5B). e best fit distribution was Advances in Meteorology observed to change spatially and seasonally with high spatial variability during spring and summer. is was an expected result due to the high spatial variability of precipitation processes that occurred during these seasons. Log-normal (LN) and Generalized Pareto (GP) were found to be the best distributions to fit to the daily data during all four seasons. LN distribution assumes that the logarithm of daily total precipitation is normally distributed and is useful when the variable of interest is skewed to the right. LN was found to be the best distribution in grids to fit with proportionally low precipitation amount. Additionally, GP was found to be the best distribution to fit data in grids with large daily precipitation amount. Figure 6 shows the mean square error (MSE) between theoretical and empirical copula functions of all 1699 grids. For each grid, the best copula function that fits the joint      distribution of precipitation and temperature is the one that minimizes MSE. In general, no copula was the best fit for all seasons and/or grids. at was expected due to differences in geographical and geophysical conditions across the region which influence temperature and precipitation variations. Less variability in terms of best copula was observed during winter and summer as compared to variability of best copula during spring and fall. Copula variation can be explained by the fact that temperature has the same positivity sign for all grids during winter and summer (i.e., positive for summer and negative for winter) which is not the case during spring and fall (whereas the temperature positivity sign depends upon the selected grid). Gumbel was found to be the best copula during winter and Clayton copula was the best performer during summer.
ese two copulas are asymmetric and belong to the Archimedean family. Gumbel copula is used for modelling heavy dependencies in right tail, which means that more precipitations in winter are generally related to warm weather in this region. In fact, warmer air can contain more water vapor than cooler air, which can result in higher probability of more intense precipitation (i.e., Clausius-Clapeyron relationship). erefore, Clayton copula is used for modelling heavy dependencies in left tail, explaining the dependency between low amount precipitation and extreme temperature in summer. Due to the mechanisms mentioned previously in this section, greater variability in temperature and precipitation during the spring and fall seasons leads to more variability in terms of the best copulas during these two seasons.
In this section, an application showing the effect of ignoring the dependence between temperature and precipitation on the estimation of risks of extreme events such as wet/cool in the winter and dry/hot in the summer is presented. Figure 7 shows the differences between likelihoods of concurrent extremes computed by assuming P and T as dependent and independent variables. We found that generally the differences are positive (varying between 8 and 14%), which means that ignoring the dependence leads to underestimating the probabilities of occurrence of wet/cool and dry/hot events up to 14% in southern Ontario. Also, we observed that likelihood of differences is larger in the case of wet/cool event compared to that in the case of dry/hot event.
is can be explained by the fact that precipitation and temperature are more correlated in the winter than in the summer for this region. In particular, the differences are more significant in the urban centers close to the Lake Ontario. In conclusion, this analysis indicates that the correlations between P and T have a direct effect on the estimation of occurrence risk of concurrent climate events.
Another example of the usefulness of modelling the risk of an extreme event based on the joint distribution of temperature and precipitation, as opposed to basing the model on these two variables separately, is shown in Figure 8. In fact, this figure shows the temperatureprecipitation space of grid 1682 (i.e., Toronto area) for the four seasons. e contours represent the bivariate quantile curve for different, simultaneous, nonexceedance events. Note that the quantile function expresses the magnitude of the event in terms of its exceedance or nonexceedance probability, which is also related to the return periods. e quantile curve is composed of two parts: the naïve part (tail) and the proper part (central). e naïve part is composed of two segments starting at the end of both extremities of the proper part. Detailed description and proprieties of the bivariate quantile function can be found in Chebana and Ouarda [77]. For the grid 1682 example in winter, for risk value of p � 0.9, the univariate quantile values of precipitation and temperature are 10 mm and 5°C, respectively. e violet-colored curve in Figure 8 represents the bivariate quantile extracted from the joint distribution. Note that the combination of the univariate values (precipitation � 10 mm and temperature � 5°C) does not belong to the bivariate quantile curve. is combination corresponds to another risk (p) smaller than the actual risk of p � 0.9. erefore, the wrong conclusion (in terms of magnitude and return period) could be made without taking into consideration the joint or confounding distribution. e generality of the bivariate quantile function, which can give several possible scenarios related to the same risk p, is not the case of univariate quantile. e univariate quantile represents the extreme points of the proper part of the bivariate quantile curve indicated as the black-dashed line in Figure 8.

Discussion
In this study, a copula-based approach was used to determine the inherent seasonal relationship between the average temperature and the total, nonzero precipitation at daily resolution across southern Ontario. We found a positive correlation between temperature and precipitation in the entire region. Our analytical results are consistent with previous studies reported in the literature, which indicate a positive correlation between precipitation and temperature in this region (e.g., [21,78]) but with a different conclusion as in the work of Trenberth and Shea [28] in which they combine all "warmer" months from May to September. Isaac and Stuart [21] computed a temperature-precipitation index that derives the percentage of precipitation based on temperatures colder than median daily temperature for 51 Canadian meteorological stations. ey report a positive correlation during winter for all stations located in southern Ontario and Quebec. e results of this study are an important step forward to characterize and quantify nonlinear dependence between precipitation and temperature in an area that is prone to the occurrence of extreme weather events and has profound social and economic significance for Canada. e proposed approach demonstrated here is flexible and assumption-free. In fact, the modelled dependence structure does not require the normality of marginal distributions of variables which provide flexibility in correlating variables and also prevents the necessity of similar marginal distributions.
We found that the Gumbel and Clayton copulas are the most suitable to fit the dependence between temperature and precipitation during winter and summer, respectively. Gumbel and Clayton are two asymmetric copulas, with right (Gumbel) and left (Clayton) tail dependence [79]. is explains the significant correlation between extreme temperature and extreme precipitation in this region during these seasons. Mechanistically, this reveals that less precipitation falls in general during cold weather conditions in winter, and more precipitation falls during warm weather in fall. Asymmetric copulas perform better to identify nonlinear correlation between temperature and precipitation in this region.

Conclusion
In this paper, a copula-based approach was presented to model the seasonal joint or confounding distribution of temperature and precipitation. e proposed approach is flexible and free of assumptions. Five copula members belonging to three families were fitted to 1699 grids (∼10 km × 10 km resolution) in southern Ontario. For each grid, an information criterion was computed using empirical and theoretical copulas and the best fit copula was selected based on the strength of nonlinear correlations. Modelling of the joint distribution of precipitation and temperature will help to produce improved simulations of weather events which may help to increase the accuracy of risk evaluations.
Results showed that no copula performed consistently as the best copula for all four seasons and for the entire region. Gumbel performed as the best fit copula for winter and Clayton copula performed as the best fit copula for summer. More variability in terms of best copula was found in spring and fall, which may be due to the variation in temperature during specific events or thresholds (i.e., around 0°C) between grid points and mixing of precipitation types which can lead to different links with temperature (as noted in [28]). By extracting the multivariate and the univariate quantiles related to a preselected risk, it was found that ignoring the joint or confounding/combined distribution of precipitation and temperature may lead to underestimation of the risk of an extreme event. is underestimation may lead to a misinterpretation and a wrong conclusion in terms of return periods. Our study also reveals that the relationships and physics of combined occurrence of precipitation and temperature events should be taken into account in interpreting climate changes and in climate risk analyses. Extremes in meteorological variables may have a significant impact on ecosystems and society through the occurrence of extreme weather events. Evaluating the risk potential of climate variable extremes is critical for resilience policy and mitigation of the negative effects of climatic change. A potentially valuable extension of this research is trivariate copula analyses to connect precipitation and temperature with crop production planning and agricultural economics and flood evaluation. Such study could be used in developing risk reduction strategies for farmers and decision-makers, which will become increasingly important in the face of climate change and its associated modification in extremes, water cycle, and hydrometeorological hazards.

Data Availability
e daily CanGrid data used in this study are obtained from the following website: ftp://ftp.nrcan.gc.ca/pub/outgoing/ canada_daily_grids.

Conflicts of Interest
e authors declare that they have no conflicts of interest.