Stochastic Modeling of Hydroclimatic Processes Using Vine Copulas

Abstract

The generation of synthetic time series is important in contemporary water sciences for their wide applicability and ability to model environmental uncertainty. Hydroclimatic variables often exhibit highly skewed distributions, intermittency (that is, alternating dry and wet intervals), and spatial and temporal dependencies that pose a particular challenge to their study. Vine copula models offer an appealing approach to generate synthetic time series because of their ability to preserve any marginal distribution while modeling a variety of probabilistic dependence structures. In this work, we focus on the stochastic modeling of hydroclimatic processes using vine copula models. We provide an approach to model intermittency by coupling Markov chains with vine copula models. Our approach preserves first-order auto- and cross-dependencies (correlation). Moreover, we present a novel framework that is able to model multiple processes simultaneously. This method is based on the coupling of temporal and spatial dependence models through repetitive sampling. The result is a parsimonious and flexible method that can adequately account for temporal and spatial dependencies. Our method is illustrated within the context of a recent reliability assessment of a historical hydraulic structure in central Mexico. Our results show that by ignoring important characteristics of probabilistic dependence that are well captured by our approach, the reliability of the structure could be severely underestimated.

1. Introduction

In the field of hydrology, the study of time series and their synthetic generation is of great importance. They are used to drive models in a wide range of applications, from reservoir design [1,2,3] and planning [4,5,6,7], ecological flow estimation [8], to flood risk [9]. There are a fair amount of models available in the literature and despite potential differences between them, they share two common goals: i) adequate modeling of the marginal distribution and ii) adequate characterization of the auto-correlation structure.

Regarding the modeling of time series of hydroclimatic variables, several challenges can be identified. Firstly, these variables typically present highly skewed behavior, especially in fine time scales [10,11]. In common practice, skewed distributions are simulated using the Thomas–Fiering models [12]. These are linear models that approximate skewness by introducing non-Gaussian white noise in the generation scheme; however, the estimation of white noise distribution parameters requires moments with an order higher than two. The estimates of higher-order moments are very uncertain [13] and have a significant impact on the reproduction of the distribution [14]. Additionally, the linear combination of the i.i.d. non-Gaussian white noise terms does not necessarily result in the target distribution [15], but rather a distribution that has equal moments up the preserved order [16]. Furthermore, hydroclimatic processes exhibit a wide range of dependence structures. Typically, these are categorized as short-range dependence (SRD) and long-range dependence (LRD) [17]. The main characteristic of SRD is the fast decay of the auto-correlation after a few lags, while LRD implies the persistence of high auto-correlation after many lags. It has been shown that LRD is present in many geophysical processes [17,18,19,20,21]. Most applications employ models that are able to reproduce only SRD with some notable exceptions [22,23].

Modeling auto-correlation with linear models (such as those typically used to account for SRD) can have some serious implications. Tsoukalas et al. [16] show that linear models result in unnatural bounded auto-correlation structures when processes other than Gaussian are simulated. This effect was termed envelope behavior. Envelope behavior is caused by the linear combination of white noise described by bounded distributions in the time series generation scheme. For example, when a Pearson-III distribution (which is bounded from below) is used to model the white noise terms in a first-order auto-regressive model (AR(1)), the resulting auto-correlation structure is bounded from below as well [16].

Tsoukalas et al. ([16]) show the importance of shifting the focus from the traditional narrow view of correlation, as expressed by the Pearson correlation coefficient, to the dependence structure in a joint distribution. This is a very important step not only in an effort to recognize the limitations of linear stochastic models, but also to understand better the dependence relationships between hydroclimatic variables. Contemporary research indicates that in many cases this relationship deviates from the Gaussian and the effects on the resulting uncertainty estimates are significant [24,25].

A promising alternative can be found in copulas. Copulas are joint distributions with uniform $[0,1]$ margins. They allow the modeling of the joint distribution through the decoupling of the dependence structure between variables and their marginal distributions. In many cases, this enables the entire description of the dependence structure between variables (or between lagged versions of the same variable) to one measure, such as Spearman’s coefficient. Moreover, copula models are able to reproduce explicitly the marginal distributions.

In comparison to fitting a linear model, copula models can be relatively more complex and time-consuming. In order to fit a variable to a copula function, knowledge about the (physical) behavior of such a variable is required. Additionally, the complexity of these models grows significantly when dimensionality increases. Consequently, this hinders mathematical and computational handling as well as conceptual understanding. Nevertheless, copulas have been widely used and it is a rapidly growing field.

Copulas have gained popularity among the hydrological community [25,26,27,28,29]. Bivariate copulas have been used to generate synthetic time series at the River Nile [30] and the Colorado River [31]. Further, they have been applied to generate multi-site precipitation time series [32]. In the field of hydraulic engineering, for example, [33] applies copulas to produce wind speed and wave height time series for the scheduling of off-shore operations.

