Nordic contributions to stochastic methods in hydrology

In 1971, 50 years ago, stochastic approaches to hydrological problems were relatively sparse in the Nordic region. By inspiration from the international literature and increasingly more frequent international cooperation combined with actual needs, a fruitful period was initiated, resulting in many important contributions to stochastic methods in hydrology from the region.

In this review, the focus is on the scientific contributions from the Nordic and Baltic countries. We have therefore concentrated the efforts on studies, where at least one of the authors has a Nordic or Baltic affiliation. The literature was screened by searching in abstracts databases (e.g. Scopus) combined with key references from the author team's own knowledge. It should be noted that stochastic approaches are used to solve practical needs for society. Therefore, a large amount of ‘grey’ literature exists in form of technical reports, which only exceptionally have been considered this review.

The main purpose/aim of this paper is to provide a broad review of the development of stochastic methods by focusing on stochastic modelling of floods, droughts, rainfall extremes, water resources management (WRM) and climate change impacts. The first publications emerged in the 1970s and the development during time is indicated in Figure 1 by the number of referenced publications distributed on decades. In the first two decades, the activity was modest followed by a factor 4 in the next two. After 2010, the activity is again significantly multiplied revealing a strongly growing interest in the subject and the need for further development. In the conclusions, we summarize emerging trends in approaches and methods applied by the scientific community as a response to societal needs.

Estimation of flood magnitudes of low annual exceedance probability or large return periods is needed for design purposes to reduce the potential risk from flood hazards. Design flood estimates are used for dam safety, area planning, bridges and levees, etc. As flood frequency analysis is a crucial part of risk management, approaches and methods for design flood estimation are published in national guidelines and reports, i.e. in the ‘grey literature’. Several of these approaches are summarized in Castellarin et al. (2012). In this review, the focus will be on the scientific publications.

Two basic approaches for statistical analysis of floods can be identified: (1) statistical flood frequency analysis based on flood data and (2) precipitation–runoff approaches where statistical analysis of extreme precipitation and possibly snowmelt and catchment moisture content are used to derive the statistical distribution of floods. Examples from both approaches will be reviewed. Since design flood estimates are used to plan land use and design infrastructure to sustain floods for the coming 50–100 years, climate change impact on floods is in recent years extensively addressed in the literature by analysing past trends and linking floods to large-scale weather patterns and flood-generating processes, as well as assessing future projections of floods. The latter is addressed in a separate section on climate change and adaptation efforts.

Most studies aim to estimate the design flood specified as a T-year flood (Q_T). Alternatively, as suggested in Rasmussen & Rosbjerg (1989), the design flood can be determined by specifying the probability that a flood is exceeding a given value in a period of L years, R_L, where L is the lifetime of a structure. In Rasmussen & Rosbjerg (1991a) it is further demonstrated that using the distribution of R_L might be more relevant than using the distribution of the T-year flood for design purposes, since this probability can directly be used in risk assessments. Jiang et al. (2019) and Haghighatafshar et al. (2020) argue that the probability over lifetime is more relevant to use in a non-stationary context, where the annual exceedance probability of floods varies with time.

To characterise flooding, most studies use flood discharge (m³/s). However, other processes than high discharges might cause large inundations. In Nordic climate, ice-jams might cause severe flooding. Pagneux et al. (2010) therefore propose to use inundation extent as a key variable to analyse magnitude and frequency of flooding events. Using flood discharge will underestimate the flood risk in areas where ice-jams create large flooding.

At locations with streamflow measurements, statistical flood frequency analysis involves fitting a distribution function to either an annual maxima series (AMS) or a series of peaks-over-threshold (POT). The latter is often denoted partial duration series (PDS) in hydrology. The POT approach is based on extracting flood peaks above a predefined threshold. To estimate a T-year flood, both the distribution of the flood magnitudes and the frequency of exceedances need to be assessed.

The first challenge when using the POT approach is to select the threshold. In Rosbjerg & Madsen (1992), a threshold level based on a frequency factor is recommended, whereas Madsen et al. (1997a) used the 2% percentile. Alternatively, a threshold resulting in a predefined average number of events per year can be used. Madsen et al. (1997a) recommend a minimum of 1.64 floods per year to get more precise design flood estimates from POT than from AMS.

The second challenge is to account for the dependence between flood peaks. In Rosbjerg (1977b, 1985, 1987), the dependency is accounted for when estimating the annual exceedance probability. In Rosbjerg & Madsen (1992), another approach is used by selecting only independent flood peaks. The independence is achieved by following the recommendation by USWRC (1976) that subsequent flood peaks should be separated by at least 5+In(A) days, where A is the catchment area in square miles, and in addition the inter-event discharge should decrease by at least 75% of the lowest of two succeeding flood peaks. A third challenge when using the POT approach is to account for seasonal variations in floods. However, in Rasmussen & Rosbjerg (1991b) it is shown that for most data sets it is not necessary to account for these variations to estimate a T-year flood.

The annual maximum approach is based on extracting the largest flood each year and fit a suitable distribution. AMS is less systematically evaluated in the Nordic–Baltic scientific literature, nevertheless it is used in several operational guidelines (Castellarin et al. 2012). Madsen et al. (1997a) compared the AMS and POT approaches and showed that for most applications in hydrology, the POT approach provides the smallest estimation uncertainty for design floods.

In both POT and AMS approaches, it is required to choose a suitable distribution for the flood magnitude and an estimation method for the distribution parameters. Usually, ordinary moments, L-moments or probability weighted moment are used to estimate the distribution parameters (e.g. Rosbjerg et al. 1992); however, in some studies, the maximum likelihood approach is also applied (e.g. Engeland et al. 2004). In more recent years, Bayesian approaches using the Markov Chain Monte Carlo (MCMC) method for estimation is increasingly being applied (e.g. Steinbakk et al. 2016; Kobierska et al. 2018). When a reasonable prior information is available, the Bayesian method is a good choice (Kobierska et al. 2018). The Bayesian approach allows to easily include non-stationarities, calculate uncertainties for design flood estimates and account for uncertainty sources in these estimates.

The generalized extreme value (GEV) distribution or the Gumbel distribution is most commonly used for the annual maximum approach. The GEV distribution has three parameters, whereas the Gumbel distribution, a special case of the GEV distribution, has two parameters. In Kobierska et al. (2018), several distributions are compared, and for Norwegian flood data, it is found that the GEV distribution is the best choice. The exponential distribution or the generalized Pareto (GP) distribution are the corresponding choices for the POT approach. In Rosbjerg et al. (1992), these two distributions are compared using a simulation-based approach. The results show that the exponential distribution generates the smallest estimation uncertainty for data sampled from the GP distribution if the shape parameter is slightly negative.

Since the design flood x_T is an estimate and thus a random variable, the associated uncertainty can be estimated as well. The classical approach is based on assuming that x_T is normally distributed. In Rasmussen & Rosbjerg (1989), Taylor series approximations are used to analytically estimate the uncertainty of the T-year design flood. The approach is further developed in Rasmussen & Rosbjerg (1991a), where a Bayesian formulation of the estimation problem is suggested. However, since x_T might have a skewed distribution, other approaches that account for the skewness have come into use. In Engeland et al. (2004), bootstrapping is applied, whereas Steinbakk et al. (2016) and Kobierska et al. (2018) use a Bayesian approach.

The uncertainty in the estimates of x_T originates from both the observation uncertainty and the sample uncertainty (Petersen-Øverleir & Reitan 2009; Steinbakk et al. 2016). Rosbjerg & Madsen (1998) analysed a wide range of methods to account for sampling uncertainty and concluded that, in general, a reasonable safety margin can be obtained by adding one standard deviation to the estimated design value. For gauging stations where the flood magnitudes mainly are based on a rating curve, high flood observation may introduce estimation biases for return levels, since the estimation is based on combining two highly skewed distributions. The sampling uncertainty, however, is dominating (Steinbakk et al. 2016).

