The probabilistic hydrological MARCS^HYDRO (the MARkov Chain System) model: its structure and core version 0.2

The question of the environmental risks of social and economic infrastructure has recently become apparent due to an increase in the number of extreme weather events. Extreme runoff events include floods and droughts. In water engineering, extreme runoff is described in terms of probability and uses methods of frequency analysis to evaluate an exceedance probability curve (EPC) for runoff. It is assumed that historical observations of runoff are representative of the future; however, trends in the observed time series show doubt in this assumption. The paper describes a probabilistic hydrological MARCS^HYDRO (the MARkov Chain System) model that can be applied to predict future runoff extremes. The MARCS^HYDRO model simulates statistical estimators of multi-year runoff in order to perform future projections in a probabilistic form. Projected statistics of the meteorological variables available in climate scenarios force the model. This study introduces the new model's core version and provides a user guide together with an example of the model set-up in a single case study. In this case study, the model simulates the projected EPCs of annual runoff under three climate scenarios. The scope of applicability and limitations of the model's core version 0.2 are discussed.

1 Introduction

Streamflow runoff serves as a water resource for humans, food production and energy generation, while the risks of water-sensitive economics are usually connected to runoff extremes. In fact, runoff extremes are always connected to human activity since they do not exist in a natural water cycle. Engineering science considers runoff extremes as critical values of runoff that lead to the damage of infrastructure or water shortages, and it introduces the extremes in terms of probability. In particular, in water engineering runoff extremes are evaluated from the tails of exceedance probability curves (EPCs) that are used in risk assessment for water infrastructure and decision-making in cost–loss situations (Mylne, 2002; Murphy, 1977, 1976). The EPC of multi-year runoff allows the estimation of the runoff extremes and supports the designing of building construction, bridges, dams and withdrawal systems, etc.

Modern hydrology uses two approaches to evaluate runoff extremes with their exceedance probability: conceptual modelling (Lamb, 2006) and frequency analysis (Kite, 1977; Benson, 1968; Kritsky and Menkel, 1946). In the conceptual modelling approach, synthetic runoff series are simulated from meteorological series in order to calculate the runoff values of a chosen exceedance probability (Arheimer and Lindström, 2015; Veijalainen et al., 2012; Seibert, 1999). In the frequency-analysis approach, historical yearly time series of runoff are used to evaluate statistical estimators, that is, the mean value, the coefficient of variation (CV) and the coefficient of skewness (CS) (van Gelder, 2006). These estimators are applied to calculate runoff values with their exceedance probability (Guidelines SP 33-101-2003, 2004; Guidelines, 1984; Bulletin 17–B, 1982) needed to support the designing of roads, dams, bridges or water-withdrawal stations. The basic assumption of this approach is that the future risks during an infrastructure's operational period are equal to the risks estimated from the past observations. The runoff extremes are simply extrapolated for the next 20–30 years on the assumption that past observations are representative of the future: the “stationarity” assumption (Madsen et al., 2013).

The number of weather extremes – including hurricanes, wind, rain- and snowstorms, floods and droughts – has increased (Vihma, 2014; Wang and Zhou, 2005; Manton et al., 2001). Historical time series of many climate variables show evident trends which are statistically significant, among which are the series of streamflow runoff (Wagner et al., 2011; Dai et al., 2009; Milly et al., 2005). Rosmann et al. (2016) applied the Mann–Kendall test to analyse a time series of daily, monthly and yearly river discharges for the last four decades. The highest number of trends was detected for the yearly time series of annual runoff. The statistically significant trends are founded on historical time series; thus the water engineers and managers are motivated to revise the basic stationarity assumption that lies behind infrastructures' risk assessment since past observations are not representative of the future (Madsen et al., 2013; Kovalenko, 2009; Milly at al., 2008).

In this paper, we described a method that combines conceptual modelling and frequency analysis in order to estimate runoff extremes in a changing climate. The method adapts the theory of stochastic systems to water engineering practice, and it was further named as advance of frequency analysis (AFA). It was introduced by Kovalenko (1993) and relied on the theory of stochastic systems (Pugachev et al., 1974). The basic idea behind the method is to simulate the statistical estimators of multi-year runoff (annual, minimal and maximal runoff) from the statistical estimators of precipitation and air temperature on a climate scale (Budyko and Izrael, 1991). The simulated statistical estimators of runoff are used to construct EPCs with distributions from the Pearson system (Pearson, 1895). Kovalenko (1993) suggested modelling the EPCs within a Pearson type III distribution based on traditional practice in water engineering (Rogdestvenskiy and Chebotarev, 1974; Matalas and Wallis, 1973; Sokolovskiy, 1968). However, the distribution can also be chosen by fitting (Laio et al., 2009), defined in accordance with local hydrological guidelines (Bulletin 17-B, 1982), or somehow more advanced (Andreev et al., 2005).

A linear “black box” (or a “linear filter model”) with stochastic components is suggested as a catchment-scale hydrological model (Kovalenko, 1993). For this linear model, the theory of stochastic systems provides methods to direct the simulation of probability distributions for a random process (Pugachev et al., 1974). The theory of stochastic systems is applied to analyse and predict runoff extremes on various timescales, ranging from days (Rosmann and Domínguez, 2017) to seasons (Domínguez and Rivera, 2010; Shevnina, 2001), and on various climate scales (Shevnina et al., 2017; Kovalenko, 2014; Viktorova and Gromova, 2008). The AFA approach is a simplification of the theory of stochastic systems on a climate scale. Kovalenko et al. (2010) gave guidelines for water engineers to estimate runoff extremes in a changing climate.

AFA was suggested about 30 years ago; however, a full description of this approach has still not been published in English. Moreover, previous publications in Russian contain many typewriting mistakes in the formulas (Kovalenko, 1993; Kovalenko et al., 2006), and this makes understanding them troublesome, even for native Russians. In this paper, the theory and assumptions of the AFA approach were formulated step by step (see Appendix 1), and the formulas behind the core of the probabilistic hydrological model MARCS^HYDRO (the MARkov Chain System) were accepted for the new version, version 0.2 (see Sect. 1). This model core allows the prediction of a skewness parameter of a Pearson type III distribution. An example of the model set-up, forcing and output for a case study of the Iijoki River is given in Sect. 2. The main features of the model and the limitations of the AFA method are formulated in the Discussions section in order to better place the MARCS^HYDRO model among other hydrological models.

2 Model description

The probabilistic hydrological MARCS^HYDRO model consists of six blocks (Shevnina, 2015). Figure 1 shows the tools for data analysis grouped into blocks: two blocks for the analysis and screening of observed data (DPB and DSB); a block with the model parametrization, cross-validation and hindcasts (“forecasts in the past”; PHP); a block to visualize the model's results (VAB); and a block with socio-economic applications (EAB). Shevnina and Gaidukova (2017) provided details about the algorithms already implemented in each block in the model. In this paper, only version 0.2, for the model's core, is introduced. The formulas behind the model's core version 0.1 are published in the Appendix of the work of Shevnina et al. (2017).