Bivariate copulas are limited on the order of auto-correlation they can preserve and the number of variables (for example in terms of geographical location) they can simultaneously simulate (some recent developments allow for approximating non-parametrically complex joint densities. See for example [34,35]). For this reason, extended copula models, i.e., vine copulas, can be used to represent a more complex joint distribution. This method, which was introduced by [36], uses bivariate copulas to decompose multivariate joint distributions. For example, ref. [37] generated cyclostationary stream flow time series considering higher-order auto-dependence. Jäger and Nápoles [38] demonstrated how to produce correlated significant wave heights and mean zero-crossing periods. Sarmiento et al. [39] successfully employed a regular vine (R-vine) model developed by [40] for the simulation of wind speed and direction time series. Vine copulas require the fitting of a large number of copula functions. For example, for the simulation of four variables (taking into consideration only the first-order auto-dependence and their cross-dependence), it would require the calculation of 22 copula parameters.

In this paper, we present a methodology that can be used to generate synthetic time series of hydroclimatic variables. Specifically, we present a case study that concerns an ancient dike that was located within an extinct lake in present-time Mexico City. In order to recreate the hydrological conditions of the lake at the time where the dike was used, time series of precipitation and evaporation are needed. A first attempt to characterize the dependence between the two variables is the model developed in [41]. In this paper, the model presented in [41] is extended to reproduce the first-order auto-dependence structure in the evaporation process. Additionally, specific considerations to better reproduce cross-correlations between the evaporation and precipitation processes are introduced. Furthermore, a novel framework to simulate from an arbitrary number of geographical locations is presented. This approach introduces the coupling of the spatial and temporal models using repetitive sampling. In this way, the spatial model and the temporal model can be freely chosen according to each data set and the required parameters for the quantification of the joint distribution are reduced.

This paper starts by describing the general methodology to generate the time series. Next, the case study is presented including the description of the data. Then, the precipitation–evaporation model is applied to simulate daily realizations for one and four stations. Finally, a discussion of the results, conclusions, and recommendations for future work are presented.

2. Methods

2.1. First-Order Univariate Processes

At the core of the modeling methodology proposed here lies the notion of copulas. Copulas are multivariate distributions with uniform marginals in $[0,1]$. They provide an elegant way to describe the dependence between variables by removing the effects of the marginal distributions. Sklar [42] showed that any multivariate joint distribution can be expressed as a function of its marginals. In a two-dimensional context his theorem states:

$$F_{X_1,X_2}(x_1,x_2)=C\{F_{X_1}\left(x_1\right),F_{X_2}\left(x_2\right)\}$$

(1)

where $X_1$ and $X_2$ are random variables with marginal distributions $F_{X_1}\left(x_1\right)$ and $F_{X_2}\left(x_2\right)$ and C the copula function. Equation (1) indicates that two components are required to define the joint distribution: the copula function and the marginal distributions.

Bivariate copulas have been used before to model time series [43]. The simplest case of temporal dependence regards only two time steps. Let $\left\{X_t\right\}$ for $t\in \mathbb{N}$ denote the time series of interest. In terms of a copula, the relationship between consecutive time steps can be written as $F_{X_t,X_{t-1}}(x_t,x_{t-1})=C\{F\left(x_t\right),F\left(x_{t-1}\right)\}$ where, $F\left(x_t\right)=F_{X_t}\left(x_t\right)$ is the marginal distribution of $X_t$. Moreover, $P(X\le x_t|X_{t-1}=x_{t-1}) =F\left(x_t\right|x_{t-1})=C_{X_t|X_{t-1}}\left\{F\left(x_t\right)\right|F\left(x_{t-1}\right),\mathsf{\Theta }\}$, where $C_{X_t|X_{t-1}}$ is the conditional copula with arguments $\in {[0,1]}^2$. Notice that the parameter vector $\mathsf{\Theta }$ would model auto-correlation of order 1 for the time series of interest. Since $\left\{X_t\right\}$ is assumed to be a stationary process $F\left(x_t\right|x_{t-1})$ corresponds to a first-order Markov model, which is characterized with a single marginal distribution and an appropriate copula. Often, the conditional distribution is expressed as [44]:

$$C_{X_t|X_{t-1}}\left\{F\left(x_t\right)\right|F\left(x_{t-1}\right)\}=\frac{\partial C\left\{F\left(x_t\right),F\left(x_{t-1}\right)\right\}}{\partial F\left(x_{t-1}\right)}$$

(2)

For simplicity, a different notation is introduced based on Aas et al. [45], according to which $h_{X_t,X_{t-1}}\left(F\left(x_t\right),F\left(x_{t-1}\right),\mathsf{\Theta }\right) =\frac{\partial C\left\{F\left(x_t\right),F\left(x_{t-1}\right)\right\}}{\partial F\left(x_{t-1}\right)}$. Moreover, let $h_{X_t,X_{t-1}}^{-1}(F\left(x_{t-1}\right),v,\mathsf{\Theta })$ be the inverse conditional distribution, where the $F\left(x_{t-1}\right)$ is the conditioning variable, v is an independent uniform [0,1] random variable and $\mathsf{\Theta }$ the parameter vector of the bivariate copula $C_{X_t,X_{t-1}}$. Solving for $x_t$ yields the definition of the first-order univariate process.