Regional flood frequency analysis is used to estimate design floods at locations with no streamflow measurements. The most common approach is based on an index-flood procedure (Madsen & Rosbjerg 1997a, 1997b; Madsen et al. 1997b; Kjeldsen & Rosbjerg 2002; Kjeldsen et al. 2002; Yang et al. 2010; Hailegeorgis & Alfredsen 2017). The first step of this approach is to define homogeneous regions. For each region a regression analysis between an index flood (i.e. the mean or median flood) is performed, and later all flood data within a region are applied to establish one growth curve. In Gottschalk & Krasovskaia (2002), the regional curve is established by calculating regional L-moments as weighted averages of the L-moments estimated at each streamflow station for annual maxima (AM) data. The weights are proportional to the record length of each time series. In Madsen & Rosbjerg (1997a), POT data are used to establish an index-flood model, where the shape parameter of the GP distribution is constant within a region, and the location parameter and flood rate are site specific. In Madsen & Rosbjerg (1997b), the model is further developed by introducing generalized least-squares (GLS) regression to estimate the parameters of the regional model combined with a Bayesian approach to estimate a posterior model that combine local data with the regional model at each location with observations. In Thorarinsdottir et al. (2018), a fully Bayesian approach to regional flood frequency analysis of AM data is developed, where all parameters in the GEV distribution depend on catchment characteristics. No best model is selected. Instead, a sample of models is established based on the posterior model probability, and the posterior distribution of design floods is achieved by sampling from both models and the posterior parameter distribution. The benefit of this approach is that it is not necessary to define homogeneous regions prior to model estimation. The uncertainty by applying the following different regional estimation methods is analysed by Rosbjerg & Madsen (1995): (i) A standard index-flood approach; (ii) a station-year approach; (iii and iv) two quantile regression approaches and (v) the US Water Resources Council (USWRC 1981) method based on a regional average of the logarithmic skewness. The index-flood and station-year methods were shown to provide the smallest uncertainty. Gottschalk (1989) estimates the empirical exceedance probabilities of regional extremes by combining order statistics with the equivalent number of independent stations based on the spatial correlation of floods.

To analyse trends in flooding, both historical data (Kjeldsen et al. 2014; Retsö 2015; Engeland et al. 2018; Xiong et al. 2020) and paleo-data in the form of lake sediments (Bøe et al. 2006; Støren et al. 2010, 2012; Rimbu et al. 2016; Johansson et al. 2020) are useful, since they increase the length of flood records with 100–10,000 years. In Kjeldsen et al. (2014) and Engeland et al. (2018, 2020) the historical and paleo-data are used to improve the design flood estimates by combining the flood information from both streamflow gauges and historical/paleo-information. The Bayesian inference provides a theoretical framework for combining the different flood data, when estimating distribution parameters and the magnitudes of design floods. In Bøe et al. (2006), Johansson et al. (2020), Støren et al. (2010), Støren et al. (2012), Rimbu et al. (2016) and Engeland et al. (2020) the paleo-data are used to investigate the connection between flood occurrences and large-scale climate variations through the Holocene period (approximately the last 10,000 years). The studies show that both precipitation and temperature explain variations in flood frequency when colder periods imply larger floods in catchments, where snowmelt is an important flood generating process.

In many statistical approaches, several aspects of floods are modelled. One direction is to model the relationship between peak flood and flood volume as a bivariate distribution (Brunner et al. 2016, 2018a; Filipova et al. 2018). In Jiang et al. (2019), this approach is extended to a four-dimensional model by modelling floods of duration 1, 3, 7 and 15 days. These approaches use copulas to model the dependency. In Filipova et al. (2018), catchment properties and flood-generating processes (rain and snow melt) are used to explain the copula parameters. In Jiang et al. (2019), the dependency structure and the marginal distribution are shown to depend on urbanization and reservoir capacity, and the estimated design floods are found to be more sensitive to non-stationarity in the marginal distributions than in the dependency structure.

Instead of modelling floods of different duration as a multivariate extreme, a parallel to intensity–duration–frequency (IDF) models for precipitation has been developed for floods. They are known as runoff–duration–frequency (qdf) models. In Breinl et al. (2021), a parametric qdf-model is fitted to all durations simultaneously without considering the possible dependency between floods of different durations. Then, the duration dependence is explained by catchment properties.

Flood hydrographs might also be needed. In Brunner et al. (2018a), an approach for modelling synthetic flood hydrographs is introduced, and in Brunner et al. (2018b), a regionalization approach is established, where it is shown that the flood peak is relatively easy to regionalize, whereas regionalizing the shape of the flood hydrograph is more challenging, since it depends on flood type as well as catchment properties.

Flooding can also be caused by large waves and/or storm surges. In Harstveit & Jenssen (2003), the relationship between extreme floods and winds is explored.

Precipitation–runoff approaches for design flood estimation has mainly been developed and applied in Norway, Sweden and Finland. Traditionally, an event-based approach based on a simplified Hydrologiska Byråns Vattenbalansavdelning (HBV) model known as PQRUT (Filipova et al. 2016) is applied in Norway. The design flood is calculated using a design rainfall event and saturated initial conditions. The design rainfall is symmetric and should represent the T-year precipitation for all durations in the compound design event. In Filipova et al. (2016), the parameters of the PQRUT model is regionalized based on catchment characteristics. Filipova et al. (2019) further develop the traditional PQRUT approach by introducing a distribution of initial conditions that is combined with the design storm. In Lawrence et al. (2014), a stochastic semi-continuous approach known as SHADEX (Paquet et al. 2013) is applied and evaluated for Norwegian catchments. This approach uses a seasonal dependent mixture distribution of extreme precipitation based on a seasonal threshold and exponential distributions (Blanchet et al. 2015) as well as the probabilities for the precipitation amounts before and after the peak precipitation event. A hydrological model using standard precipitation (and temperature) data is used to establish a continuous simulation. In a first step, flood peaks and the corresponding rainfall events are identified. Subsequently, new precipitation values are systematically sampled for each precipitation event, and new flood peaks are calculated using the hydrological model with initial conditions from the continuous simulation. Thousands of simulations are carried out for each flood event. The seasonal dependent probability of each sampled precipitation event is used to establish an empirical distribution for extreme floods. In the end, SHADEX offers the possibility to analyse the underlying factors driving large floods, i.e. the importance of snow melt and soil moisture, as well as the seasonality of extreme floods (Lawrence et al. 2014).

Slightly different approaches are used for design flood estimation in Sweden (Bergström et al. 1992; Harlin & Kung 1992; Lindström & Harlin 1992; Harlin et al. 1993) and in Finland (Veijalainen & Vehviläinen 2008). In Sweden, the approach uses a 14-day design precipitation event with the largest precipitation amounts on day 9 and constructed so that they roughly contain the largest amounts experienced over 1- to 14-day periods of the 100 years of observations available in Sweden. Such sequences were defined for five regions in Sweden. This design precipitation sequence is systematically inserted into a 10-year climatological record, where a statistical 30-year snowpack has replaced the initial snowpack. In Finland, some modifications to the Swedish approach are introduced by using a 40-year climatological record, where the design precipitation represents a 1,000–10,000 precipitation event and depends on season (Veijalainen & Vehviläinen 2008). A significant advantage of this approach over event-based methods is that the catchment response to extreme precipitation sequences over a range of catchment conditions corresponding to differing saturation states and contributions from snowmelt is sampled during the simulation process. This approach, however, does not provide a probability of the estimated flood magnitudes.

In Gottschalk & Weingartner (1998), an analytic expression for the flood frequency distribution is derived by combining frequency analysis of rainfall volumes with a beta distribution for the runoff coefficient. As pointed out in Filipova et al. (2019), both the SHADEX approach and the stochastic version of PQRUT described above can also be considered as derived distribution approaches, since they aim at sampling from the flood distribution.

A key topic in research has been to identify trends in observed floods and linking trends to changes in river regulations, land use and/or climatic drivers.