Figure 1 The MARCS^HYDRO model structure and core versions.

The MARCS^HYDRO model simulates three non-central statistical moments of multi-year runoff based on the means of precipitation calculated over a period of 20–30 years. Now, the model's application is only limited by a prediction on the climate scale. The development of a socio-economic infrastructure also needs the climate-scale prediction of river runoff (Milly et al., 2008) because water extremes, such as floods and droughts, lead to economical losses. The AFA approach has found practical applications in building construction (Shevnina et al., 2017; Kovalenko, 2009). The MARCS^HYDRO model allows a quick analysis of the runoff extremes under different climate scenarios. The model needs fewer computational resources because it simulates the parameters of the distribution, while conceptual hydrological models simulate the runoff time series.

The MARCS^HYDRO model parametrization, cross-validation and hindcasts need observations of the river water discharges of a hydrological network for a period in the past (Kovalenko, 1993). For the cross-validation, the yearly time series of river runoff are split into two sub-periods, namely the training period and the control period (Shevnina et al., 2017). The splitting year corresponds to the year when the statistically significant difference in observations within two periods is detected by the Student and Kolmogorov–Smirnov tests (Kovalenko, 1993; Kovalenko et al., 2006). The description of the analysis and screening of the observed river runoff time series, as well as the model cross-validation procedure, fell outside the topics of this paper. We focused on the equations behind the model's core version 0.2 and its limitations.

2.1 Model input

Two blocks of the MARCS^HYDRO model are needed to analyse and screen the observations. The time series of river runoff and precipitation are required for a period as long as possible. However, the length of yearly time series on water discharges does not usually exceed 80–90 years. Hydrological yearbooks or runoff data sets provide observations at sites of national hydrological networks, and the river runoff is expressed as a volumetric flow rate (water discharge, m³ s⁻¹). In the data preparation block of the model, the volumetric flow rate (m³ s⁻¹) is converted to a specific water discharge (ARR, mm yr⁻¹):

$$\mathrm{ARR}=\mathrm{1000}QT/A,$$

where Q is the yearly average water discharge (m³ s⁻¹), T is the number of seconds in a year and A is the catchment area (m²). In the data screening block of the model, the yearly time series of ARR are used in the analysis of homogeneity and trends (Dahmen and Hall, 1990) and to define a period for the model parametrization (called the “reference period” by Shevnina et al., 2017). Then, the reference three non-central moments $m_k(m_k=\mathrm{1}/n\sum _{i=\mathrm{1}}^n{\mathrm{ARR}}_i^k$ for k=1, 2, 3) are estimated from time series of ARR using the method of moments (van Gelder et al., 2006).

The observations on precipitation are collected from meteorological sites, and they may be interpolated into grids in order to better estimate a precipitation rate over a river basin area. In the data preparation block of the model, the mean annual precipitation rate (mm yr⁻¹) is calculated from the observed yearly time series for the reference period. The mean annual precipitation rate for the future period can be calculated from an output of any global/regional climate model or even a set of models. In a study on the catchment scale, the time series of water discharges can be extracted from the Global Runoff Data Centre (GRDC), while the precipitation rate can be estimated from gridded data sets (Willmott and Robeson, 1995). These two data sets were used to perform an example of the model application on the Iijoki River basin.

2.2 Model cross-validation

The MARCS^HYDRO model allows the simulation of the non-central moments of runoff that can be used for the construction of probability distribution (or an EPC); in other words, it provides a probabilistic form of prediction. The end product of the model is the probability density function (PDF) (or the EPC), and there are no simulated time series of runoff to compare with the observations. Kovalenko (1993) suggested comparing the simulated PDF with an empirical PDF by using known statistical tests such as the Kolmogorov–Smirnov test (Smirnov, 1948). In the PHP block of the MARCS^HYDRO model, a specific cross-validation procedure allows conclusions to be drawn about the model's validation and the quality of hindcasts. For the model's cross-validation, the observed time series of river runoff is divided into two sub-periods, namely the training period and the control period. The splitting year corresponds to the year when a statistically significant difference in mean values is estimated over two periods. In this study, we did not pay much attention to the cross-validation procedure since it is described in detail in Shevnina et al. (2019) for the model version 0.2.

2.3 The MARCS^HYDRO model core

In our study, the core version 0.2 of the probabilistic MARCS^HYDRO model was suggested instead of version 0.1 (Shevnina et al., 2017). Version 0.2 allows the evaluation of the skewness parameter of the Pearson type III distribution. In the new core, the non-central statistical moments of the ARR were calculated as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }{m}_{\mathrm{1}}=a-b_{\mathrm{1}},\text{(2)}& {\displaystyle }{m}_{\mathrm{2}}=-b_{\mathrm{0}}-\mathrm{2}{m}_{\mathrm{1}}{b}_{\mathrm{1}}+m_{\mathrm{1}}a,\text{(3)}& {\displaystyle }{m}_{\mathrm{3}}=-\mathrm{2}{m}_{\mathrm{1}}{b}_{\mathrm{0}}-\mathrm{3}{m}_{\mathrm{2}}{b}_{\mathrm{1}}+m_{\mathrm{2}}a,\end{array}$$

where m₁, m₂ and m₃ are the moment estimates of the non-central statistical moments of the ARR; and a, b₀, b₁ and b₂ are the parameters of the distributions of the Pearson equation (Andreev et al., 2005).

To set up the MARCS^HYDRO model, observations of water discharges are needed. For the reference period (denoted by low index r) the moments' estimates for the non-central moments (m_1r, m_2r, m_3r) were first calculated from observed times series of runoff (mm yr⁻¹); then the non-central moments were used to evaluate the parameters of the Pearson equation a, b₀ and b₁:

$$\begin{array}{}\text{(4)}& {\displaystyle }a=\mathrm{0.5}(\mathrm{5}{m}_{\mathrm{1}r}{m}_{\mathrm{2}r}-\mathrm{4}{m}_{\mathrm{1}r}^{\mathrm{3}}-m_{\mathrm{3}r})/(m_{\mathrm{2}r}-m_{\mathrm{1}r}^{\mathrm{2}}),\text{(5)}& {\displaystyle }{b}_{\mathrm{0}}=\mathrm{0.5}(m_{\mathrm{1}r}^{\mathrm{2}}{m}_{\mathrm{2}r}-\mathrm{2}{m}_{\mathrm{2}r}^{\mathrm{2}}+m_{\mathrm{1}r}{m}_{\mathrm{3}r})/(m_{\mathrm{2}r}-m_{\mathrm{1}r}^{\mathrm{2}}),\text{(6)}& {\displaystyle }{b}_{\mathrm{1}}=\mathrm{0.5}(\mathrm{3}{m}_{\mathrm{1}r}{m}_{\mathrm{2}r}-\mathrm{2}{m}_{\mathrm{1}r}^{\mathrm{3}}-m_{\mathrm{3}r})/(m_{\mathrm{2}r}-m_{\mathrm{1}r}^{\mathrm{2}}).\end{array}$$