$$x_t=F^{-1}\left[h_{X_t,X_{t-1}}^{-1}\left(F\left(x_{t-1}\right),v,\mathsf{\Theta }\right)\right]$$

(3)

where $F^{-1}$ is the inverse marginal distribution. The process, briefly described in this section, has been used previously, for example, to model traffic loads for bridge reliability in [46,47].

2.2. First-Order Bivariate Processes

In this section, the univariate model briefly described in Section 2.1 is extended to account for dependence between two variables. Let us consider two dependent processes $\left\{X_t\right\}$ and $\left\{Y_t\right\}$ that present additionally serial correlation. Similarly, as in the univariate case, $t\in \mathbb{N}$. A good alternative to model these kinds of processes are vine copulas or simply vines. Vines were originally introduced in Joe [44], Cooke [36], and Bedford and Cooke [48] (see [49]). Vines are a graphical way of representing multivariate joint distributions. They are a generalization of dependence trees. Roughly, a vine on n elements $V\left(n\right)=\{T_1,\dots ,T_n\}$ is a nested set of trees where the edges of tree j are nodes on the tree $j+1$ for $j=1,\dots ,n-1$. In particular, regular vines are of interest. These are vines whose edges in tree j are connected as nodes in tree $j+1$ only if they share a common node in tree j. Vine copulas as dependence models assign a bivariate copula to each edge in the first tree of the vine and conditional bivariate copulas to the edges of every tree $>1$. For a formal definition and their statistical treatment see for example [45]. A graphical representation of a regular vine on 3 nodes is presented in Figure 1. The edges in Figure 1 are assigned copulas that are parametrized by rank correlation coefficients and a conditional rank correlation of $X_t,Y_t|Y_{t-1}$. Rank correlations are denoted by r followed by a subscript indicating the respective variables and time step. For example, $r_{Y_t{Y}_{t-1}}$ denotes the rank correlation between process Y at times t and $t-1$.

Notice that the essential correlations that the model needs to capture are the temporal dependencies $\left\{X_t\right\}$ and $\left\{Y_t\right\}$ and the cross-correlation $\{X_t,Y_t\}$. If the process $\left\{X_t\right\}$ is generated independently with the univariate conditional copula model discussed in Section 2.1, then all that is left is to induce the correct dependence $\{X_t,Y_t\}$ and $\left\{Y_t\right\}$. This can be achieved by modeling $Y_t$ conditional on both $Y_{t-1}$ and $X_t$ with a vine as the one presented in Figure 1. $Y_t$ can be sampled according to Equation (4) where $F_X$ and $F_Y$ are the marginal distributions underlying $X_t$ and $Y_t$ ($F_X^{-1}$ and $F_Y^{-1}$ their inverses). Notice that the time index t is omitted since both processes are stationary. Moreover, h is the conditional copula (and $h^{-1}$ its inverse) operator, as it was established earlier, and v is as before, an independent uniform random variable on [0,1].

$$y_t=F_Y^{-1}\left[h_{Y_t,Y_{t-1}}^{-1}\left(F_Y\left(y_{t-1}\right),h_{Y_t,X_t|Y_{t-1}}^{-1}(F_X\left(x_t\right),v,{\mathsf{\Theta }}_1),{\mathsf{\Theta }}_2\right)\right]$$

(4)

The modeling procedure for the generation of a continuous bivariate stochastic process may be summarized as:

i

Fit the appropriate marginal distributions characterising the random variables (or use the empirical distribution).

ii

Select suitable copulas $C_{X_t,X_{t-1}}$, $C_{Y_t,Y_{t-1}}$, and $C_{Y_t,X_t|Y_{t-1}}$.

iii

Generate the first variable ($X_t$) using the univariate copula model from the previous subsection.

iv

Generate the second variable ($Y_t$) using Equation (4).

Similar models to the one briefly described in this section have been used in [33,38].

2.3. First-Order Intermittent-Continuous Bivariate Processes

In this section, the algorithm introduced in Section 2.2 is expanded to take into account intermittency. Results are based on the methodology presented in Torres-Alves et al. [41]. The basic idea is the coupling of the copula models with Markov chains. Specifically, the intermittent process is split into two sub-processes, (i) a discrete state-time Markov Chain accounting for the occurrence of a certain environmental condition, and (ii) copula models accounting for “amounts”. For example, in the case of daily precipitation, the Markov chain model would simulate days with or without rain while the copula models would account for the amount of rain (mm) per day. In an intermittent process, a Markov Chain model accounts for two states: one for zero amount (i.e., a dry day) and one for larger than zero amount (i.e., a wet day). The first-order discrete state-time Markov process is described by:

$$P(X_t=s|X_{t-1},X_{t-2},\dots ,X_1) =P(X_t=s|X_{t-1})$$

(5)