Most approaches for detecting trends in floods are based on analysing either POT or AM flood series. Both POT and AM series are used to analyse trends in flood sizes, whereas POT series are also used to analyse trends in flood occurrence (e.g. Mangini et al. 2018). Many studies find trends that are more significant in flood occurrence than in magnitudes (e.g. Wilson et al. 2014; Mangini et al. 2018). Changes in flood statistics can be detected by comparing the empirical and theoretical distributions for different sub-periods (e.g. Sarauskiene & Kriaučiūnienė 2011). To detect the most likely time for a change point in flood statistics, a Cumulative Sum of Departures of Modulus Coefficient (CSDMC) (Chen et al. 2013) or the Pettitt-test (Zhang et al. 2014) can be applied. For analysing more gradual changes, the non-parametric Mann–Kendall test is used to identify significance of trends (e.g. Meilutytė-Barauskienė & Kovalenkovienė 2007), whereas the Sen's slope method is used to quantify the trends in flood data (e.g. Bian et al. 2020). Trends can be also identified in the moments of the flood data (Chen et al. 2019), in parameters of the statistical distribution (e.g. Kallache et al. 2011; Zhang et al. 2014) or in the flood quantiles (e.g. Dahlke et al. 2012). In Kallache et al. (2011), the trends are modelled as linear functions, whereas in Zhang et al. (2014), generalized additive models (GAMs) are used to allow for more flexible and non-linear trends. In Yan et al. (2017), trends are detected by using a two-component mixture distribution, where both the mixture coefficient and the parameters depend on time. Hisdal et al. (2006) point out that due to the large natural variability in floods, detected trends are sensitive to the period that is analysed. Thus, too short flood time series might lead to wrong conclusions.

Variations and trends in flood statistics can be linked to changes in climate and weather like large-scale weather patterns, precipitation, temperature or the underlying flood generating processes (i.e. snowmelt and heavy rain). Several studies link frequency of floods to large-scale circulation patterns, using data outside the Nordic and Baltic regions. Krasovskaia et al. (1999) show that the frequency of floods in Costa Rica depends on the El Niño/Southern Oscillation (ENSO) index, whereas the magnitude does not. Ward et al. (2014) perform a global study showing that ENSO influences the annual flood in river basins covering one-third of the land surface, and in Ward et al. (2016), it is shown that ENSO influences even more the flood duration.

Several studies link non-stationarity in floods to changes in temperature and precipitation. Due to the importance of snow accumulation and snowmelt in large parts of the Nordic region, the interactions between changes in temperature, precipitation and floods are complex. In a review of European trend studies, Madsen et al. (2014) conclude that even if there is an increase in extreme precipitation, no clear trends in floods are detected. However, smaller and earlier snowmelt floods are observed due to increasing temperatures. In a pan-European study, Blöschl et al. (2019) identified two opposite trends in the Nordic and Baltic regions. In areas where the floods are driven by snowmelt (i.e. northern and eastern areas), floods tend to decrease, whereas in areas where the floods are driven by rain (in the western areas), they tend to increase due to increasing precipitation. Similar results are found in national studies. In the Baltic region, a decrease in the spring flood caused by snowmelt is identified (Sarauskiene & Kriaučiūnienė 2011), whereas no trends are found when analysing annual maximum floods (Meilutytė-Barauskienė & Kovalenkovienė 2007). In Norway, Vormoor et al. (2016) analysed trends in the flood-generating processes and detected a transition where in some regions rainfall floods become more frequent and snowmelt floods is decreasing. Slightly different conclusions are drawn in Krasovskaia & Gottschalk (1993). They find that floods are more frequent and slightly larger during warm and/or wet years. In Dahlke et al. (2012), the importance of catchment properties in terms of glacier-covered areas is demonstrated. Similar changes in climate forcing may cause fundamentally different responses of the hydrological systems and flood magnitudes. The floods increased in the glaciated catchment, whereas it decreased in the snow-free ones.

Trends and variations in floods can be linked to changes in land use, in particular urbanization. Bian et al. (2020) demonstrate that urbanization explains most of the changes in floods for a semi-urban catchment in Qinhuai River Basin, Southeast China.

Several approaches are used to define drought events. The earliest studies use the annual minimum flows over accumulation periods from 1 to 15 days (e.g. Gottschalk & Perzyna 1989). As a direct parallel to POT, pits under threshold (PUT), where the minimum value of a low flow sequence below a selected threshold is used in several studies (e.g. Pacheco et al. 2006). To make drought analysis relevant for WRM, drought events can be defined using a threshold level approach (Zelenhasic & Salvai 1987; Tallaksen et al. 1997). This approach defines drought events as periods when the streamflow is below a predefined threshold. These drought events can be characterized by volume, duration and the minimum flow. Tallaksen et al. (1997) point out that using this approach might result in many minor and mutually dependent droughts, when streamflow is oscillating around the threshold. Applying a moving average filter to the original streamflow observations prior to extracting drought events is a simple and efficient way to eliminate such droughts (Tallaksen et al. 1997; Fleig et al. 2006a). Alternative approaches are pooling dependent droughts based on inter-event times and volume (Madsen & Rosbjerg 1995) and the sequent peak algorithm (Tallaksen et al. 1997). The threshold level approach can be extended to specifying seasonal thresholds to account for seasonal variations in climate and runoff generating processes (Fleig et al. 2006b). This can be particularly useful in a Nordic climate, where summer and winter droughts are caused by different processes (Fleig et al. 2006b). Winter droughts are caused by precipitation stored as snow, whereas precipitation deficits and evapotranspiration losses cause summer droughts. Quesada-Montano et al. (2018) propose to use consistent definitions of flood and drought events based on the threshold level method. Then, a joint analysis of drought/flood durations, volumes and trends can be carried out.

Another set of drought indices is based on defining droughts as anomalies from the climatology, e.g. at a monthly time scale. One approach to analyse anomalies is to extend the threshold level approach explained above to a continuous seasonal variation in thresholds (Fleig et al. 2006b) Another approach is to use the standardized streamflow index (SSI), where the monthly streamflows are translated into a probability either by using empirical or parametric distributions (Tijdeman et al. 2020).

The difference between the traditional drought indices and the anomaly-based indices is discussed in Stahl et al. (2020). They state that the threshold-based approach more easily can be linked to low flow events, impacts on ecosystems and triggering of management actions, whereas the anomaly-based approach more frequently is used for monitoring. They demonstrate that the severity of specific drought events depends on the drought definition and points out that the anomaly-based indices are increasingly used and that a more systematic comparison is needed.

Statistical modelling of extreme droughts follows similar procedures as extreme value modelling of floods. Both annual maximum droughts and POT (e.g. Engeland et al. 2004) are used. Engeland et al. (2004) show that using the maxima from block sizes of 2 years might be useful, if there are many years with minor droughts. Also, for the POT method, minor droughts and a few major droughts might cause problems when fitting a distribution. Kjeldsen et al. (2000) and Hisdal et al. (2002) show that applying a threshold to remove minor droughts from the POT dataset is useful. In Clausen & Pearson (1997), annual maximum deficit volumes were used to assess the exceedance probability of a drought event in New Zealand in 1993/4. Gottschalk & Perzyna (1989, 1993) and Pacheco et al. (2006) incorporate streamflow recession characteristics to improve local fits of probability functions to low flows and show that in particular the length of the dry spells is a key factor to define homogeneous regions.

For predicting droughts of long return periods, the standard extreme value distributions like the GEV, Gumbel, generalized logistic (GL), GP, exponential, log Pearson type II and gamma distributions are used. The Tweedie distribution is also suggested as a distribution to use for drought indices (Stagge et al. 2015).

Both Clausen & Pearson (1997) and Kjeldsen et al. (1999) apply the index-flood approach for regional modelling of droughts. Clausen & Pearson (1997) analyse annual maximum droughts for New Zealand, and they show that a baseflow index is a key covariate for predicting the mean annual drought deficit volume in ungauged catchments. Kjeldsen et al. (1999) apply regional peak over threshold approach with a two-component exponential distribution in Zimbabwe, where L-moment diagrams are used to identify homogeneous regions. In Tallaksen & Hisdal (1997), the catchment-specific weight coefficients of empirical orthogonal functions (EOFs) are used to define homogeneous regions.

Gottschalk & Perzyna (1993) interpolate low flow distribution along a river by interpolating the recession characteristics, and in Pacheco et al. (2006) these interpolated distributions are used to regionalize low flow statistics in Costa Rica. Yu et al. (2014) further develop this approach by accounting for dependencies between the recession parameters and dry spell durations in the derived distributions.

Drought modelling by severity–area–frequency curves is developed in Hisdal & Tallaksen (2003). The model establishes the probabilities of an area to experience a severe drought. Due to limited amount of data, a simulation-based approach is used to establish these curves. Such approach is useful for water resource management, where droughts are handled at a regional level.