Then, the parameters of the linear filter model (see Appendix 1 for details), $\stackrel{\mathrm{‾}}{c}$, $G_{\stackrel{\mathrm{̃}}{N}}$ and $G_{\stackrel{\mathrm{̃}}{c}\stackrel{\mathrm{̃}}{N}}$, denoted with a low index r, were calculated:

$$\begin{array}{}\text{(7)}& {\displaystyle }{\stackrel{\mathrm{‾}}{c}}_r={\stackrel{\mathrm{‾}}{N}}_{\mathrm{r}}/\left(a-b_{\mathrm{1}}/\mathrm{2}\right),\text{(8)}& {\displaystyle }{G}_{\stackrel{\mathrm{̃}}{N}r}=-\mathrm{2}{b}_{\mathrm{0}}{\stackrel{\mathrm{‾}}{N}}_{\mathrm{r}}/\left(a-b_{\mathrm{1}}/\mathrm{2}\right),\text{(9)}& {\displaystyle }{G}_{\stackrel{\mathrm{̃}}{c}\stackrel{\mathrm{̃}}{N}r}=b_{\mathrm{1}}{\stackrel{\mathrm{‾}}{N}}_{\mathrm{r}}/\left(a-b_{\mathrm{1}}/\mathrm{2}\right),\end{array}$$

where ${\stackrel{\mathrm{‾}}{N}}_{\mathrm{r}}$ is the mean annual precipitation rate (mm yr⁻¹) estimated from the observed time series as an average over any chosen reference period.

To force the MARCS^HYDRO model, the outputs from global-/regional-scale climate models are needed. The Coupled Model Intercomparison Project 5 (CMIP5; Taylor et al., 2012) is one collection of data sets that is available for climate-scale hydrological studies. At present, the model only needs to be forced by mean precipitation (mm yr⁻¹), evaluated for a future period of 20–30 years. A low index pr indicated that the values were estimated for the future, and ${\stackrel{\mathrm{‾}}{N}}_{\mathrm{pr}}$ is estimated from climate scenarios. Following the assumption that $\stackrel{\mathrm{‾}}{c}$, $G_{\stackrel{\mathrm{̃}}{N}}$ and $G_{\stackrel{\mathrm{̃}}{c}\stackrel{\mathrm{̃}}{N}}$ are constant for both periods, ${\stackrel{\mathrm{‾}}{c}}_r={\stackrel{\mathrm{‾}}{c}}_{\mathrm{pr}}$, ${\mathrm{G}}_{\stackrel{\mathrm{̃}}{N}r}=G_{\stackrel{\mathrm{̃}}{N}\mathrm{pr}}$ and $G_{\stackrel{\mathrm{̃}}{c}\stackrel{\mathrm{̃}}{N}r}=G_{\stackrel{\mathrm{̃}}{c}\stackrel{\mathrm{̃}}{N}\mathrm{pr}}$ (a “basic parametrization scheme” according to Kovalenko, 1993); new parameters of the Pearson equation are calculated from ${\stackrel{\mathrm{‾}}{N}}_{\mathrm{pr}}$:

$$\begin{array}{}\text{(10)}& {\displaystyle }a=\left(G_{\stackrel{\mathrm{̃}}{c}\stackrel{\mathrm{̃}}{N}\mathrm{pr}}+\mathrm{2}{\stackrel{\mathrm{‾}}{N}}_{\mathrm{pr}}\right)/\left(\mathrm{2}{\stackrel{\mathrm{‾}}{c}}_{\mathrm{pr}}\right),\text{(11)}& {\displaystyle }{b}_{\mathrm{0}}=-G_{\stackrel{\mathrm{̃}}{N}\mathrm{pr}}/\left(\mathrm{2}{\stackrel{\mathrm{‾}}{c}}_{\mathrm{pr}}\right),\text{(12)}& {\displaystyle }{b}_{\mathrm{1}}=G_{\stackrel{\mathrm{̃}}{c}\stackrel{\mathrm{̃}}{N}\mathrm{pr}}/{\stackrel{\mathrm{‾}}{c}}_{\mathrm{pr}}.\end{array}$$

Finally, the non-central moments of runoff are calculated for the projected period (denoted by a low index pr):

$$\begin{array}{}\text{(13)}& {\displaystyle }{m}_{\mathrm{1}\mathrm{pr}}=a-b_{\mathrm{1}},\text{(14)}& {\displaystyle }{m}_{\mathrm{2}\mathrm{pr}}=-b_{\mathrm{0}}-\mathrm{2}{m}_{\mathrm{1}\mathrm{pr}}{b}_{\mathrm{1}}+a{m}_{\mathrm{1}\mathrm{pr}},\text{(15)}& {\displaystyle }{m}_{\mathrm{3}\mathrm{pr}}=-\mathrm{2}{m}_{\mathrm{1}\mathrm{pr}}{b}_{\mathrm{0}}-\mathrm{3}{m}_{\mathrm{2}\mathrm{pr}}{b}_{\mathrm{1}}+a{m}_{\mathrm{2}\mathrm{pr}}.\end{array}$$

It should be noted that in core version 0.2 the linear filter model includes the multiplicative stochastic component (see Appendix 1 for details). It may lead to unstable solutions for the Fokker–Planck–Kolmogorov equation (m_k→∞) for statistical moments of high orders. Two methods for getting stable solutions for the Fokker–Planck–Kolmogorov equation are suggested by Kovalenko (2004), and one of them is already implemented in core version 0.1 (Shevnina et al., 2017).

2.4 Model output

In our study, the EPC of runoff was modelled within a Pearson type III distribution. This distribution is commonly used by water engineers to estimate water extremes (Koutrouvelis and Canavos, 1999; Rogdestvenskiy and Chebotarev, 1974; Matalas and Wallis, 1973). The water engineering guidelines provide the ordinates of EPCs from lookup tables (Guidelines, 1984) depending on the CV and CS. These coefficients are calculated from non-central moments' estimates (Rogdestvenskiy and Chebotarev, 1974):

$$\begin{array}{}\text{(16)}& {\displaystyle }\mathrm{CV}=\sqrt{\left(m_{\mathrm{2}}-m_{\mathrm{1}}^{\mathrm{2}}\right)}/m_{\mathrm{1}},\text{(17)}& {\displaystyle }\mathrm{CS}=\left(m_{\mathrm{3}}-\mathrm{3}{m}_{\mathrm{2}}{m}_{\mathrm{1}}+\mathrm{2}{m}_{\mathrm{1}}^{\mathrm{3}}\right)/{\mathrm{CV}}^{\mathrm{3}}{m}_{\mathrm{1}}^{\mathrm{3}}.\end{array}$$