Similar to previous subsections, $t\in \mathbb{N}$ while $s\in S$ for a countable set S. The probabilities that express the chance of transition between states are called transition probabilities. Continuing with our example, the two states of precipitation could be $wet\Rightarrow 1$ and $dry\Rightarrow 0$, thus $S=\{0,1\}$. The probabilities of interest are defined as follows:

$$\begin{array}{c}\hfill P(X_t=1|X_{t-1}=1)=P_{1|1}\\ \hfill P(X_t=0|X_{t-1}=1)=P_{0|1}\\ \hfill P(X_t=1|X_{t-1}=0)=P_{1|0}\\ \hfill P(X_t=0|X_{t-1}=0)=P_{0|0}\end{array}$$

(6)

For the above probabilities it holds that $P_{1|1}=1-P_{0|1}$ and $P_{1|0}=1-P_{0|0}$. These probabilities comprise the so-called transition-probability matrix:

$$P=\left[\begin{array}{cc}{P}_{0|0}& P_{1|0}\\ P_{0|1}& P_{1|1}\end{array}\right]$$

(7)

This matrix characterizes the Markov chain, providing all the information necessary to reproduce the process. In this way, the probability of each state is explicitly simulated. Another way to approach intermittency could be by adopting mixed type distributions that consider the dry state as an atom of probability. These have been adopted in literature [50,51] but they will not be discussed in this text; however, they are theoretically compatible with the proposed copula models.

The way that intermittency has been approached, by splitting the main process into two sub-processes of occurrence and amount dictates that the generation algorithm cannot be continuous. Specifically, the generation of the ”occurrence” process must precede the “amount” process. In other words, the Markov chain must precede the copula time series.

Consequently, to generate a synthetic time series, the first step is to generate realizations of the Markov chain. This will define the different state blocks (i.e., sequences of wet and dry periods of time) of the time series. For the dry (zero amount state) block of the intermittent process, no further action is needed. What follows is the generation of the amount process for the wet blocks. This can be performed for each block by employing the univariate model of Equation (3). Then the continuous process can be generated conditionally on the intermittent one of the same time step using the vine copula model of Equation (4); however, dry blocks require special treatment. Naturally, zero amount does not allow the expression of dependence via a copula. To overcome this, the cross-correlation is treated implicitly by conditionalizing the continuous process distribution according to the two states, wet and dry.

On a final note, to preserve correctly the first-order auto-dependence structure of the continuous process the generation of blocks must be sequential. This means that the first realization of each block should be conditional on the final realization of the previous block. The generation procedure is schematized in Figure 2. Superscripts W and D denote the conditional wet and dry continuous process marginal distributions, respectively.

The generation algorithm is summarized below:

i

Fit the appropriate marginal distributions characterizing the RV’s.

ii

Select and fit the suitable copulas $C_{X_t,X_{t-1}}$, $C_{Y_t,Y_{t-1}}$, and $C_{Y_t,X_t|Y_{t-1}}$.

iii

Calculate the transition-probability matrix of the Markov chain for the intermittent process $X_t$. Simulate a desired length of the Markov sequence.

iv

Split the time series into dry and wet blocks.

v

Generate “blocks” representing $\left\{X_t\right\}$ using the univariate copula model (Equation (3)).

vi

Generate the wet block $Y_t$ using the vine model of Equation (4) and the appropriate marginal distribution. Use as seed for the first value of the wet block the last value of the previous dry block.

vii

Generate the dry block $Y_t$ using the vine model (Equation (3)) and the appropriate marginal distribution. Use as seed for the first value of the dry block the last value of the previous wet block.

2.4. Multivariate Processes

Equations (3) and (4) define univariate and bivariate models, respectively; however, using the vine decomposition, higher-order models can be can be defined. As it was discussed during the introduction, while such models that account for temporal and spatial dependence exist, they require the calculation of a great number of parameters and are limited in the dependence structures they can represent. The methodology described herein aims to approximate dependent processes in time and space (for example precipitation and evaporation from an arbitrary number of stations) while remaining parsimonious and flexible. The key of this methodology is to approximate the conditional distribution of the spatial and temporal dependence through repetitive sampling instead of inferring it theoretically. It should be noted that the use of repetitive sampling in the context of hydrology has been used in [52,53] for the coupling of different temporal scales.

For simplicity, let us first assume a two-dimensional case. Let us define a two-dimensional process $\{Y_t^1,Y_t^2\}$ where the superscripts denote, for example, a corresponding station and the subscripts the time step. Moreover, let us assume that the temporal and spatial dependence can be described by the first-order model of Equation (3). The methodology can be described in four steps:

i

Generate one realization of $\left\{Y_t^1\right\}$ from $Y_{t-1}^1$ based on the temporal model.

ii

Generate n realizations of $Y_t^2$ based on $Y_t^1$ according to the spatial model. These are denoted with a tilde and the superscript $S2$ as $[{\tilde{Y}}_{t,1}^{S2},\cdots ,{\tilde{Y}}_{t,n}^{S2}]$.

iii