Several research papers are motivated by detecting trends in droughts. Hisdal et al. (2001) study trends in droughts all over Europe using both annual maximum and peak over threshold data. The Mann–Kendall test is used to assess the significance. They find some trends within the period 1962–1990. Increasing drought-deficit volumes were found in Spain, the eastern part of Eastern Europe and large parts of the UK, whereas decreasing drought deficit volumes occurred in large parts of Central Europe and the western part of Eastern Europe. Changes in precipitation or artificial influences might explain these trends. Wilson et al. (2010) find a tendency towards large water deficit in south-eastern Norway. Yu et al. (2018) used PUT data to analyse non-stationarities in the Loess Plateau in China and found that non-stationarities are caused by changes in recession duration and can be linked to changes in land use, particularly check dams, forest cover and grass land cover.

Extreme rainfall can be viewed in terms of virtually all scales in time and space. We can talk about an extremely high 1-min intensity as well as an extremely wet year; we can talk about an extreme rainfall observed in a specific point as well as over a 50,000 km² area. Generally, when considering extreme rainfall in a societal context, e.g. for infrastructural design, the focus is on short-duration extremes, at sub-daily time scales, within a small, local area, where a single rainfall gauge adequately represents that. This is the focus of the following review, although daily and areal extremes are also covered to some extent.

Local, sub-daily rainfall extremes are best observed using rainfall gauges with a short time step, generally a tipping-bucket or weighing type gauge. Examples of single gauges that have been operational for many decades, up to a century, exist in all Nordic countries. In terms of national sub-daily observations, the main Nordic networks are described in Table 1.

An additional key source of high-resolution rainfall is weather radar, which has a high resolution in time and space. Through the BALTRAD initiative, composites with near-full coverage of the study domain are generated at resolutions 15 min and 2 km × 2 km from a network of C-band radars since the early 2000 s (Michelson et al. 2000). However, radar-observed rainfall is affected by a range of uncertainty sources limiting the space–time accuracy, especially of high intensities (e.g. Schleiss et al. 2020). Besides C-band radars, Denmark furthermore has a network of X-band radars with even higher-resolution, and single X-band radars also exist in Finland and Sweden (e.g. Thorndahl & Rasmussen 2012). New opportunities for high-resolution rainfall observation include attenuation in commercial microwave link networks (e.g. van de Beek et al. 2020) and private weather stations online (citizen science).

In a hydrological engineering context, the observations are commonly processed to estimate IDF (or depth–duration–frequency; DDF) statistics. IDF statistics provide the rainfall intensity or depth associated with a specific combination of duration (i.e. accumulation period) and frequency (i.e. return period). The estimation typically includes the following steps: (1) extraction of extreme values from a time series, either AM or PDS ( a.k.a. POT); (2) fitting of a theoretical frequency distribution, typically GEV to AM values or GP to PDS values and (3) fit a generalized equation that estimates intensity or depth as a function of duration and return period. Also, high percentiles of empirical frequency distributions from rainfall observations are sometimes used to quantify extremes, e.g. in climate analyses.

Pioneering work in the field of rainfall extremes include Rosbjerg (1977a), who investigated event arrival assuming a Poisson process and established relations between return periods based on AM and PDS, respectively. In Sweden, Dahlström (1979) proposed a national model for sub-daily return levels with dependence on summer mean temperature and total precipitation. In the city of Lund, Sweden, a network of 12 tipping-bucket gauges was established and was operational from 1979 to 1981. The observations were used to characterize small-scale dynamic and aerial rainfall properties (Bengtsson & Niemczynowicz 1986). A key concept in rainfall extremes is Probable Maximum Precipitation (PMP), i.e. the greatest depth of precipitation for a given duration meteorologically possible for a design watershed or a given storm area at a particular location at a particular time of year, with no allowance made for long-term climatic trends (WMO 2009). Førland & Kristoffersen (1989) applied PMP in Norway using both meteorological and statistical methods.

Since the mid-1990s, when ∼15 years of observations from the national gauge network were available, substantial and significant research and development with respect to short-duration extremes have been performed in Denmark. Arnbjerg-Nielsen et al. (1994) derived IDF and DDF statistics, including uncertainty, and investigated differences between the used 56 stations as well. Another statistical analysis of Danish data, including geographical variations, was performed by Harremoës & Mikkelsen (1995) as the first in a three-part paper sequence. In the second part, the uncertainty in IDF statistics was divided into sampling variability and regional variability, where the latter was found to be statistically significant (Mikkelsen et al. 1995). In the final third part of the paper, a method for estimating the spatial inter-site correlation structure was developed and applied (Mikkelsen et al. 1996). The evidence of spatial variability found in these studies spawned several further investigations focusing on regionalization, see the section Regionalization for key papers in this effort.

A fair number of studies have focused on properties of the methods used in IDF analysis, including parameter values and estimation methods. (Madsen et al. 1997a) performed a theoretical study to identify the optimal combination of IDF/DDF-approach (AM/GEV and PDS/GP) and parameter estimation method (maximum likelihood, moments and probability weighted moments). The optimal combination was found to depend on the shape parameter. Ragulina & Reitan (2017) compared GEV shape parameters for daily time series in Norway with estimates for other international regions. In another international (global) study of daily observations, PCA was used to develop and evaluate a relationship between wet-day mean precipitation and high percentiles, e.g. 95% (Benestad et al. 2012).

As demonstrated above, the Fennoscandian countries apply different methods to establish national IDF statistics. The latest available estimates for present climate are seen in Table 2.

IFD curves are typically used in design of urban drainage systems. Traditionally, this has been done on a regulatory basis by requesting certain T-year design values depending on construction type. Alternative design methods using economic optimality have recently attracted more attention (e.g. Rosbjerg 2017).

Locally estimated IDF statistics based on observations from one single gauge will inevitably become highly uncertain. More data, and thus less uncertainty, can be included by pooling observations from nearby gauges according to the station-year method (e.g. Buishand 1991). However, the extremes will be statistically different at some distance depending on differences in the rainfall regime caused by, e.g. temperature and topography. Therefore, many investigations have focused on developing methods for estimating IDF statistics at a given location, even if ungauged.

As previously mentioned, a lot of work on regionalization has been performed in Denmark. Madsen et al. (1994) concluded that spatial homogeneity could not be justified and used the Bayesian theory to perform regional PDS modelling with both exponential and generalized Pareto distributions. The approach was shown to significantly reduce uncertainties as compared with using at-site data (Madsen et al. 1995). In Madsen et al. (1998), least-squares regression related PDS model parameters to climatic and physiographic characteristics. A clear dependence on mean annual rainfall was found as well as a metropolitan effect in the Copenhagen area. The approach was further developed and evaluated in Madsen et al. (2002). Gregersen et al. (2013a, 2017) used generalized linear modelling to assess the value of spatial and temporal explanatory variables for the IDF estimation and found mean summer temperature and precipitation particularly useful. Recently, non-stationarity has been introduced in the regionalization for Denmark due to indications of increasing extremes (Arnbjerg-Nielsen 2006). Gregersen et al. (2017) extended the method of Gregersen et al. (2013a) to allow for non-stationarity also considering more parameters. Madsen et al. (2017) related PDS parameters to daily gridded rainfall statistics, including non-stationarity as an additional source of uncertainty.

In Norway, Dyrrdal et al. (2015) used a Bayesian hierarchical model to characterize the spatial variation of hourly precipitation extremes including uncertainty. A range of physiographical and climatic covariates was used to estimate GEV parameters based on generalized linear modelling. For estimation of sub-daily rainfall, Førland & Dyrrdal (2018a) divided Norway into seven regions. Olsson et al. (2019) performed a clustering of locally estimated GEV parameters to divide Sweden into four regions. For each region, DDF statistics including uncertainty were estimated by the station-year method. The regionalization of Swedish DDF statistics by Dahlström (1979), based on summer temperature and precipitation as covariates, has been further developed into an updated, generalized model and applied to other European countries as well (Svensson et al. 2020; Dahlström 2021). Recently, a consistent GEV analysis of sub-daily rainfall extremes in the Nordic–Baltic region was made by Olsson et al. (2022). Interpolated 1° × 1° fields with estimated depths for selected return periods and durations were produced and made public in open access.