The MARCS^HYDRO model output includes the estimates of the mean values of CV and CS, calculated for the reference period from observations, as well as these estimates simulated from mean precipitation for the projected period. The ordinates of the EPC available from lookup tables then allow the calculation of the runoff values together with their exceedance probability.

3 Model application: a case study

In our study, we chose the basin of the Iijoki River at the Raasakka gauge (lat 25.4^∘, long. 65.3^∘) in order to give an example of the application of the MARCS^HYDRO model on the catchment scale. The Iijoki River is located in north-west Finland, and the Raasakka gauge outlines a watershed area of over 14 191 km². The catchment has a small population, and there are no hydropower plants of multi-year regulation to affect the natural regime of the annual cycle. Thus, one can expect that historical yearly time series of the annual runoff rate do not contain trends connected to artificial regulation. This case study shows an example of the set-up and output of the probabilistic MARCS^HYDRO model.

3.1 The MARCS^HYDRO model set-up: the reference period

The yearly time series of volumetric water discharge of the Iijoki River were extracted from a data set of the GRDC (GRDC 56068 Koblenz, Germany). The observations at the Raasakka gauge (ID = 6 854 600) cover the period 1911–2014, and they do not contain gaps. This period was considered as the reference period. The annual specific water discharge (ARR, mm yr⁻¹) was calculated from the average volumetric water discharge for each year in the reference period. Then, the non-central moments were calculated from the yearly time series of the ARR with the method of moments (see Table 1). The reference climatology (the means of precipitation and air temperature) was evaluated from the data set of NOAA (NOAA/OAR/ESRL PSD, Boulder, Colorado, USA) at a grid node nearest to the watershed centroid (this technique will be discussed in a separate paper, as will the methods of a forcing pre-analysis).

Table 1 The MARCS^HYDRO model set-up: the Iijoki River at Raasakka as a case study.

The values of m_1r, m_2r and m_3r, the moments of runoff, and the mean of precipitation (${\stackrel{\mathrm{‾}}{N}}_{\mathrm{r}})$ were evaluated from observations. The mean air temperature (${\stackrel{\mathrm{‾}}{T}}_{\mathrm{r}}{)}^{\ast }$ was not used in the model set-up in the case of the Iijoki River; however this value allows advancement of the model parametrization (Shevnina et al., 2017).

3.2 The MARCS^HYDRO model forcing: the projected period

Climate scenarios provide a range of projections for temperature and moisture regimes in the future. This range is produced by different assumptions about climate scenarios as well as specific climate models. However, the climate projections include precipitation and air temperature, and they give a forcing to hydrological models in order to simulate projections of runoff. In the case study of the Iijoki River, the data from CMIP5 (Taylor et al., 2012) for three representative concentration pathways (RCPs) were used to force the MARCS^HYDRO model. For each RCP scenario, the projections of annual precipitation rate were applied to test how the MARCS^HYDRO model simulates the EPC under different forcing trajectories. For the period of 2020–2050 (considered the projected period), the mean values of the precipitation rate (mm yr⁻¹) were calculated based on four world-leading global climate models. We used the outputs from the global models CaESM2 (Chylek et al., 2011), HadGEM2-ES (Collins et al., 2011), INM-CM4 (Volodin et al., 2010) and MPI-ESM-LR (Giorgetta et al., 2013) (see Table 2). The mean values of the precipitation rate varied by 2–5 % of the model's average over the RCP scenarios; however, these values alter substantially between the climate models. Among the outputs considered, the MPI-ESM-LR model projects the highest changes in the mean values of the precipitation rate compared to the reference period (see Tables 1 and 2). The HadGEM2-ES model gives the lowest values for the mean values of the precipitation rate. The projected means of the precipitation rate varied slightly between the scenarios. At the same time, they exhibited a significant range of changes among the climate models (the mean values of the precipitation rate ranges from 619 to 737 mm yr⁻¹) for the case of the Iijoki River at Raasakka.

Table 2 The forcing of the MARCS^HYDRO model for the case study of the Iijoki River at Raasakka.

Projected mean air temperature (${\stackrel{\mathrm{‾}}{T}}_{\mathrm{pr}}{)}^{\ast }$ is needed for a regional parametrization scheme (see details in Shevnina, 2011), and these values were not used in the model forcing in the case of the Iijoki River at Raasakka. ${\stackrel{\mathrm{‾}}{N}}_{\mathrm{pr}}$ is the projected mean annual precipitation amount.

3.3 The MARCS^HYDRO model output: the projected period

The projected non-central moments' estimates were simulated for the scenarios and models listed in Table 2. These estimates were used to calculate the mean values of CV and CS (see Eqs. 16–17) that were included in the output of the MARCS^HYDRO model. Table 3 shows the modelling results for the HadGEM2-ES and MPI-ESM-LR global models, for which the water discharges with 10 % and 90 % exceedance probabilities are given. The ordinates of the Pearson type III distribution were extracted from the lookup tables used in hydrological engineering (Druzhinin and Sikan, 2001), and they allow runoff to be expressed as water discharge (m³ s⁻¹). For the Iijoki River at Raasakka, the mean values of ARR and CV vary under the RCP scenarios by over 7 % and 5 % correspondingly. The maximum alteration in the projected mean values of ARR was obtained under RCP8.5 (619 to 737 mm yr⁻¹). Under the projections of the MPI-ESM-LR model, the mean ARR increases by over 17 %.

Table 3 The projected climatology and statistics of annual runoff: the case of the Iijoki River.

In the case of the Iijoki River at Raasakka gauge, the 10 % water discharge exceedance probability will increase in the future under the scenarios and models considered (see Table 3). It may lead to risks of energy spills at hydropower stations located within the catchment of the Iijoki River in the period 2020–2050. At the same time, risks connected with water shortages may be fewer since they are connected to a 90 % water discharge exceedance probability, which is predicted to increase. Figure 2 shows another way in which the model performs the EPC of the annual runoff rate for the Kyrönjoki River at Skatila gauge (GRDC ID: 6854900). The set of EPCs was simulated under three RCP scenarios using a similar set-up to the MARCS^HYDRO model (Shevnina et al., 2019). In the further development of the visualization block, it would be important to involve water managers and decision makers in order to better outline practical applications for the probabilistic hydrological model.

Figure 2 The variability of the tails of the EPCs for annual runoff (DR, mm yr⁻¹) for the reference period (black) and projected period (other colours).

4 Discussions

