Stochastic modeling of spatial dependency structures of extreme precipitation in the Northern Great Plains using max-stable processes

Extreme events can be defined as events that are sudden, unprecedented, and unexpected. In most cases, these events have both spatial and temporal characteristics (Zhang et al. 2001). For instance, extreme precipitation can lead to flash floods with associated devastation, and when this occurs it has a significant spatial extent. Also, from a temporal characteristic's perspective, each event may occur once in a decade, once every two decades, and so on (Villarini et al. 2011; De Paola et al. 2018; Gründemann et al. 2023); therefore, it is crucial to adequately quantify these events, especially in infrastructural design. Water resource engineers or hydrologists are interested in the magnitude, the interval between the events and the spatial structures of extreme events. Statistical hydrologists are interested in metrics such as 25-, 50- ,and 100-year flood levels or return periods (i.e., at least 1 in 100 years events) for design and management purposes. Furthermore, the modeling of precipitation patterns is important for water security and sustainable agriculture. The impacts of extreme precipitation due to a combination of excess snowmelt and rainfall in North America have been well documented (Curriero et al. 2001; Kunkel 2003; Novotny & Stefan 2007; MI 2023). For instance, in 2014, the Assiniboine River was flooded due to excess precipitation affecting both Manitoba and Saskatchewan in Canada. Ahmari et al. (2016) reported that the damages in Manitoba and Saskatchewan were more than CAD 200 million and CAD 360 million, respectively. Szeto et al. (2015) further reported that 2014 flooding in the southeast Prairies may have been influenced by anthropogenic factors (such as landscape modifications) may have also contributed to the increased runoff.

One of the largest river basins in North America is the Nelson Churchill River Basin (NCRB). This basin covers four states in the United States and four provinces in Canada. Westra et al. (2013) conducted a study on nonstationary generalized extreme value (GEV) analysis for each of the 8,326 stations across the entire globe. The authors concluded that it is difficult to visually discern any clear spatial pattern in the daily annual extreme precipitation in North America. Similarly, Zhang et al. (2019) concluded that there appear to be no detectable trends in extreme precipitation in Canada. Zhang et al. (2019) further stated that more weather stations had experienced an increase than a decrease in extreme daily precipitation. Therefore, there is a need to do a regional trend analysis to detect these changes, as the NCRB is an important basin for Canada's water resource supply. Furthermore, the southern Prairie of Canada, which is an intensive agricultural region, covers a substantial part of the NCRB. Considering this, extreme events due to excess precipitation can be catastrophic to Canada's food security and water management. For the United States, the NCRB covers important agricultural states such as North Dakota and Minnesota. The NCRB is also crucial to the US society and economy.

Mathematical models developed through geostatistic techniques are insufficient in resolving the spatial characteristic of extremes (Ribatet 2017). In other words, these geostatistic models focus on mean process levels (Padoan et al. 2010). Extreme value theory (EVT) is a statistical method commonly used to study extreme-event behavior in the univariate case. EVT depends on asymptotic theory to model the tail, but extremes are spatially dependent; therefore, multivariate series of rainfall intensities from several locations were modeled by means of max-stable processes (MSP, Coles & Tawn, 1991; Coles, 1993; Schlather & Tawn, 2003; Ribatet 2013). In other words, MSP has the same conceptualization to the univariate GEV distribution (Padoan et al. 2010). MSP can, therefore, be defined as an adapted EVT that can generalize extremal dependence from a continuous space standpoint. The Spatial GEV model provides a pointwise estimate of the return period or return level, but it does not consider spatial dependence. In other words, it makes a wrong assumption of independence of extreme weather stations as, in reality, closer weather stations will have similar values. When defining the spatial GEV model parameters, simple linear functions of geographical coordinates and the topography of the weather stations can be used. Other environmental covariates such as the minimum/maximum temperature, average annual soil moisture index, the normalized difference vegetation index (NDVI), and so on can be used as additional variables that can help resolve the spatial variability of the extremal dependence. It is also possible to have more than one competing spatial GEV model. Takeuchi information criterion (TIC) computes TIC values that can be used to select the best model among competing spatial GEV models (Sakamoto et al. 1986; Varin & Vidoni 2005; Gao & Song 2009). This TIC estimate is equivalent to the Akaike information criterion (AIC) when the model is misspecified (analysis of variance (ANOVA)) or the Fisher Test can also be used when evaluating if any additional covariates are necessary to be included in the model configuration.