Generate n realizations of $Y_t^2$ based on $Y_{t-1}^2$ according to the temporal model. These are denoted with a tilde and the superscript $T2$ as $[{\tilde{Y}}_{t,1}^{T2},\cdots ,{\tilde{Y}}_{t,n}^{T2}]$. The spatial and temporal realizations are all plausible realizations of $Y_t^2$, which implies that a section between the two sets exists.

iv

Identify the common space. This is performed by identifying the realizations that minimize the root mean squared error (RMSE).

$$\underset{i}{min}\{\sqrt{{\left({\tilde{Y}}_{t,i}^{T2}-{\tilde{Y}}_{t,i}^{S2}\right)}^2}\}$$

(8)

Notice that RMSE in step iv above is not used as a goodness of fit measure but rather as a selection criterion between the spatial and temporal components of the model. Another way to approach this could be a Monte Carlo procedure where a convergence criterion is targeted; however, in the context of time series generation this would result in non-practical computational demand. Other measures for selection may be very well possible; however, investigating them is out of the scope of this paper since we focus mainly on the general methodology for simulating spatially diverse data sources.

The above procedure can easily be extended to more dimensions. Let us consider another case where a set of four correlated variables $[Y_t^1,Y_t^2,Y_t^3,Y_t^4]$ is generated. It is reminded that the vine representation in four variables is not unique [54]. Two very common vines are the canonical or a C-Vine and the D-Vine. Their graphical representations are found in Figure 3.

It can be seen that the C-Vine representation has a unique variable, in our case $Y_t^1$, from which the other variables depend. On the contrary in the D-Vine representation, the variables depend serially. Both could be probable representations of spatial dependence; however, the algorithm depends on the vine’s sampling order. Let us assume that a C-Vine is selected. Further, let $v_1,v_2,v_3,v_4$ denote independent random draws from the uniform $(0,1)$ distribution. The sampling algorithm is given by the following equations:

$$y_t^1=F_Y^{-1}\left(v_1\right)$$

$$y_t^2=F_Y^{-1}\left[h_{Y_t^2,Y_t^1}^{-1}\left(F_Y\left(y_t^1\right),v_2,{\mathsf{\Theta }}_1\right)\right]$$

(9)

$$y_t^3=F_Y^{-1}\left[h_{Y_t^3,Y_t^1}^{-1}\left(F_Y\left(y_t^1\right),h_{Y_t^2,Y_t^3|Y_t^1}^{-1}\left(F_Y\left(y_t^2\right),v_3,{\mathsf{\Theta }}_2\right),{\mathsf{\Theta }}_3\right)\right]$$

$$y_t^4=F_Y^{-1}\left\{h_{Y_t^4,Y_t^1}^{-1}\left[F_Y\left(y_t^1\right),h_{Y_t^2,Y_t^4|Y_t^1}^{-1}\left(F_Y\left(y_t^2\right),h_{Y_t^3,Y_t^4|Y_t^1,Y_t^2}^{-1}\left(F_Y\left(y_t^3\right),v_4,{\mathsf{\Theta }}_4\right),{\mathsf{\Theta }}_5\right),{\mathsf{\Theta }}_6\right]\right\}$$

Furthermore, let us assume that the temporal dependence can be described by the first-order model of Equation (3). The first step is to generate a basis realization on which the repetitive sampling will be based. This decision depends on the selected spatial dependence model. It is noted that according to Equation (9), the C-Vine sampling order is 1-2-3-4. Thus, in this case, the algorithm should start by computing $Y_t^1$. If a D-Vine was selected instead, the algorithm should start by computing $Y_t^3$ since the sampling order of a D-Vine is 3-2-4-1 [54]. The rest of the steps remain the same. n realizations are generated according to the spatial and temporal models of which the one that minimizes the RMSE (Equation (10)) is selected.

$$\underset{i}{min}\left\{\sqrt{\frac{{\sum }_{j=1}^m{\left({\tilde{Y}}_{t,i}^{Tj}-{\tilde{Y}}_{t,i}^{Sj}\right)}^2}{m}}\right\}$$

(10)

where $m=4$ in our case. This procedure is schematized in Figure 4.

The methodology above can be extended for m stations. The required number of parameters of such a model is:

$$N_{\mathsf{\Theta }}=m(m+1)/2$$

(11)

The multivariate generation algorithm is summarized below:

i

Fit the appropriate marginal distributions characterizing the RV’s.

ii

Identify suitable models to describe spatial and temporal dependence.

iii

Select and fit $N_{\mathsf{\Theta }}$ suitable copulas to model dependencies.

iv

Select the number of trials n for the repetitive sampling.

v

Generate the first temporal realization according to the sampling order of the selected spatial model.

vi

Generate n realizations from the spatial and temporal models independently.

vii

Select the realization which minimizes the error (Equation (10)).

2.5. Admissible Marginal Distributions and Copula Fitting

One of the major benefits of modeling stochastic time series with copulas is the direct use of the marginal distribution, instead of moment approximations. This allows the utilization of classical fitting techniques such as maximum likelihood estimation [55] and L-Moments [56] and state-of-the-art methods such as K-Moments [14] and the metastatistical extreme value framework (MEV) [57]. Moreover, any finite variance distribution is admissible, which increases the flexibility of the models.