Nowadays, the future vision of the climate is changing continuously. Climate projections are updated almost every 5–6 years, and many climate models generate meteorological projections for variables such as precipitation and air temperature. Hydrological models are needed to perform an express analysis of future changes in water resources and water extremes (floods and droughts) on a climate scale. The climate scale means that the express analysis is provided for a period of 20–30 years. Lumped or semi-distributed physically based hydrological models are traditionally used on a short-term or seasonal scale to simulate a runoff time series from a time series of meteorological variables (Seibert, 1999). In many catchment-scale hydrological studies, these models are driven by the outputs of climate models or their ensembles in order to evaluate water resources and extremes in the near future (Arheimer and Lindström, 2015; Veijalainen et al., 2012; Yip et al., 2012). The simulation of the runoff time series from a time series of meteorological variables (see Fig. 2 in Veijalainen et al., 2012) leads to high computational costs for such estimations that need to be provided in terms of probability in economical applications (Murphy, 1976). The probabilistic MARCS^HYDRO model is computationally cheaper when compared to lumped or semi-distributed physically based hydrological models. It can easily be coupled with global and regional climate models, and it can provide the express analysis of water resources under a modern version of the future climate.

In this paper we described the structure for the probabilistic hydrological MARCS^HYDRO model, together with the AFA method that lies behind the new model's core version 0.2. The AFA method has a more than 25-year-long history; however, most of the studies are published in Russian (Kovalenko, 1993, 2004, 2009; Kovalenko et al., 2010). The AFA method is based on the statistical theory of automatic systems (Pugachev et al., 1974), which is an outsider among the classical hydrological disciplines. The AFA method is one simplification of the Fokker–Planck–Kolmogorov equation approach that has been developed in the Russian State Hydrometeorological University. It has been tested in many case studies on river basins located in Russia, Colombia, Bolivia and Mali, etc. There are also a number of publications in English (Rosmann and Domínguez, 2017; Shevnina et al., 2017; Kovalenko, 2014; Domínguez and Rivera, 2010; Viktorova and Gromova, 2008). In this paper we formulated the theory logically in an attempt to provide the equations for the new core 0.2 of the MARCS^HYDRO model; however, it also needs to describe the AFA method that lies behind it.

The probabilistic hydrological MARCS^HYDRO model includes the core versions 0.1 and 0.2. In both cores, only three non-central moments are evaluated to construct the EPC within the theoretical distribution the Pearson III type, which is among the traditional distributions of the frequency and risk analysis in hydrology (Kite, 1977; Rogdestvenskiy and Chebotarev, 1974; Sokolovskiy, 1968; Elderton, 1969; Benson, 1968). The model simulates three estimates of non-central moments of runoff instead of a runoff time series, and this circumstance makes the computations by the MARCS^HYDRO model low-cost compared to conceptual hydrological models (Arheimer and Lindström, 2015; Veijalainen et al., 2012). The MARCS^HYDRO model allows the projections of runoff to be put in terms of probability; that is, they appear as runoff values together with their exceedance probability.

The MARCS^HYDRO model includes six modules, and each module allows improvements by including new methods. In this paper, the new model – core version 0.2, extended to simulate the third statistical estimator (skewness) – is presented. The applicability of core version 0.2 is limited by the assumptions behind the AFA approach. Among others, there is the “quasi-stationarity” assumption for the expected climate change. In this case, the climate is described by the statistical estimators (i.e. mean value, variability, etc.) of precipitation, air temperature, evapotranspiration and river runoff, etc., for the period of 20–30 years. It is assumed to consider two time period periods with statistically different climates, namely the reference period and the projected period. Another limitation is connected to the linear filter stochastic model (for details, see Appendix 1) used in core version 0.2. It should be noted that there is a multiplicative component in the model core, and it may lead to unstable solutions of the Fokker–Planck–Kolmogorov equation. Kovalenko (2004) suggests two solutions that result in stable solutions of the Fokker–Planck–Kolmogorov equation. One of the solutions was given by Kovalenko et al. (2010) and is coded in the model core version 0.1 (Shevnina et al., 2017). However, a checking procedure needs to be applied before using this core version. In the checking procedure we plan to use the “beta criterion” method suggested by Kovalenko (2004) to further develop the MARCS^HYDRO model.

Further improvements of the MARCS^HYDRO model are going to be implemented in the block of parametrization and hindcasts. Recently, only the basic parametrization scheme (Kovalenko, 1993) has been included. This basic scheme gives over 70 %–80 % successful hindcasts (forecasts in the past) in the model cross-validation (Shevnina et al., 2017), and the implementation of a regionally oriented parametrization scheme (Shevnina, 2011) is the next step. It needs to include a mean value of the air temperature of the parameter, connected to “noised” watershed physiography in Eq. (A4), the inverse of the runoff coefficient in the work of Kovalenko (1993). It is also important to study the role of the spatial resolution of meteorological forcing in affecting the modelling uncertainties for the simulated mean values of the CV and CS of runoff.

To place the probabilistic MARCS^HYDRO model among other hydrological models, its practical applications need to be better outlined. The model serves as a probabilistic form of long-term hydrological projections, and they require adaptation to the needs of water engineers and water managers as a tool for risk analysis under the expected climate change. The projected EPCs of multi-year river runoff can be applied in designing bridges, pipes and dams, etc., in order to minimize the future risks connected to extreme floods (Shevnina et al., 2017; Kovalenko et al., 2014; Kovalenko, 2009) and to water shortage due to droughts (Viktorova and Gromova, 2008). It is important to define informative forms for the outputs of the MARCS^HYDRO model that can be adapted to the needs of practice, and the development of the block of economic application is among the others studies that are to be continued in close cooperation with water managers and decision makers.

5 Conclusions

The paper describes the theory and assumptions of the AFA approach, as well as the probabilistic hydrological MARCS^HYDRO model's structure and core version 0.2. The features of the model are the close connection to water engineering due to providing the runoff projection in terms of probability, cheapness in terms of computational cost and a wide range of techniques allowing model improvement. In the new core, the third moment linked to the location parameter of the Pearson type III distribution (or asymmetry) was implemented for simulation. In the previous version of the model core, a constant CS∕CV ratio is used to calculate the location parameter of the distribution.

To give a practical example of how to set up the MARCS^HYDRO model, the case of the Iijoki River at Raasakka located in Finland was considered. The model simulated the tailed values of 10 % and 90 % of annual water discharge from the outputs of global climate models. We showed two forms of the probabilistic projections of runoff: an EPC and the runoff values with their exceedance probability. This case study of the Iijoki River at Raasakka shows that the MARCS^HYDRO model gives reasonable results for the meteorological projections considered. The practical applications in water management and decision-making should be clarified in further studies in close co-operation with water engineers.

Appendix A: The theoretical basis for core version 0.2

A1 The assumptions behind advance of frequency analysis