The fractal theory contains a concept of self-similarity or scale invariance that proved useful in, for example, analysis of rainfall extremes. One key manifestation of this concept is power laws, i.e. statistical log-log linearity. In this context, a well-known power law is Jennings' scaling law of global rainfall maximum depths vs. duration (Jennings 1950; WMO 1994).

Olsson et al. (1993) investigated 1-min observations from Lund, Sweden, and found evidence of scale invariance (or scaling). The scaling properties were dependent on intensity level (i.e. multiscaling or multifractal) as well as aggregation time scale interval. Daily observations in two contrasting climates, monsoon (China) and temperate (Sweden) were analysed by Svensson et al. (1996). A multifractal behaviour was found with properties depending on both time scale, climate and rainfall-generating mechanism.

Multifractal statistical properties imply that the process can be described by a random cascade model, which can be viewed as a gradual zoom-in to higher resolutions (in time or space) with associated redistribution of the considered quantity (rainfall depth, in this case). Olsson (1998) developed a random cascade model for temporal rainfall disaggregation, where redistribution parameters are dependent on rainfall volume and position in the rainfall sequence. In Güntner et al. (2001), the model was evaluated in two different climates (UK, Brazil) with particular focus on extremes as well as parameter transferability in time and space. The model has later been further developed and used for various applications related to extreme rainfall, e.g. soil erosion (Jebari et al. 2012).

Despite their obvious importance for society, areal rainfall extremes have so far received rather limited attention in the study domain. This is probably related to limited availability of appropriate empirical data, i.e. observations with a high spatial resolution (e.g. a dense gauge network) and accurate intensity. Improved accuracy in radar-based rainfall products is expected to spawn more investigations and development focusing on areal extremes from now on.

A concept relating point-scale to areal extremes is areal reduction factors (ARFs), where a point-scale intensity is reduced to estimate the associated mean intensity within a surrounding area (e.g. review by Svensson & Jones 2010). Bengtsson & Niemczynowicz (1986) estimated ARFs using 1-min observations from 12 gauges in Lund, Sweden in a small-scale investigation. Recently, Thorndahl et al. (2019) estimated storm-centred ARFs from Denmark's 15 years of radar data. A generic relationship was developed between ARF, duration (1 min–1 day) and area (0.1–100 km²).

In Norway, estimation of larger-scale, daily areal extremes has been the focus. Førland & Kristoffersen (1989) applied the growth-factor method (NERC 1975) to characterize large-scale extreme events in catchments up to 5,250 km². Skaugen et al. (1996) developed a method to simulate areal extremes by means of fractionally divided rainy areas in catchments up to 48,000 km². Dyrrdal et al. (2016) used the GEV distribution to describe areal extremes from a gridded 1 km × 1 km dataset of daily precipitation. The GEV shape parameter was found to depend on the precipitation-generating mechanism.

A category of investigations has focused on identifying and analysing relationships between rainfall extremes and the large-scale atmospheric circulation, often represented by climatic patterns or indices. These efforts aim to understand the physical causes and large-scale spatial dynamics behind the extremes, including analyses of temporal variability and trends, which are relevant for assessing climate variability and change.

Hellström (2005) related extreme and non-extreme daily precipitation in Sweden to weather types as well as atmospheric variables. Extremes were found to be associated with the cyclonic weather type as well as higher vertical velocity and specific humidity. Gustafsson et al. (2010) analysed moisture trajectories associated with summer extremes in Sweden and concluded that Europe and the Baltic act as important source regions. Azad & Sorteberg (2017) investigated moisture fluxes over the North Atlantic Ocean in connection with extreme precipitation events at the Norwegian west coast. Almost all events proved to be generated by atmospheric rivers, i.e. narrow plumes of intense high-level moisture.

In some investigations, climatic teleconnection patterns have been used for establishing links to rainfall extremes. The East Atlantic (EA) pattern was found to be a significant explanatory variable in the model of Gregersen et al. (2015) for short-duration extremes in Denmark. Irannezhad et al. (2017) investigated spatio-temporal changes in a suite of precipitation indices in Finland during 1961–2011 as well as their relation to teleconnection patterns. Significant correlations to several patterns were found, including EA, Scandinavia and Polar patterns. Marshall et al. (2020) studied the spatio-temporal variability of precipitation extremes in Artic Fennoscandia during 1968–2017. In addition, distinct relationships with several teleconnection patterns were evident, with both seasonal and geographical dependencies.

Hydrosystems are typically designed and operated based on the stochastic nature of triggering water cycle components. Therefore, different aspects of WRM, such as water supply and demands, groundwater issues or coping with extreme events, must be treated stochastically for making a sustainable decision (Smith 1973). In addition to the stochastic nature of hydrologic elements, a complex network of demands, supplies and stakeholders, which involves economic and social factors, leads to stochastic approaches capable of considering uncertainties germane to each component (Dralle et al. 2017). In this section, we briefly review some of the studies that have utilized stochastic methods to solve different issues in WRM. The studies are classified into five main WRM categories of (i) water supply, allocation and demand; (ii) droughts and floods; (iii) reservoirs and hydropower plants; (iv) groundwater management and (v) risk assessment and data assimilation (DA).

Water distribution between different beneficiaries is a challenging task in WRM due to the stochastic features of available resources and uncertainties of users' demand (Gaivoronski et al. 2012b). The problem gets more complicated in arid and semi-arid regions suffering from water scarcity. Several stochastic methods have been used to tackle water supply, allocation and demand problems. The associated research relies on both classic and data-driven techniques and covers different issues from both quantitative and qualitative perspectives.

From a quantitative point of view, Knudsen & Rosbjerg (1977) demonstrated that the general dynamic programming could be utilized for optimal scheduling of water supply considering the cost of operation. Selek et al. (2013) applied a consecutive high dimensional stochastic non-linear mixed-integer method to assess water distribution systems under water-demand uncertainties and develop an optimal control policy. Sveinsson (2014) showed that stochastic synthetic time series could be beneficial to estimate an irrigation system's performance under uncertainties of water deliveries. The author also highlighted the potential use of stochastics models for water supply and demand forecasting. Sechi et al. (2019) implemented the stochastic quasi-gradient optimization technique to optimize multi-user and multi-reservoirs in water supply systems and develop a robust pumping decision strategy in water supply systems. Most recently, Golmohamadi & Asadi (2020) suggested a multi-stage stochastic approach to schedule a time-oriented agricultural demand. Hatsey & Birkie (2020) showed that failure in pumping systems could be simulated via stochastic models.

Water quality has always been a targeted concern in WRM. Access to high-quality water is vital for human health and is an essential component of sustainable social and economic development and the protection of the environment. The WRM literature has witnessed several stochastic approaches to solve water quality issues. Some examples include the use of (i) stochastic data-driven models, such as artificial neural network, support vector machine and decision tree-based models to predict the concentration of phosphate (Latif et al. 2021), (ii) stochastic model predictive control to optimize the energy-demanding aeration process for nitrogen removal (Stentoft et al. 2019), (iii) stochastic mobility model to monitor contamination in water distribution networks (Du et al. 2016) and (iv) stochastic failure process to assess the aggregated risk of multiple contaminant hazards in groundwater wells (Enzenhoefer et al. 2015).

Droughts and floods are of paramount importance stochastic natural events that can be considered benchmarks for WRM (Van Loon et al. 2016; Quesada-Montano et al. 2018). These events may significantly affect the quality and quantity of water resources and pose a great challenge for WRM in many regions. The relevant literature shows different approaches implemented to forecast, monitor and alleviate their negative impacts. Regarding groundwater vulnerability, Amundsen & Jensen (2019) considered drought events as a stochastic threat and showed that an increase in the probability of drought might yield counteracting effects. The authors recommended a drought-based policy for optimal groundwater extraction and thus dynamic groundwater management. Brunner & Tallaksen (2019) conducted a ground-breaking study on the probability of multi-year droughts across Europe. The authors demonstrated that the proneness of European catchments to droughts could be realized by applying phase randomization and a flexible four-parameter kappa distribution in the stochastic simulation of historical streamflow records. It is worth mentioning that statistical justification for persistent droughts is often tricky in typical stochastic approaches. This can result in a dramatically underestimated risk of failure when forecasting the reliability of reservoirs. The Hurst phenomenon offers a consistent basis to remedy this. The Hurst coefficient and Monte Carlo simulation can be used to quantify the amount of persistence in a time series for reservoir operation (Cox et al. 2006). Similarly, stochastic regression models and spectral analysis were also suggested to simulate teleconnections between large-scale atmospheric circulation patterns and catchment hydrology to predict extreme weather events (Räsänen & Kummu 2013).