The same flexibility can be demonstrated in the admissible copula fitting techniques. In this study, the estimation of copula parameters is performed via the Spearman or Rank correlation coefficient accompanied by measures to test the goodness of fit, such as the Cramer–Von Mises statistic [58] and semi-correlations [43]; however, other methods such as maximum likelihood, Bayesian information criteria, or Kendall’s tau can be supported. A more comprehensive review can be found in [25,59].

3. Case Study

To demonstrate the effectiveness of the proposed vine copula model, two scenarios are presented: (i) simulation of daily evaporation and precipitation at one station and (ii) simulation of daily evaporation at four stations, all located at the Valley of Mexico. Additionally, the first scenario illustrates the influence of the dependence structure on the reliability of flood defenses.

The proposed methodology is validated on the accuracy of the reproduced historical characteristics, such as the approximation of the historical marginal distributions, the first-order auto-correlation dependence, and the cross-correlation dependence structure.

3.1. Area of Study

The case study described herein is presented in Torres-Alves & Morales-Nápoles (2020) [41]. It concerns the reliability analysis due to overflow of the ancient Nezahualcoyotl dike that was once located in the Valley of Mexico (Figure 5). The Aztecs built this structure around 1450 AD to protect their capital, the city of Tenochtitlan, from rising water levels in the lacustrine system.

Around 1519, the Valley of Mexico was covered by a lacustrine system comprising of six interconnected lakes. These lakes connected during high water and because of the salinity of the lakes, the rising water posed a threat to agriculture and freshwater supplies in Tenochtitlan. Nowadays, there are no remains of this dike and the lakes were almost completely drained by the end of the 19th century. Ref. [41] reconstructed the geometry and position of the ancient structure based on historical sources. The dike was 16 km long with a height of 8 m and a width of 3.5 m. It was made out of wood, stone, and mud. Moreover, ref. [41] estimated the extent of the lacustrine system and proposed a simple hydrological balance equation (Equation (12)) to compute the water level fluctuation at the lake. Finally, an assessment of the reliability of the dike due to overflow was presented.

$$X_t{A}_{tr}-Y_t{A}_l=\frac{d{V}_t}{dt}$$

(12)

where:

• $X_t$ is the precipitation at time t (m/day).
• $Y_t$ is the evaporation at time t (m/day).
• $A_{tr}$ is the tributary area of the basin (m${}^2$).
• $A_l$ is the surface area of the lake (m${}^2$).
• $\frac{dV}{dt}$ is the daily change in volume (m${}^3$/day).

The hydrological balance was simulated only for the wet seasons. For this reason, an initial water level ($L_0$) needed to be assumed. The dike’s probability of failure due to overflow was computed for six initial water levels (1 to 6 m).

3.2. Data

To simulate the precipitation and evaporation time series, records from five stations situated within the Valley of Mexico are analyzed. This analysis was originally performed in [41]. Therefore, only a brief overview is provided in the present work. For a detailed presentation of the historical records the reader is referred to the original publication. The data are available from Mexico’s national database [60]. Information about the stations is given in

Table 1. Summary of precipitation and evaporation information for each station.

Code	Name	Coordinates		Time Period
		Longitude	Latitude
15008	Atenco	−98.92	19.54	1 January 1961–22 August 2013
15044	La Grande	−98.92	19.58	1 January 1964–31 August 2014
15083	San Andres	−98.92	19.53	1 May 1967–31 December 2014
15167	El Tejocote	−98.92	19.44	1 December 1957–1 January 2007
15170	Chapingo (DGE)	−98.9	19.48	13 March 1957–1 January 2000

.

In this region, two seasons are identified, (i) wet season (May–October), and (ii) dry season (November–April). The data from the stations were analyzed and appropriate distributions were fitted using Maximum likelihood (

Table 2. Selected distributions for all variables and stations.

Station	Probability Distribution Function
	Precipitation	Evaporation
Atenco	Weibull	Generalized Extreme Value
La Grande	Weibull	Generalized Extreme Value
San Andres	Weibull	Generalized Extreme Value
El Tejocote	Weibull	Generalized Extreme Value
Chapingo (DGE)	Weibull	Generalized Extreme Value

). The respective cumulative distribution functions are presented in Appendix A.

Furthermore, bivariate copulas were fitted for each data pair (precipitation–evaporation). The results are presented in

Table 3. Selected copula pairs for all stations.

Pairs	Atenco	La Grande	San Andres	El Tejocote	Chapingo (DGE)
$X_t-X_{t-1}$	Gauss	Gauss	Gauss	Gauss	Gauss
$Y_t-Y_{t-1}$	Gumbel	Gumbel	Gumbel	Gumbel	Gumbel
$X_t-Y_t$	Gumbel	Gumbel	Gumbel	Gumbel	Gumbel

. In [41] the authors use different goodness of fit (GOF) techniques to show that the pairs in Table 3 are adequate choices for the (conditional) copulas of interest. For further details about the GOF techniques used, the reader is referred to [41] and references therein. Their functional representations are given in Appendix B.