Advance of frequency analysis (AFA) is based on the theory of stochastic systems, specifically, the Fokker–Planck–Kolmogorov equation, which is simplified into a system for three non-central statistical moments (Pugachev et al., 1974). The time series of annual runoff is considered as a realization of a random-process Markov chain type that is assumed to be “stationary”. It means that the statistical estimators (mean, variance and skewness) do not change over the period considered. The statistical estimators are used to model an exceedance probability curve (EPC) of the annual runoff with a Pearson type III distribution. The AFA approach is developed with an assumption of quasi-stationarity (Kovalenko et al., 2010; Kovalenko, 1993). The quasi-stationarity assumption suggests that the statistical estimators of multi-year runoff are different for two periods (the reference period and the projected period). For the reference period, the statistical estimators are evaluated from historical yearly time series of runoff. For the projected period, the statistical estimators of runoff are simulated based on the outputs of global- or regional-scale climate models under any climate scenario.

A2 The linear filter stochastic model

In this context, models replace a complicated hydrological system using maths abstractions and aim to reveal the spatial and temporal river runoff features which are important depending on the goals of study. Among other models, black box hydrological models consider a river basin as a dynamic system with lumped parameters. These models are “based on analysis of concurrent inputs and temporal output series” (WMO-No. 168, 2009) and transform series of meteorological variables (precipitation, air temperature) into series of river runoff. Both input and output series are functions of time (WMO-No. 168, 2009):

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}& a_n\left(t\right){\displaystyle \frac{d^{n}Q}{\mathrm{d}{t}^n}}+a_{n-\mathrm{1}}\left(t\right){\displaystyle \frac{d^{n-\mathrm{1}}Q}{\mathrm{d}{t}^{n-\mathrm{1}}}}+\mathrm{\dots }+a_{\mathrm{1}}\left(t\right){\displaystyle \frac{\mathrm{d}Q}{\mathrm{d}t}}+a_{\mathrm{0}}\left(t\right)Q=\\ & b_n\left(t\right){\displaystyle \frac{d^{n}P}{\mathrm{d}{t}^n}}+b_{n-\mathrm{1}}\left(t\right){\displaystyle \frac{d^{n-\mathrm{1}}P}{\mathrm{d}{t}^{n-\mathrm{1}}}}+\mathrm{\dots }+b_{\mathrm{1}}\left(t\right){\displaystyle \frac{\mathrm{d}P}{\mathrm{d}t}}+b_{\mathrm{0}}\left(t\right)P,\end{array}\end{array}$$

where Q is the runoff in volumetric flow rate, P is the precipitation in volumetric flow rate (rain, snowmelt), and the coefficients a_i and b_i are the empirical parameters of a translating system. These coefficients are the lumped parameters of the black box model. The solution to Eq. (A1) for zero initial conditions gives

$$\begin{array}{}\text{(A2)}& Q\left(t\right)=\underset{\mathrm{0}}{\overset{t}{\int }}h(t,\mathit{\tau })P\left(\mathit{\tau }\right)\mathrm{d}\mathit{\tau },\end{array}$$

where the function h(t,τ) represents the response of a river basin at time t to a single portion of precipitation P(τ) at time (τ). In the AFA approach, a river basin is considered as a linear system, transforming the annual precipitation into the annual runoff:

$$\begin{array}{}\text{(A3)}& a_{\mathrm{1}}\left(t\right){\displaystyle \frac{\mathrm{d}Q}{\mathrm{d}t}}+a_{\mathrm{0}}\left(t\right)Q=b_{\mathrm{0}}\left(t\right)P.\end{array}$$

On the other hand, a river basin can be considered as a linear system with stochastic components in the input function and the model parameter (Kovalenko, 1993):

$$\begin{array}{}\text{(A4)}& \mathrm{d}Q=[-(\stackrel{\mathrm{‾}}{c}+\stackrel{\mathrm{̃}}{c}\left(t\right))Q+(\stackrel{\mathrm{‾}}{N}+\stackrel{\mathrm{̃}}{N}\left(t\right)\left)\right]\mathrm{d}t,\end{array}$$

where $a_{\mathrm{0}}\left(t\right)=\stackrel{\mathrm{‾}}{c}+\stackrel{\mathrm{̃}}{c}\left(t\right)$ is the stochastic parameter of the system (a noised watershed physiography, the inverse of runoff coefficient), $b_{\mathrm{0}}\left(t\right)P=\stackrel{\mathrm{‾}}{N}+\stackrel{\mathrm{̃}}{N}\left(t\right)$ is the stochastic input for the system (a noised precipitation) and a₁=1. The stochastic components of $\stackrel{\mathrm{̃}}{c}\left(t\right)$ and $\stackrel{\mathrm{̃}}{N}\left(t\right)$ are the Gaussian “white noise” with zero means, and their intensities are $G_{\stackrel{\mathrm{̃}}{c}}$ and $G_{\stackrel{\mathrm{̃}}{N}}$. The intensities are mutually correlated as $K_{\stackrel{\mathrm{̃}}{c}\stackrel{\mathrm{̃}}{N}}\left(\mathit{\tau }\right)=E\left(\stackrel{\mathrm{̃}}{c}\right(t\left)\stackrel{\mathrm{̃}}{N}\right(t+\mathit{\tau }\left)\right)=G_{\stackrel{\mathrm{̃}}{c}\stackrel{\mathrm{̃}}{N}}\mathit{\delta }\left(\mathit{\tau }\right)$. It should be noted that the multiplicative parameter $\stackrel{\mathrm{‾}}{c}+\stackrel{\mathrm{̃}}{c}\left(t\right)$ in Eq. (A4) is the sum of the constant $\stackrel{\mathrm{‾}}{c}$ and Gaussian white noise $\stackrel{\mathrm{̃}}{c}\left(t\right)$, and it may lead to unstable solutions of the Fokker–Planck–Kolmogorov equation (i.e. it may lead to infinite statistical moments of high orders). It limits the application of the AFA method (Kovalenko, 1993). Kovalenko (2004) suggests two solutions, and we will introduce them in a further paper.

A3 The Fokker–Planck–Kolmogorov equation and simplifications

The Fokker–Planck–Kolmogorov equation can be applied to simulate the probability density function (PDF) for the stochastic Q(t) in Eq. (A4) (Kovalenko, 1993; Pugachev et al., 1974):

$$\begin{array}{}\text{(A5)}& \begin{array}{rl}{\displaystyle \frac{\partial p(Q,t)}{\partial t}}& =-{\displaystyle \frac{\partial }{\partial Q}}\left(A\right(Q\left)p\right(Q,t\left)\right)+\mathrm{0.5}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial Q^{\mathrm{2}}}}\\ & \left(B\right(Q\left)p\right(Q,t\left)\right),\end{array}\end{array}$$

where p(Q,t) is the PDF of Q at time t; and the drift coefficient (A(Q)) and diffusion coefficients (B(Q)) are calculated as follows (Kovalenko, 1993; Pugachev, 1974):