Energy and water resources are strongly tied; energy is necessary to deliver water, and water is required to generate energy (Pereira-Cardenal et al. 2016). Both amounts of water and the associated generated energy in hydropower deal with the stochastic behaviour of key hydrological variables, such as precipitation, streamflow and groundwater (Li et al. 2019). Therefore, standalone stochastic methods and, in some cases, their hybridized version with optimization techniques, such as stochastic dynamic programming (SDP), water value method and stochastic time series modelling, have been frequently applied to optimal design and operation of reservoirs (e.g. Sveinsson 2014; Davidsen et al. 2015a; Pereira-Cardenal et al. 2015). For instance, Pereira-Cardenal et al. (2015) applied the water value method and power market models through SDP to minimize the cost of hydropower and maximize irrigation profit. In another study, Pereira-Cardenal et al. (2016) demonstrated that SDP could also be used to attain optimum water allocation by considering the inflow, reservoir capacity, hydropower generation and irrigation. Bauer-Gottwein et al. (2016) showed that SDP could be applied to solve multi-objective optimization problems in reservoir systems. The SDP combined with the hybrid linear programming optimization technique was applied for cost-optimal water quantity and quality management and water allocation. Reservoir releases, groundwater pumping, wastewater treatments and water curtailments have been considered by Davidsen et al. (2015b). In a more recent study, SDP was utilized to estimate the value of water and, therefore, simulate and forecast hydropower generation (Pérez-Díaz et al. 2020). Hydrological and meteorological forecasting could support reservoir operation for drought management, flood warning and water allocation (Gelati et al. 2011; Acharya et al. 2020). Optimization techniques (e.g. the shuffled complex evolution (SCE) algorithm) can be applied for reservoir operation (Le Ngo et al. 2007, 2008). A Bayesian probabilistic forecasting model was used to assess the uncertainty of the flood limiting water level control for reservoir operation in real-time dynamic (Liu et al. 2015). A compound model has been applied to optimize the reservoir rule curve; for this purpose, a simulation was coupled with a genetic algorithm and inner linear programming (Taghian et al. 2014).

Regarding stochastic time series models, several studies are available in the literature. For example, a Markov-switching model using exogenous inputs and autoregressive models was applied to forecast monthly inflow to perform stochastic optimization of reservoir release (Gelati et al. 2011). Grey-box is another robust stochastic modelling technique used for runoff prediction and quantification of the associated uncertainties in urban drainage systems (Löwe et al. 2013). In case of a shortage of historical measurements, generation of high-resolution stochastic time series of key hydrological variables was recommended by Sveinsson (2014). The author showed that synthetic but high-resolution time series could be used to assess the probability distribution of the variables and predict short-term operation rules and hydropower generation, which are crucial for securing WRM. The climate change impacts on these time series were also considered in recent studies (Pereira-Cardenal et al. 2014; Sørup et al. 2017). Sechi et al. (2019) used the stochastic quasi-gradient method to optimize and simulate a pumping system's often conflicting water deficit and energy consumption.

Groundwater resources constitute nearly 30% of the Earth's fresh water and are critical to supporting aquifers' long-term viability and protecting their nearby surface waters. Likewise, surface and atmospheric water resources, several studies have used stochastic methods to cope with quantitative and qualitative groundwater management issues. Examples include (but are not limited to) the use of (i) a risk-cost minimization model based on groundwater value to develop a decision support system and find alternatives in order to reduce the contamination in drinking water (Rosén et al. 1998), (ii) a discrete-time stochastic model to tackle the problem of managing groundwater under random recharge in a single-cell aquifer (Krishnamurthy 2017), (iii) a meta-analytical method based on a stochastic simulation of water/solute transport to assess the impact of land use change on groundwater quality and (iv) applying multi-objective optimization for well field management and operation (Dorini et al. 2012; Hansen et al. 2012, 2013).

The complicated interaction between different risk-related mechanisms in WRM cannot be easily modelled to characterize the system response under different risk control settings. Flood management is usually addressed from a risk viewpoint to quantify the hazard and depth of inundation and to consider the susceptibility to flooding (Oliver et al. 2019; Andaryani et al. 2021; Darabi et al. 2021). Various stochastics techniques have been used to assess the risk in WRM and related concerns, particularly for extreme events during floods and droughts. For example, Jonsdottir et al. (2005) applied a stochastic simulation model to produce synthetic runoff records with an acceptable length that can be used to estimate occasional droughts and associated risk. Gaivoronski et al. (2012a) showed that coupled cost optimization and risk management models could support the water resources manager to minimize the risk of wrong decisions in WRM. Furthermore, the resilience of the water resources system can be assessed through stochastic programming (Najafi et al. 2020).

DA is a powerful tool for WRM by contributing to real-time hydrological modelling and forecasting (Vrugt et al. 2006; Ridler et al. 2014a; He et al. 2019). The most commonly used approach in DA is the ensemble Kalman filter (EnKF) that has been widely applied for different purposes, such as the interaction of groundwater and surface water (He et al. 2019), integrated hydrological modelling (Rasmussen et al. 2015; Zhang et al. 2015; Ridler et al. 2018), reservoir water level measurement (Pereira-Cardenal et al. 2011), modelling of river hydrodynamics (Schneider et al. 2018), flood forecasting (Butts et al. 2007), groundwater modelling (Drécourt et al. 2006), soil moisture and groundwater interaction (Zhang et al. 2016), risk assessment (Borup et al. 2015) and downscaling satellite soil moisture (Ridler et al. 2014b) and forecasting flows and overflows in urban drainage systems (Lund et al. 2019).

IPCC (2021) states that the frequency and intensity of heavy precipitation events have increased over a majority of land regions with good observational coverage since 1950, and that the increase in frequency and intensity will continue under future global warming.

Evidence that extreme rainfall intensity is increasing at the global scale has strengthened considerably in recent years (Westra et al. 2014). Myhre et al. (2019) concluded that increase in the frequency of heavy rainfall events is the main reason for an increase in total precipitation, and that the increase of intensity is less significant. For Europe and the USA, Benestad et al. (2019), however, found positive trends in rain intensity over the period 1961–2020 at most locations with observation series longer than 50 years.

For Northern and Central Europe, a high number of studies have detected changes in heavy daily rainfall, all adding to the debate on anthropogenic climate change and its potential impact on rainfall extremes (Gregersen et al. 2015). For most parts of Finland, Irannezhad et al. (2017) indicated significant increases in the frequency and intensity of precipitation extremes during 1961–2011. Sorteberg et al. (2018) found mainly positive trends for Norway, when studying changes in the highest measured daily precipitation for the summer season during 1968–2017. Over the Nordic–Baltic region (Finland, Sweden, Norway, Denmark, Estonia, Latvia and Lithuania), Dyrrdal et al. (2021) found that for a majority of stations, there had been a statistically significant increase in the intensity of annual maximum 1-day precipitation during the last 50 years, and also for a majority of the long-term stations (from 1901 to till date). Their results also indicated that the annual maximum 1-day precipitation events now occur somewhat later in the year compared to the beginning of the last century.

In addition to the general observed increase in heavy rainfall, it is also important to understand natural variations imposed on the past (and future) changes in heavy rainfall. In an analysis of return periods, Heino et al. (1999) found that there were high frequencies of ‘extraordinary’ 1-day rainfalls in the Nordic countries in the 1930s and later in the 1980s. Based on smoothed series of daily rainfall extremes from Denmark and southern Sweden, Gregersen et al. (2015) concluded that the frequency of the extreme events shows both a general increase from 1874 to present and an oscillation with a cycle of 25–40 years. The magnitude of the extreme events also oscillates, but with a cycle of 15–30 years and with a smaller amplitude. Similarly, Førland & Dyrrdal (2018a) demonstrated that the positive trend in maximum 1-day rainfall is not monotonic; i.e. decadal variability is superposed on the long-term trend. The long-term centennial variability of heavy 1-day rainfall reveals high values in the 1930s and low values in the 1960 and 1970s. In addition, for the Nordic-Baltic region Dyrrdal et al. (2021) found that smoothed series of annual 1-day rainfall maxima have a cyclic development.