3.3. Simulation of Daily Evaporation and Precipitation

The simulation of precipitation $\left\{X_t\right\}$ and evaporation $\left\{Y_t\right\}$ poses a particular challenge due to the intermittent behavior exhibited by precipitation. Moreover, these variables are the basic drivers for many hydrological models (i.e., [61]). Producing correlated time series enables the quantification of uncertainty in these drivers and a probabilistic description of the results. For this scenario, the data from Atenco station were analyzed (for further details on why this station is selected, refer to [41]). From these data, 5000 hypothetical wet seasons were generated, each with a length of 185 days.

The ability of the model to reproduce the historical distributions is exhibited in Figure 6 with the use of QQ plots. These depict the relationship between the theoretical and empirical quantiles of the cumulative distribution function. In the presented case, the plotted quantiles form an almost straight line; thus, good agreement is found between historical and simulated CDFs for both variables. Figure 7 provides a comparison between historical and synthetic dependence structure in the real domain. Precipitation and evaporation exhibit a moderate correlation in the range of 0.3 to 0.4. Additionally, between the variables, a low negative correlation is observed. These dependencies are quantified in terms of rank correlations in

Table 4. Comparison of historical and simulated rank correlations.

	${\mathit{X}}_{\mathit{t}}-{\mathit{X}}_{\mathit{t}-1}$	${\mathit{Y}}_{\mathit{t}}-{\mathit{Y}}_{\mathit{t}-1}$	${\mathit{X}}_{\mathit{t}}-{\mathit{Y}}_{\mathit{t}}$
Historical	0.27	0.4	−0.22
Synthetic	0.27	0.37	−0.18

. Overall, the model reproduces adequately the first-order auto-dependence as well as the dependence between variables.

The Effect of the Choice of Copula

For the description of the historical first-order auto-dependence, a Gumbel copula was employed; however, in non-copula models addressing directly the characteristics of the observed dependence structure is not straightforward. To demonstrate the importance of a correct dependence description, a reliability calculation will be performed using two sets of the synthetic time series. The first one will be the one generated during the previous section. The second one is a precipitation and evaporation pair with the auto-dependence of evaporation described by a Gaussian copula. The latter is selected for comparison purposes because it represents a simple dependence structure that is reproduced by commonly used linear Gaussian models and more recent models such as [50,51,53,62]. The results for both cases are summarized in

Table 5. Probability of failure for different initial water levels ($L_o$) as calculated with the Gaussian and Gumbel time series.

Water Level at the Foot of the Dike (m) (${\mathit{L}}_{\mathit{o}}$)	m.a.s.l.	Gaussian Copula		Gumbel Copula
		${\mathit{P}}_{\mathit{f}}$	Return Period (Years)	${\mathit{P}}_{\mathit{f}}$	Return Period (Years)
1	2231	0.00171	586	0.00099	1006
2	2232	0.02641	38	0.02068	48
3	2233	0.16189	6	0.14956	7
4	2234	0.45131	2	0.44657	2
5	2235	0.74568	1	0.75102	1
6	2236	0.91468	1	0.92038	1

.

Given an initial water level of $L_0$ = 1 m. the probability of failure computed for the Gaussian copula case is $73\%$ higher than the Gumbel case.

Note that the Gaussian copula reproduces the same first-order rank auto-correlation as the Gumbel. The difference can be entirely attributed to the change in the shape of the dependence structure (Figure 8). In the absence of tail dependence, there is no clustering of high evaporation events. This moderates the outgoing flux of the lake resulting in more frequent high-water levels and thus, a larger probability of failure. This can be illustrated by comparing the two evaporation time series (Figure 9).

In Figure 9, the time series generated by the Gumbel copula presents clusters of high evaporation values around time steps 40 to 70, 100 and 160. On the contrary, the time series generated by the Gaussian copula demonstrate a more symmetrical dependence pattern. This is in accordance with the scatter plots provided in Figure 8 where the difference of tail dependence between both copulas is demonstrated. Additionally, it is observed that the differences between the computed probabilities decrease as the initial water level $L_0$ increases (Table 5). This is because larger probabilities (smaller return periods) correspond to events further from the tail, where the differences between the copula functions are found. These results provide evidence on the importance of selecting an appropriate model that characterizes the dependence structure of extreme events.

3.4. Simulation of Evaporation across Multiple Stations

This section presents the application of the multivariate methodology. The developed multivariate algorithm is implemented for the simulation of daily evaporation in four of the five stations at the Valley of Mexico. Evaporation data from Atenco $\left\{Y_t^1\right\}$, La Grande $\left\{Y_t^2\right\}$, San Andres $\left\{Y_t^3\right\}$, and Chapingo (DGE) $\left\{Y_t^4\right\}$ were used to generate 100 seasons (each of 185 days). The spatial dependence is described by a Gaussian copula while the temporal dependence is described by a Gumbel copula. For the latter, only the first-order auto-dependence is taken into account. The sampling number is set to $n=100$. A comparison between historical and simulated first-order auto-dependence is provided in Figure 10.