$$\begin{array}{}\text{(A6)}& {\displaystyle }A\left(Q\right)=-(\stackrel{\mathrm{‾}}{c}-\mathrm{0.5}{G}_{\stackrel{\mathrm{̃}}{c}})Q-\mathrm{0.5}{G}_{\stackrel{\mathrm{̃}}{c}\stackrel{\mathrm{̃}}{N}}+\stackrel{\mathrm{‾}}{N},\text{(A7)}& {\displaystyle }B\left(Q\right)=G_{\stackrel{\mathrm{̃}}{c}}{Q}^{\mathrm{2}}-\mathrm{2}Q{G}_{\stackrel{\mathrm{̃}}{c}\stackrel{\mathrm{̃}}{N}}+G_{\stackrel{\mathrm{̃}}{N}}.\end{array}$$

The analytical solution of Eq. (A5) is difficult and not always needed for practical applications in water engineering since the PDFs of runoff are modelled from a set of statistical estimators, and the moments are from, among others, van Gelder et al. (2006). The PDFs are described with the set of moments $m_k=\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}{Q}^{k}p(Q,t)\mathrm{d}Q$ (where k is the number of the moment, k→∞). To obtain the equations for m_k, both sides of Eq. (A5) were multiplied by a differentiable function ψ(Y) and then integrated within limits from −∞ to +∞ by Q (however, it is supposed that Q>0):

$$\begin{array}{}\text{(A8)}& \begin{array}{rl}& {\displaystyle \frac{d\left(\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}\mathit{\psi }\right(Q\left)p\right(Q,t\left)\mathrm{d}Q\right)}{\mathrm{d}t}}=\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}p(Q,t)A\left(Q\right){\displaystyle \frac{\partial \mathit{\psi }\left(Q\right)}{\partial Q}}\mathrm{d}Q\\ & +\mathrm{0.5}\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}p(Q,t)B\left(Q\right){\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\psi }\left(Q\right)}{\partial Q^{\mathrm{2}}}}\mathrm{d}Q.\end{array}\end{array}$$

Then, ψ(Q) was replaced with ψ(Q)=Q^k, and Eq. (A8) was written as

$$\begin{array}{}\text{(A9)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{m}_k\left(t\right)}{\mathrm{d}t}}=\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}p(Q,t)A\left(Q\right){\displaystyle \frac{\partial \left(Q^k\right)}{\partial Q}}\mathrm{d}Q+\mathrm{0.5}\\ & \underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}p(Q,t)B\left(Q\right){\displaystyle \frac{{\partial }^{\mathrm{2}}\left(Q^k\right)}{\partial Q^{\mathrm{2}}}}\mathrm{d}Q.\end{array}\end{array}$$

For a stationary random process, $\mathrm{d}{m}_k\left(t\right)/\mathrm{d}t=\mathrm{0}$, and the drift and diffusion coefficients are constant. Thus, Eq. (A9) was simplified as follows:

For k=1,

$$\begin{array}{}\text{(A10)}& -(\stackrel{\mathrm{‾}}{c}-\mathrm{0.5}{G}_{\stackrel{\mathrm{̃}}{c}})m_{\mathrm{1}}-\mathrm{0.5}{G}_{\stackrel{\mathrm{̃}}{c}\stackrel{\mathrm{̃}}{N}}+\stackrel{\mathrm{‾}}{N}=\mathrm{0}\end{array}$$

For k≥2,

$$\begin{array}{}\text{(A11)}& \begin{array}{rl}& -k(\stackrel{\mathrm{‾}}{c}-\mathrm{0.5}{\mathrm{kG}}_{\stackrel{\mathrm{̃}}{c}})m_k+k\stackrel{\mathrm{‾}}{N}{m}_{k-\mathrm{1}}-k(k-\mathrm{0.5})G_{\stackrel{\mathrm{̃}}{c}\stackrel{\mathrm{̃}}{N}}{m}_{k-\mathrm{1}}\\ & +\mathrm{0.5}k(k-\mathrm{1})G_{\stackrel{\mathrm{̃}}{N}}{m}_{k-\mathrm{2}}=\mathrm{0}.\end{array}\end{array}$$

Further, the summands in Eqs. (A10)–(A11) were divided by $(\mathrm{2}\stackrel{\mathrm{‾}}{c}+G_{\stackrel{\mathrm{̃}}{c}})$, and new notation was introduced as suggested in the work of Kovalenko (1993) and Pugachev et al. (1974):

$$\begin{array}{rl}& a={\displaystyle \frac{G_{\stackrel{\mathrm{̃}}{c}\stackrel{\mathrm{̃}}{N}}+\mathrm{2}\stackrel{\mathrm{‾}}{N}}{\mathrm{2}\stackrel{\mathrm{‾}}{c}+G_{\stackrel{\mathrm{̃}}{c}}}};b_{\mathrm{0}}=-{\displaystyle \frac{G_{\stackrel{\mathrm{̃}}{N}}}{\mathrm{2}\stackrel{\mathrm{‾}}{c}+G_{\stackrel{\mathrm{̃}}{c}}}};b_{\mathrm{1}}={\displaystyle \frac{\mathrm{2}{G}_{\stackrel{\mathrm{̃}}{c}\stackrel{\mathrm{̃}}{N}}}{\mathrm{2}\stackrel{\mathrm{‾}}{c}+G_{\stackrel{\mathrm{̃}}{c}}}};\\ & b_{\mathrm{2}}=-{\displaystyle \frac{G_{\stackrel{\mathrm{̃}}{c}}}{\mathrm{2}\stackrel{\mathrm{‾}}{c}+G_{\stackrel{\mathrm{̃}}{c}}}}.\end{array}$$

Then, for k= 1, 2, 3, 4 the system of Eqs. (A10)–(A11) includes

$$\begin{array}{}\text{(A12)}& {\displaystyle }& {\displaystyle }{m}_{\mathrm{1}}(\mathrm{2}{b}_{\mathrm{2}}+\mathrm{1})-a+b_{\mathrm{1}}=\mathrm{0},\text{(A13)}& {\displaystyle }& {\displaystyle }(\mathrm{3}{b}_{\mathrm{2}}+\mathrm{1})m_{\mathrm{2}}+(\mathrm{2}{b}_{\mathrm{1}}-a)m_{\mathrm{1}}+b_{\mathrm{0}}=\mathrm{0},\text{(A14)}& {\displaystyle }& {\displaystyle }(\mathrm{4}{b}_{\mathrm{2}}+\mathrm{1})m_{\mathrm{3}}+(\mathrm{3}{b}_{\mathrm{1}}-a)m_{\mathrm{2}}+\mathrm{2}{b}_{\mathrm{0}}{m}_{\mathrm{1}}=\mathrm{0},\text{(A15)}& {\displaystyle }& {\displaystyle }(\mathrm{5}{b}_{\mathrm{2}}+\mathrm{1})m_{\mathrm{4}}+(\mathrm{4}{b}_{\mathrm{1}}-a)m_{\mathrm{3}}+\mathrm{3}{b}_{\mathrm{0}}{m}_{\mathrm{2}}=\mathrm{0}.\end{array}$$