Willems et al. (2012) stated that trend testing has to account for temporal clustering of rainfall extremes. Thus, it is likely that the present high level of heavy rainfall events will be followed by decades with lower values (Førland & Dyrrdal 2018a). Consequently, rainfall design-estimates based on observations from the most recent decade(s) may tend to be biased compared to estimates for 30-year reference periods, e.g. 1971–2000.

Arnbjerg-Nielsen (2006) studied 41 Danish sub-daily rainfall series with an observation period close to 20 years. A statistically significant trend was found for the 10-min maximum intensity towards more extreme and more frequently occurring rainstorms. For the 6-h maximum intensity and total event volume the trends were less pronounced. The findings were confirmed by comparison to physically based climate models and studies based on large regions.

A comparison of the two periods 1979–1997 and 1997–2005 showed a general increase in extreme rainfall characteristics for Denmark (Madsen et al. 2009). For the durations of 30 min to 3 h and return periods in the order of 10 years being typical for most urban drainage designs, the increase in intensity is in the order of 10%. The analysis revealed that the changes are not statistically significant compared with the uncertainties of the regional estimation model, but the increases in design intensities were large and may have significant consequences to the costs of engineering designs.

For Norway, Førland & Dyrrdal (2018b) concluded that during the last 50 years the frequency and intensity of sub-daily rainfall have increased at a majority of the measuring sites. This positive trend also affects rainfall design values using IDF estimates.

Global climate models (GCMs) are the most comprehensive and widely used models for simulating the response of the global climate system to large-scale changes in greenhouse gas emissions. GCMs are not designed to represent local precipitation statistics, but their ability to simulate large-scale features can be utilized to infer local precipitation changes through downscaling (Benestad et al. 2008; Benestad 2010; Olsson et al. 2015; Sunyer et al. 2015a, 2015b). There are two main approaches for downscaling GCM simulation: (1) dynamical downscaling (regional climate models), and (2) empirical-statistical downscaling (ESD).

Dynamical downscaling involves running an area-limited high-resolution regional climate model (RCM) with large-scale variables from a GCM as boundary conditions. The RCMs tend to be expensive to run and may not provide realistic local conditions at small spatial scales (Benestad & Haugen 2007). Several analyses have applied RCMs to study precipitation extremes. Olsson et al. (2012a) outlined a framework for downscaling precipitation from RCM projections to high resolutions in time and space required for urban hydrological climate change impact assessment. Their basic approach was use of the delta change approach, developed for both continuous and event-based applications. Benestad & Haugen (2007) used RCMs to study shifts in the frequency of complex extremes in Norway. Mayer et al. (2018) used the high-resolution EURO-CORDEX ensemble to estimate changes in future sub-daily rainfall in Norway by fitting the GEV distribution to annual maximum precipitation from the simulations. Both stationary and non-stationary methods were applied.

In ESD, a wide range of models and approaches has been used, and this approach may be divided into several subcategories (Benestad et al. 2008): (a) linear models (e.g. regression or canonical correlation analysis), (b) non-linear models or (c) weather generators. The choice of ESD model type should depend on which climatic variable is downscaled, as different variables have different characteristics that make them more or less suitable in a given model. RCMs and ESD complement each other.

Benestad (2007) used ESD to establish (1) a relationship between coefficients in the frequency function of daily rainfall and (2) local projections of mean temperature and precipitation to infer changes in the 95th percentiles of future rainfall for 2050. In a later study, Benestad (2010) explored a new ESD-method for predicting the upper tail of the precipitation distribution. Skaugen et al. (2004) used time-slices of 20 years for 1-day precipitation for present and future climate as training data for a precipitation simulation model. This model was used to generate a time series of 1000 years to assess possible changes in the extreme precipitation regime due to climate change.

Sunyer et al. (2015a) explored three statistical downscaling approaches: A delta change method for extreme events, a weather generator (WG) combined with a disaggregation method and a climate analogue method. All three methods relied on different assumptions and used different outputs from the RCMs. The results of their study highlighted the need to use a range of statistical downscaling methods and RCMs to assess changes in extreme precipitation.

For Denmark, Sørup et al. (2016) modelled precipitation at 1 h temporal resolution on a 2-km grid using a WG. Precipitation time series used as input to the WG were obtained from a network of 60 tipping-bucket rain gauges irregularly placed in a 40 km × 60 km model domain. The WG simulated precipitation time series were comparable to the observation-based extreme precipitation statistics. The climate change signals from six different RCMs were used to perturb a WG. All perturbed WGs resulted in more extreme future precipitation at the sub-daily to multi-daily level, and these extremes exhibit a much more realistic spatial pattern than observed in RCM precipitation output. Overall, the WG produced robust results and is seen as a reliable procedure for downscaling RCM precipitation output for use in urban hydrology.

RCMs and GCMs have shown severe problems with their sub-grid scale parameterizations of convective processes, which affect their ability to reproduce, for example, the diurnal cycle of rainfall intensity, the peak storm intensities and extreme hourly intensities. It is therefore questionable to which extent such RCMs are capable of describing short-duration extremes in the present and future climate (Berg et al. 2019). Westra et al. (2014) stated that present-day GCMs have limited ability to simulate sub-daily precipitation extremes correctly, as they do not explicitly resolve convective processes. This casts strong doubts on future projections of changes in sub-daily precipitation extremes derived from these models. However, the authors stressed that RCMs run at convection-permitting resolutions (CPRCMs) show promising improvements to key attributes of sub-daily rainfall.

Willems et al. (2012) made a critical review of methods for assessing climate change impacts on rainfall in urban areas and concluded that estimation of extreme, local and short-duration rainfall is highly uncertain. Downscaling results from GCMs or RCMs to urban catchment scales is needed because these models cannot accurately describe the rainfall process at suitable high temporal and spatial resolution for urban drainage studies. A review made by Arnbjerg-Nielsen et al. (2013) concluded that many limitations exist in our understanding of how to describe precipitation patterns in a changing climate for design and operation of urban drainage infrastructure

Gregersen et al. (2013b) compared temporal and spatial characteristics from three different high-resolution RCMs to sub-daily rainfall extremes estimated from observed data. All analysed RCM-derived rainfall extremes showed a clear deviation from the observed correlation structure for sub-daily rainfalls, partly because RCM output represents areal rainfall intensities and partly due to well-known inadequacies in the convective parameterization of RCMs. The paper takes the first step towards a methodology by which RCM performance and other downscaling methods can be assessed in relation to the simulation of short-duration rainfall extremes.

For Denmark, Sunyer et al. (2012) analysed five different statistical downscaling methods using results from four RCMs driven by different GCMs. Special focus was given to the changes of extreme events, since downscaling methods mainly differ in the way extreme events are generated. Sunyer et al. (2013) used both gauge data and gridded data to rank different RCMs according to their performance using two different metrics. A set of indices ranging from mean to extreme precipitation properties was calculated for all the data sets. The papers by Sunyer et al. (2012, 2013) highlighted the need to be aware of the properties of observational data chosen in order to avoid overconfident and misleading conclusions with respect to climate model performance, and to acknowledge the limitations and advantages of different statistical downscaling methods as well as the uncertainties in downscaling climate change projections for use in hydrological models.

Berg et al. (2019) evaluated summertime IDF-values of a subset of the EURO-CORDEX 0.11 ensemble against observations for several European countries for durations of 1–12 h. Most of the model simulations strongly underestimated 10-year rainfall depths for durations up to a few hours but performed better at longer durations. Projected changes were assessed by relating relative depth changes to mean temperature changes. A strong relationship with temperature was found across different sub-regions of Europe, emission scenarios and future time periods.