Starting from the work of de Haan (1984), MSP has been applied in various applications dealing with spatial extremes (Smith & Stephenson 2009; Ribatet 2017). Davison et al. (2012) described various spatial extremes based on latent variables, copulas, and spatial max-stable processes using an extreme rainfall dataset from Switzerland as a case study. Their study concluded that the model based on the variable modeling provided a better fit to marginal distributions but underperformed in fitting the joint distributions of the extremes; therefore, copula or max-stable models were recommended as the most plausible models. Furthermore, several parametric max-stable models have been used to model specific stochastic spatial extremes. Ribatet (2013) demonstrated the application of MSP in the modeling of extreme wind gusts in the Netherlands using ‘SpatialExtremes’ in R statistical package (Ribatet 2013). The authors found that the Extremal-t model appears to be the most competitive models. Padoan et al. (2010) used MSP to analyze US precipitation data using 45 weather stations. The authors concluded that the likelihood-based inferential approach was well suited for the joint modeling of marginal and dependence parameters. In addition, the computation was at a moderate cost and was a good estimator for a limited sample. Reich & Shaby (2012) used MSP to model yearly maximum extreme precipitation from a regional climate model. The authors concluded that there would be a statistically significant increase in the upper quantiles of the precipitation under future climate scenarios. Their proposed model is also closely related to a Gaussian extremal value process. Westra & Sisson (2011) used the MSP process to analyze 30 sub-daily gauges in east Australia from 1965 to 2005. The authors found that MSP improves the precision of inference and an increase of 18% for the 6-min rainfall was noted. Since the focus of this study is not to provide a comprehensive review of MSP applications, other applications can be found in Erhardt & Smith (2012), Gaume et al. (2013), Ribatet (2013), etc. To develop a reliable MSP model, reliable and complete time series datasets are needed. One major challenge facing this kind of study is a lack of good-quality and complete data, especially in a sparse weather network area. Simulations from the fitted MSP provide a solution to this challenge because realizations can be generated from the MSP model after the diagnostic process. Different types of MSP model families have been developed over the years. Starting from the Smith model (Smith 1990), also called the Gaussian extreme value process (Schlather 2002), there have also been the extremal Gaussian process (Schlather 2002), Brown–Resnick process (Brown & Resnick 1977), and Extremal-t process (Nikoloulopoulos et al. 2009). More information on each of these MSP model families can be found in Ribatet et al. (2016).

After an exhaustive literature search, to the best of the author's knowledge, this study represents the first attempt to characterize the spatial dependency and trend of annual daily maxima precipitation in the NCRB. Therefore, this paper aims to develop a stochastic model based on the MSP that will capture the spatial dependency structure of the daily annual extreme (total) precipitation in the NCRB from 1970 to 2020 using the daily annual maximum precipitation datasets extracted from Environment and Climate Change Canada and Global Historical Climatology Network daily repositories. Furthermore, this study will provide new insight into the impacts of both spatial and temporal covariates in the extreme spatial modeling of annual maxima precipitation in the NCRB. In other words, the temporal covariates can help to quantify the impact of climate change on extreme precipitation. This study is crucial to understanding extreme spatial precipitation for flood protection and preparedness, water management, and hydrologic regimes of a cold river basin such as the NCRB.

The NCRB is a basin inside the North Great Plains (NGP) (Figure 1). The NCRB covers four states and four provinces in the United States and Canada. Vegetation is generally composed of tall grass, mixed, and short grass, and the NGP is regarded as one of only four remaining intact temperate grasslands globally (WWF 2023). It is also the agricultural heartland of Canada, extending from the Precambrian Shield in Winnipeg to the Rocky Mountains in the west (U.M. 2023). This region experiences a heterogeneous climate, climatic extremes such as droughts, and severe storms (NASA 2021). On the Canadian side, the three climatic regions have been identified as the Prairies, Northwestern Forest, and Northeastern Forest. Westmacott & Burn (1997) further reported that the ‘Prairies, Northwestern Forest, and Northeastern Forest climatic regions have statistically significant increasing temperature changes corresponding to 0.9 °C, 1.3 °C, and 0.5 °C over the last century’ (Gullett & Skinner, 1992). This implies that these regions have been impacted by climate change, the consequences of which could be long dry spells or frequent extreme precipitation. Under average conditions, the western portion of the NCRB is generally wetter compared to the Prairie portion. The northern portion of the basin is generally colder, with an average temperature below −5 °C, whereas the southern portion can have an average temperature of more than +5 °C. A heterogeneous rainfall pattern characterizes the NCRB. This is due to the enormous quantities of orographic rainfall in the west and the convective precipitation in the Prairies (which is dominated by short-duration, sporadic, extreme rain), causing millions of dollars in damages. From an average precipitation pattern perspective, the Rocky Mountain areas (western portion) show high average precipitation, whereas the eastern portion of NCRB (southern Prairies: Alberta, Saskatchewan, and Manitoba) is generally dry with an average annual precipitation of approximately 320 mm (Bajracharya et al. 2020; Boluwade 2023).