The set of four moments (m₁, m₂, m₃, m₄) is sufficient to model distributions from the Pearson equation (Andreev et al., 2005; Elderton and Johnson, 1969). However, in water engineering we usually only use three-parameter probability distributions fitted to observations (Guidelines, 2004, 1984; Bulletin 17-B, 1982). In this case, $G_{\stackrel{\mathrm{̃}}{c}}\ll \stackrel{\mathrm{‾}}{c}$ is assumed; thus it leads to $b_{\mathrm{2}}=-G_{\stackrel{\mathrm{̃}}{c}}(\mathrm{2}\stackrel{\mathrm{‾}}{c}+G_{\stackrel{\mathrm{̃}}{c}})\approx \mathrm{0}$ and $(\mathrm{4}{b}_{\mathrm{2}}+\mathrm{1})\approx \mathrm{1}$, $(\mathrm{3}{b}_{\mathrm{2}}+\mathrm{1})\approx \mathrm{1}$, $(\mathrm{2}{b}_{\mathrm{2}}+\mathrm{1})\approx \mathrm{1}$. To model the PDFs (or EPCs) of annual river runoff within the Pearson type III distribution, the system of Eqs. (A12)–(A15) is simplified as follows:

$$\begin{array}{}\text{(A16)}& {\displaystyle }& {\displaystyle }-a+b_{\mathrm{1}}=-m_{\mathrm{1}},\text{(A17)}& {\displaystyle }& {\displaystyle }{b}_{\mathrm{0}}+\mathrm{2}{m}_{\mathrm{1}}{b}_{\mathrm{1}}-a{m}_{\mathrm{1}}=-m_{\mathrm{2}},\text{(A18)}& {\displaystyle }& {\displaystyle }\mathrm{2}{m}_{\mathrm{1}}{b}_{\mathrm{0}}+\mathrm{3}{m}_{\mathrm{2}}{b}_{\mathrm{1}}-a{m}_{\mathrm{2}}=-m_{\mathrm{3}}.\end{array}$$

Denoting $\mathit{l}\mathit{k}=\left(\begin{array}{l}-m_{\mathrm{1}}\\ -m_{\mathrm{2}}\\ -m_{\mathrm{3}}\end{array}\right)$, $\mathit{x}=\left(\begin{array}{l}{b}_{\mathrm{1}}\\ b_{\mathrm{0}}\\ a\end{array}\right)$ and $\mathbf{A}=\left(\begin{array}{ccc}\mathrm{1}& \mathrm{0}& -\mathrm{1}\\ \mathrm{2}{m}_{\mathrm{1}}& \mathrm{1}& -m_{\mathrm{1}}\\ \mathrm{3}{m}_{\mathrm{2}}& \mathrm{2}{m}_{\mathrm{1}}& -m_{\mathrm{2}}\end{array}\right)$, the parameters a, b₀ and b₁ are calculated as ${\mathit{x}}_i=D_i/D$, where D is the determinant of matrix A, and D_i is the determinant of the matrix obtained by replacing of the column i (1, 2, 3) in matrix A by the vector lk. Finally, the parameters a, b₀ and b₁ are calculated as follows:

$$\begin{array}{}\text{(A19)}& {\displaystyle }{b}_{\mathrm{1}}& {\displaystyle }=\mathrm{0.5}(\mathrm{3}{m}_{\mathrm{1}}{m}_{\mathrm{2}}-\mathrm{2}{m}_{\mathrm{1}}^{\mathrm{3}}-m_{\mathrm{3}})/(m_{\mathrm{2}}-m_{\mathrm{1}}^{\mathrm{2}}),\text{(A20)}& {\displaystyle }{b}_{\mathrm{0}}& {\displaystyle }=\mathrm{0.5}(m_{\mathrm{1}}^{\mathrm{2}}{m}_{\mathrm{2}}-\mathrm{2}{m}_{\mathrm{2}}^{\mathrm{2}}+m_{\mathrm{1}}{m}_{\mathrm{3}})/(m_{\mathrm{2}}-m_{\mathrm{1}}^{\mathrm{2}}),\text{(A21)}& {\displaystyle }a& {\displaystyle }=\mathrm{0.5}(\mathrm{5}{m}_{\mathrm{1}}{m}_{\mathrm{2}}-\mathrm{4}{m}_{\mathrm{1}}^{\mathrm{3}}-m_{\mathrm{3}})/(m_{\mathrm{2}}-m_{\mathrm{1}}^{\mathrm{2}}).\end{array}$$

A4 Notation

There is too much notation used to describe the model's core version 0.2; thus the secondary parameters of equations were grouped by the model behind it. Table A1 shows the notation and description of the secondary parameters for the linear filter stochastic model. Equation (A3) is a simplification of Eq. (A1) that limits the first-order ordinal differential equation. It includes three parameters, a₀, a₁ and b₀, and two of them are assumed to be noised. These noised parameters include a constant component (indicated with a bar) and a Gaussian white noise component (indicated with a tilde) with their own intensities.

Table A1 The notation and description of the parameters for a linear filter stochastic model.

Table A2 gives a description of the parameters of the Fokker–Planck–Kolmogorov equation and the Pearson system. It should be noted that we do not solve the Fokker–Planck–Kolmogorov equation, and only its simplification for the system of three non-central moments is applied. These non-central moments are estimated from runoff observations for the reference period. For the projected period the moments are calculated from the mean of precipitation.

Table A2 The notation of the Fokker–Planck–Kolmogorov equation and the Pearson equation.

Appendix B: A short user guide for the MARCS model

To set up the model for a single river catchment, the non-central moments should be calculated from historical time series of the annual river runoff rate as well as from a mean value of annual precipitation rate. These values should be placed manually (lines 45–48 in model_core.py, located at https://zenodo.org/record/1220096\#.WyTXxxxRVhw, last access: 4 May 2019) as should the ID number of the catchment (line 51 of model_core.py). To force the model, the projected mean value of the annual precipitation rate should be evaluated from an output of a climate model, and then the model_core.py can be run in the Unix command line: ./model_core.py XXX (where XXX is the mean of the annual precipitation rate for the projected period). The output of model_core.py is stored in the output file model_GPSCH.txt and included in line with the following format: the ID of the catchment, the first non-central moment estimate of the annual runoff rate (mm yr⁻¹) for a reference period, the mean value of the annual precipitation rate (mm yr⁻¹) for a reference period, the coefficient of variation for a reference period, the coefficient of skewness for a reference period, the model parameters $\stackrel{\mathrm{‾}}{c}$, $G_{\stackrel{\mathrm{̃}}{N}}$ and $G_{\stackrel{\mathrm{̃}}{c}\stackrel{\mathrm{̃}}{N}}$, the first non-central moment estimate of annual runoff rate (mm yr⁻¹) for a projected period, the mean value of the annual precipitation rate (mm yr⁻¹) for a projected period, the coefficient of variation for a projected period and the coefficient of skewness for a projected period.