On the temporal level, evaporation in all stations demonstrates a rank correlation in the range of 0.35 to 0.45. Moreover, the historical data exhibit upper tail dependence between consecutive time steps. In all cases, the simulated data reproduce the shape of the historical dependence structure. The results are quantified in

Table 6. Comparison of historical and simulated rank auto-correlations.

	${\mathit{Y}}_{\mathit{t}}^1-{\mathit{Y}}_{\mathit{t}-1}^1$	${\mathit{Y}}_{\mathit{t}}^2-{\mathit{Y}}_{\mathit{t}-1}^2$	${\mathit{Y}}_{\mathit{t}}^3-{\mathit{Y}}_{\mathit{t}-1}^3$	${\mathit{Y}}_{\mathit{t}}^4-{\mathit{Y}}_{\mathit{t}-1}^4$
Historical	0.41	0.36	0.45	0.44
Synthetic	0.41	0.37	0.46	0.40

where a very good agreement is found between historical and simulated rank correlations.

Furthermore, evaporation exhibits a high spatial correlation, in the order of $r\approx 0.6$ among the stations. From a physical point of view, this can be explained by the proximity of the stations and their position within the same watershed. The proposed algorithm was able to reproduce these dependencies to a good extent (Figure 11); however, a small deviation is observed between pairs $Y_t^1-Y_t^4$ and $Y_t^2-Y_t^3$ (

Table 7. Comparison of historical and simulated rank cross-correlations.

	${\mathit{Y}}_{\mathit{t}}^1-{\mathit{Y}}_{\mathit{t}}^2$	${\mathit{Y}}_{\mathit{t}}^1-{\mathit{Y}}_{\mathit{t}}^3$	${\mathit{Y}}_{\mathit{t}}^1-{\mathit{Y}}_{\mathit{t}}^4$	${\mathit{Y}}_{\mathit{t}}^2-{\mathit{Y}}_{\mathit{t}}^3$	${\mathit{Y}}_{\mathit{t}}^2-{\mathit{Y}}_{\mathit{t}}^4$	${\mathit{Y}}_{\mathit{t}}^3-{\mathit{Y}}_{\mathit{t}}^4$
Historical	0.57	0.70	0.70	0.50	0.48	0.63
Synthetic	0.60	0.70	0.62	0.56	0.50	0.63

). This difference can be attributed to the copula choice.

Finally, the model’s ability to reproduce the marginal distributions is presented in Figure 12 where the empirical and synthetic quantiles of the cumulative distribution function are compared. Overall, a good agreement is observed for all variables; however, small differences are found near the tails. The size of this disagreement is not atypical in hydrological research or other fields where uncertainty in modeling or observations is, in general, significant.

Exploration of the Repetitive Sampling Length

As with any repetitive sampling approach, the accuracy of the method depends upon the sampling length. To determine the relationship between accuracy and sampling length, four simulations were performed with $n=[25,50,100,250]$ and the empirical cumulative distributions of their RMSE (Appendix B) were computed. In Figure 13, the relationship between RMSE and sampling length for 90%, 95%, and 99% percentiles is plotted. As the sampling length increases, the RMSE decreases. The largest decrease can be found between $n=25$ and $n=50$. For $n=100$, the RMSE of the 99% percentile is equal to 0.1. The previous simulation suggests that this level of error yields adequate results. Furthermore, increasing the sampling length does not offer significant improvements in accuracy while at the same time it drastically increases the computational cost. From this perspective, a sampling length between 50 to 100 is a reasonable choice.

4. Conclusions

This paper presents a methodology to simulate hydroclimatic variables through copula-based models. The proposed methodologies focus on the use of vine copulas to characterize complex temporal and spatial probabilistic dependence. Two cases are presented: (i) the generation of correlated precipitation and evaporation on a daily scale and (ii) the generation of correlated daily evaporation time series from four stations. The first case consists of a trivariate vine copula to handle temporal dependence of evaporation and cross-dependence with precipitation. To capture intermittency, Markov chains were coupled with a copula-based model. In the second case, the methodology relied on repetitive sampling to couple the temporal and spatial vine copula models. The result is a parsimonious and flexible model capable of simulating accurately observations from an arbitrary number of stations.

In both case studies discussed herein, the models proved their capability to capture the underlying dependence structures, reproduce the marginal distributions, and intermittency. In fact, for the dependence structure, important asymmetries (such as tail dependence) may be incorporated. These asymmetries are often found but overlooked in the modeling of hydroclimatic variables. One of our case studies shows that by ignoring these asymmetries, the reliability of hydraulic structures, for example, could be underestimated by a factor of 2 in some cases (Table 5).

We have shown the ability of our approach combining Markov chains and vine copula models to reproduce short-range dependence (SRD) structures [63]. Many hydroclimatic processes however exhibit long-range dependence (LRD) [64,65,66]. The performance of our approach for these kind of dependence may be a subject for future research.