From a hydrologic standpoint, the NCRB drains most of the northern United States and Central Canada (Boluwade et al. 2018). It is more than 1.4 million km². In addition, the NCRB has eight sub-basins: the Saskatchewan River Basin (SRB, ∼335,900 km²), Lower Nelson River Basin (LNRB, ∼90,500 km²), Assiniboine River Basin (ARB, ∼182,000 km²), Churchill River Basin (CRB, ∼281,300 km²), Red River Basin (RRB, ∼100,000 km²), Lake Winnipeg River Basin (LWRB, (∼190,000 km²), and Winnipeg River Basin (WRB, ∼106,500 km²). Several of these sub-basins have experienced significant extreme flooding due to excess rainfall, especially during the spring and summer periods. For the RRB and ARB, the highest flooding was recorded in 1950, 1997, 2009, and 2011. For instance, in 1997, the RRB was recorded as having the most severe flooding since 1852 (Rannie 2016; MI, 2023). The 1997 flood event started with a dry summer in 1996. This was followed by heavy rainfall in the fall that saturated the soil. There was heavy snowfall the following winter. MI (2023) further reported that the total precipitation from the beginning of winter to early May at the crest of Red River was around 221 mm. This is well above the average norm of 130 mm. The resulting flood affected all the areas in the Red River Valley covering the cities of Winnipeg in Canada and Fargo in the United States. Around 27,000 people were relocated in Manitoba with associated damages close to $1 billion. Similarly, the 2011 flooding due to extreme precipitation started in October 2010. Around this period, southern Manitoba experienced one of its wettest years. This was due to excess snowfall and rainfall of up to 50–100 mm. At soil freeze-up, the soil moisture was already at full saturation, meaning there was no space for additional moisture infiltration to the soil. Furthermore, a delayed snowmelt in the spring compounded this phenomenon. In other words, there was still a significant snowpack on the land surface by mid-April. When spring precipitation arrived, it exacerbated the flooding problem, leading to a significant increase in heights in several rivers in the NCRB such as the Assiniboine River, Red River, Souris River, Qu'Appelle River, and Pembina River. From the foregoing, total precipitation (snow water equivalent plus rainfall) was considered in this study because of the dynamics in the hydrologic processes and interconnectedness between rainfall in the summer/spring and snowfall in the winter. The connection between these two processes was responsible for the extreme precipitation event recorded in this basin. For instance, heavy rainfall in the fall (which saturates the soil profile) usually precedes heavy snowfall in the winter and then is compounded by heavy rainfall in the spring.

This study extracted daily maximum precipitation from two main repositories. The Global Historical Climatology Network daily (GHCNd) provides daily climate summaries from weather stations across the globe. These datasets come from several sources and have been quality controlled. The database contains information from more than 100,000 weather stations across 180 countries and territories (NOAA 2023). Variables such as daily temperature, min–max temperature, total precipitation, and snow depths are provided by the database. For this study, daily annual maximum total precipitation (rainfall plus snow water equivalent) was extracted from 67 weather stations located within the NCRB. An R package (FedData, Bocinsky 2022) code (https://docs.ropensci.org/FedData/) was used to extract this dataset from the repository. These are the stations that passed the quality control (i.e., without missing annual values). The other database used was the Historical Climate Data from the Environment and Climate Change Canada (ECCC) website (https://climate.weather.gc.ca). An R package (weathercan:LaZerte & Alber 2018) were used to extract total precipitation from an additional 13 weather stations, mostly in the Canadian portion of the NCRB. Total daily extreme precipitation information from January to December from 1970 to 2020 is considered.

Topographic features of landscapes such as elevation values are correlated with orographic precipitation. The NCRB has the Rocky Mountains at the western side. It is important to consider the impact of this variable as a driving factor for the extreme annual precipitation in the basin. The topographic data were extracted from the GTOPO30 elevation dataset (http://eros.usgs.gov/products/elevation/gtopo30.html) developed by the United States Geologic Survey (USGS). Using the extract function in R, corresponding values of the USGS elevation data were extracted using both GCHNd and ECCC weather station geographical coordinates.

The block maxima (B.M.) approach was used to select the extreme value from the original data. One year was selected as the block B.M., and maximum values were extracted for each year (from 1970 to 2020). Spatial dependence was tested using the extremal coefficient for each pair of stations used in the study.

The best model would be selected based on Takeuchi's information criterion (TIC) (also called the composite likelihood information criterion (CLIC)). TIC is defined as ‘the linearization of maximum likelihood estimator bias which shrinks the model parameters towards the maximum entropy distribution, even when the model is misspecified’ (Dixon & Ward 2018)

The corresponding standard errors for the location, scale, and shape coefficients are: ⁠, ⁠, {0.03511}, respectively.

The corresponding standard errors for the location, scale, and shape coefficients are: ⁠, {0.03434}. The standard error for the time variable is 0.0007569. According to Padoan et al. (2010), statistical analysis tests such as ANOVA can be used to evaluate whether time as a covariate is necessary to be included in the model, as stated above. ANOVA analysis for the best model (TRS₁₀) with and without time shows the p-value of the test is 0.6559. This indicates that the best model with time (as an additional variable) is not significantly different from the one without time at the level of α (significance level) = 0.05. In other words, the ANOVA analysis shows that time is not significant as a predictor. This may be due to both downtrend and uptrend seen in the western and southern sections of the basin. However, the Mann–Kandall test revealed that there is a temporal trend in the extreme data; therefore, time would be included in the spatial GEV model parameters definition as specified above. Furthermore, the TIC of the TRS₁₀ with and without time is 16,760.72 and 16,753.37, respectively. Further analysis and discussion in the following sections will show the impact of including time as a temporal covariate.

The implication of the above is that there is a reduction of 0.0000206 mm of extreme precipitation per year across the entire NCRB. Again, the relatively small precipitation rate may be due to both a reduction and increase in the various sections of the basin as shown in Figure 3. This overall relatively small reduction does not negate the fact that the southern portion of the basin has shown an increase in extreme trends, which is consistent with other findings (C2CES, 2023) and the western portion of the basin shows a large negative deviation from the mean of the annual extreme precipitation. Infact, this relatively small reduction is consistent with Zhang et al. (2019)’s study. The estimated values for the dependence parameters of the Extremal-t MSP model for range, smooth, and Degree of Freedom (DoF) are 106,500 km, 1.083, and 6.539, respectively. This shows that precipitation is relatively smooth at a large range of approximately 106,500 km. It also means that there is a spatial dependency of weather stations up to 106,500 km. In other words, the range value implies that after this distance, there seems to be no correlation or similarities in the magnitude of extremes in the NCRB. This is crucial information for management activities such as flood control.

This study incorporates relevant spatial and temporal covariates to model a heterogeneous-extreme watershed. The derived 25- and 50-year return levels for the daily annual maxima precipitation were modeled through considerations across the weather stations inside the basin. Interestingly, the occurrence of annual precipitation amounts is well captured within the basin. In other words, studies have consistently observed increases in the frequency and intensity of heavy rainfall events at the weather stations in N.D. and M.N. (inside the RRB). Still considering the RRB and from a climatic projects standpoint, it has been estimated that under a higher emissions scenario (A2), for 2071–2099 projections, General Circulation Model (GCMs) have projected average winter and spring precipitation to increase by 10–20% relative to the base years (1971–2000) (Pryor et al. 2014). These base years are within the period considered in this study.

Spatial modeling of daily annual extreme precipitation is challenging in theoretical and computational aspects and is the main reason why it is not as popular as conventional variogram-based geostatistics. Therefore, there is still room for more research from both a theory and applications standpoint. This research shows a temporal trend in the extreme daily precipitation. The lack of statistically significant inclusion of time in the final model is interesting. There is an apparent reduction and increase in extreme precipitation in the various sections of the basin. Future studies should focus on these two different sections (regions) individually. More on future research, it would be interesting to perform this analysis at different extreme precipitation accumulation rates (1-, 6- and 12-h). This study used 24-h accumulated extreme precipitation. Furthermore, it would be interesting to compare the results from this study with those based on the Bayesian hierarchical modeling technique, where spatial information can be pooled from the latent Gaussian processes specified on the distributional parameters. Another key issue this study has revealed is the lack of meteorological stations (or substantial stations with significant missing values) in many of the sub-basins in the NCRB. Therefore, there should be caution in interpreting the return levels in the sub-basins with relatively low weather stations. This calls for more collaborative efforts to improve data collection and establish more weather stations in northern Canada. This effort is required from both federal, provincial, and territorial jurisdictions.

The steps and modeling processes used in this study can be reproduced at other river basins in Canada or beyond. However, care should be taken when defining the covariates that are important in determining the trend surface since underlying climatic driving factors may be different. Further research should consider the impacts of other covariates such as the North Atlantic Oscillation, NDVI, Melt Index and Timing Factor (Boluwade & Rasumessen 2015). It would be of high interest in further research to quantify the inclusion and impacts of these covariates in improving the extremal behavior of daily annual extreme precipitation in the NCRB. Finally, under the United Nations Framework Convention on Climate Change Paris Agreement and Net-Zero at 2050, there will be more frequent and intense precipitation leading to severe flooding and associated damages.

This study presents stochastic modeling of the spatial extreme of annual precipitation maxima of the NCRB. Annual maxima data were extracted from the Global Historical Climatology Network Daily and Environment and Climate Change Canada. In this study, the GEV parameters are expressed as linear combinations of geographical coordinates (e.g., longitudes and latitudes) and topography. Additional covariates, such as time, were considered and included in defining the GEV parameters. Takeuchi's Information Criterion (TIC), which is most relevant when the model is misspecified, was used to select the most preferred model among nine others. Whenever extreme events occur, rainfall values from weather stations that are in close proximity geographically will have similar values, meaning dependency structures integration in extreme precipitation predictions is paramount. This study applied the Extremal-t max-stable process (MSP) model in predicting spatial extremes in the NCRB. The results show that this MSP model captured the spatial dependency structure of the extreme precipitation in the NCRB using the Quantile–Quantile (Q–Q) plots of the return levels and other metrics.

The following conclusions are drawn from this study:

• (1)

  The precipitation extremes in the NCRB show both spatial trends and dependency. The southern portion of the basin (i.e., Red River Basin (RRB) and portions of Winnipeg River Basin (WRB)) has shown an increase in extreme trends, which is consistent with other findings about the increase in extreme events in this region.

• (2)

  The model GEV parameters were expressed as simple linear combinations of geographical coordinates (i.e., longitude and latitude) and topography. Including topography shows that the Rocky Mountains have some influence on the extreme precipitation in the NCRB. Further, the temporal covariate (i.e., time) was included in the MSP model which captured well the trend and dependency of spatial extremes in the NCRB.

• (3)

  The Q–Q diagnostic plots of return level for some randomly selected weather stations show that the predictions are within the 95% confidence bound, meaning that the fitted Extremal-t MSP adequately captured the spatial dependency structures of the annual daily extreme precipitation in the NCRB.

• (4)

  Considering the eight sub-basins inside the NCRB, 25- and 50-year return-level scenarios for the daily annual extreme precipitation show spatial variability in all the sub-basins, except the RRB and WRB, where the return level is relatively uniform and maximum across the entire weather stations. This is consistent with other studies of increased extreme precipitation in these two sub-basins.

Overall, this study is critical in understanding the environmental statistics of a large watershed in a cold climate. It also extends the research on understanding the impacts of spatial–temporal covariates in capturing and reproducing the daily maximum extreme precipitation in a large basin. Furthermore, for an optimal environmental risk assessment, it is crucial to understand both local and regional spatial extreme precipitation dependency and trends from the standpoint of several weather stations; therefore, the results from this study will be useful in the study of extreme areal rainfall for flood protection, climate change adaptation, and resilience in North America.

A.B. did all the work that produced this paper and approved the manuscript.

Datasets are publicly available at: https://www.ncei.noaa.gov/products/land-based-station/global-historical-climatology-network-daily and https://climate.weather.gc.ca/historical_data/search_historic_data_e.html repositories.

Datasets are publicly available at: https://www.ncei.noaa.gov/products/land-based-station/global-historical-climatology-network- daily and https://climate.weather.gc.ca/historical_data/search_historic_data_e.html repositories.

The authors declare there is no conflict.