Olsson et al. (2012b) formulated and tested a stochastic model for downscaling short-term rainfall from RCM grid scale to local (i.e. point) scale using data from Stockholm, Sweden. The results suggested that the model might effectively reproduce observed IDF curves. Olsson et al. (2015) explored to what degree sub-hourly IDF statistics in an RCM converge to observed point statistics for Swedish stations, when gradually increasing the resolution from 50 to 6 km. At 50 km, the intensities were underestimated by 50–90%, but at 6 km they were nearly unbiased, when averaged over all locations and durations. In addition, the reproduction of short-term variability and less extreme maxima were overall improved with increasing resolution. At 6-km resolution, a parameterized RCM (RCA3) is in approximate agreement with hourly gauge observations in Sweden (Olsson et al. 2015). Recently, the first historical high-resolution (3 km × 3 km) CPRCM simulations over Fenno-Scandinavia have been made and analysed (Lind et al. 2020; Olsson et al. 2021). The results indicate a substantially improved representation of precipitation, compared with lower-resolution RCM simulation, in terms of both average properties (e.g. diurnal cycle) and extreme properties (e.g. IDF statistics). These results imply that we can anticipate a higher confidence in future precipitation changes estimated by CPRCM projections, which currently are being analysed.

As temperatures increase, the atmosphere can hold more moisture and the potential for more frequent and intense rainfall is present. Most projections indicate that the greatest increases are likely to occur in short-duration storms lasting less than a day, potentially leading to an increase in the magnitude and frequency of flash floods (Westra et al. 2014). Thus, changes in extreme precipitation are expected to be one of the most important impacts of climate change in cities. Based on observations and climate models, Myhre et al. (2019) projected that on a global scale the most intense precipitation events observed today are likely to almost double in occurrence for each degree centigrade of further global warming. Changes to extreme precipitation of this magnitude are dramatically stronger than the more widely communicated changes to global mean precipitation.

For large parts of the Nordic countries, Benestad (2007, 2010) pointed towards an increase in the number of extreme precipitation events, except for the most extreme percentiles for which sampling fluctuations gave rise to high uncertainties. In projections for Denmark, Sunyer et al. (2012) found that three of the four studied RCMs showed an increase in extreme events in the future. The increases in extreme precipitation were higher for higher spatial resolutions and shorter temporal aggregations (Sunyer et al. 2013). In a later study, Sunyer et al. (2015a) assessed the changes and uncertainties in extreme precipitation at hourly scale over Denmark. The results of the three methods applied pointed towards an increase in extreme precipitation, but the magnitude of the change varied depending on the RCM used and the spatial location. In general, a similar mean change was obtained for the three methods.

For Sweden, Olsson et al. (2012b) concluded that the future increase of local-scale short-term rainfall extremes might be higher than the predicted increase of short-term extremes at the RCM grid scale. Olsson & Foster (2014) analysed extreme precipitation for durations between 30 min and 1 day in simulations with the RCA3-RCM for Sweden. The increase in extreme precipitation decreased with increasing duration, and at the daily scale, the percentage values are approximately halved.

For Norway, Hanssen-Bauer et al. (2017) found that events with heavy 1-day rainfall should be more intense and occur more frequently in all regions. By fitting the GEV distribution to annual maximum precipitation from the ensemble of EURO-CORDEX simulations, Mayer et al. (2018) found that also sub-daily extreme precipitation is projected to increase for most areas in Norway. The largest increase was found for higher return periods and shorter precipitation duration.

Cities are becoming increasingly vulnerable to flooding because of rapid urbanization, installation of complex infrastructure, and changes in the precipitation patterns caused by anthropogenic climate change (Willems et al. 2012). Design and optimization of urban drainage infrastructure considering climate change impacts and co-optimizing these with other objectives will become ever more important to keep our cities habitable into the future (Arnbjerg-Nielsen et al. 2013). For urban floods, robust information on changes in extreme precipitation at high-temporal resolution is required to design climate change adaptation measures. Urban hydrological climate change impact assessment requires an assessment of how short-term local precipitation extremes, e.g. expressed as IDF statistics, are expected to change in the future.

Projections of future extreme precipitation are regularly updated. For the Fennoscandian countries the present recommended climate factors (change from present to future climate) are outlined in Table 3.

A last group of studies aims to predict the change in flood sizes in a future climate. A review of both modelling approaches and key findings is provided in Madsen et al. (2014). Most studies use a series of linked models and analyses. The basis is the climate change projections from a Global circulation model (GCM). The outputs from these models need downscaling using either or both RCMs and statistical downscaling to produce precipitation, temperature and other weather variables needed to run a hydrological model. The output from the hydrological models can then by analysed for a reference period and future periods to provide estimates of expected changes. There are many options at each of these steps (several GCMs combined with several RCMs combined with several statistical downscaling methods that again is combined with an ensemble of hydrological model or model parameters), resulting in an ensemble of projections. In Lawrence (2020), the uncertainty is attributed to each of these steps showing that the uncertainty in the statistical flood frequency analysis is of the same magnitude as the uncertainty originating from the different GCM/RCM combinations.

The output from the Nordic and Baltic studies shows that both increase and decrease in flood sizes can be expected. In Vormoor et al. (2014,) trends in seasonality and flood generating processes are analysed, demonstrating that a shift towards more rain floods is expected. In Lawrence (2020), it is shown that for Norway the expected change is linked to the flood generating processes, where a decrease in flood size can be expected in several inland catchments, where snow melt is the main flood generating process, whereas increase in flood sizes is expected in catchments where rain is the main flood generating process. Particularly in small catchments, an increase in flood sizes is expected (Tsegaw et al. 2020). Similar results are shown in Sweden (Bergström et al. 2012), Finland (Veijalainen et al. 2010) and the Baltic countries (Meilutytė-Barauskienė et al. 2010).

When analysing climate change impact on droughts, a modelling chain that includes climate models, dynamical downscaling using a RCM, statistical downscaling and finally a hydrological model is used. In large-scale studies, the downscaling parts might be omitted. Expected changes in drought in Norway are analysed by Hisdal et al. (2006) and Wong et al. (2011) finding substantial increases in hydrological drought duration and drought affected areas in a future climate, especially in the southern and northernmost parts of the country. Reduced summer precipitation is a major factor that affects changes in drought characteristics in the south, while temperature increases play a more dominant role for the rest of the country. Chan et al. (2021) also finds increasing drought severity and frequency in Denmark in a future climate based on emission scenario RCP8.5. Samaniego et al. (2017) compare the uncertainty contribution from GCMs and hydrological models for projection of drought severity for continental catchments in a future climate. They show that the largest uncertainty contribution comes from the GCMs.

In the Nordic countries, the application of stochastic methods in hydrology was initiated in the 1970s and has been gradually growing since them. As seen in this review, Nordic contributions have been manifold, not just fulfilling local needs, but also playing an important role in the international arena. The stochastic methods have been essential for the development of modern hydrology. Today, Nordic researchers are breaking new avenues for development and application of methods to solve problems related to flooding, droughts, heavy precipitation, WRM, and climate change impacts on the water cycle.

We see some emerging trends in both the addressed problems and in the methods that are used. Some key trends we want to highlight are:

• Non-stationarity: At the beginning of the period, stationary approaches were completely dominant. However, gradually it has been recognized that we must also include climate variability and change in the analyses. This sets challenging demands and acts as a vehicle for the future development of stochastic hydrology.

• Design: New design criteria based social cost-benefit analyses and economic optimization should be developed to supplement traditional regulatory-based design.

• Risk analysis: There is an increasing demand from the society to address the risk related to WRM and extremes. New approaches are needed to connect hydrology to impacts and address the probability of adverse impacts.

• Big data/crowd sourcing: Automatic weather stations, radar mapping, smartphones and new measurement technologies, e.g. using satellites or drones, open up for new insights and approaches.

• Machine learning: Due to increasing computer power as well as new data sources, machine learning is becoming increasingly used for hydrological applications and statistical approaches can be seen as special cases of machine learning.

K.E. has acted as the main author of the section Flood frequency and low flow analyses, J.O. has acted as the main author of the section Assessment of rainfall extremes, A.D.M. and A.T.H. have acted as the main authors of the section Stochastic approaches to WRM, E.F. has acted as the main author of the section Approaches to climate change and adaption efforts, while D.R. has undertaken the overall coordination and editing of the paper.

All relevant data are included in the paper or its Supplementary Information.