2.1 Observational data:
Daily precipitation records from 86 meteorological stations over SESA were considered during the period 1979-2017 (Figure 1a), which have been used and quality-controlled in previous studies by the research group (Olmo and Bettolli, 2021b). These records were provided by the National Weather Services of Argentina, Brazil, Paraguay and Uruguay and the National Institute of Agricultural Technology of Argentina and presented less than 20% of missing data.
Based on the evaluation of multiple cross-validated ESD models performed by Olmo and Bettolli (2021b), the present experiment selected specific models from this previous work following the perfect-prognosis approach, which considers reanalysis predictors as quasi-observations to train the statistical models (Maraun et al. 2015). ESD models were trained during the observational period 1979-2017 using ERA-Interim daily mean fields (Dee et al. 2011) of the following variables interpolated at 2° grid resolution over 50-67°W and 20-40°S (blue box in Figure 1a): geopotential height at 500 and 1000 hPa (Z500 and Z1000, respectively); meridional wind at 850 hPa (V850); specific humidity at 700 and 850 hPa (Q700 and Q850, respectively) and air temperature at 700 and 850 hPa (T700 and T850, respectively). The choice of predictors was proposed considering previous downscaling experiments and the evaluation of the synoptic environment associated with extreme precipitation over SESA (Rasmussen and Houze 2016; Bettolli and Penalba 2018; Olmo and Bettolli 2021a). Thus, the predictor sets used in this work involved different configurations with pointwise predictors - using the values from the four and sixteen nearest grid points to the target station point - and spatial-wide predictors information by using the principal components (PCs) - explaining 95% of the total variance (Table 1). The first approach was considered to address local influences, whereas the last one was used to retain the large-scale atmospheric structures. Although a detailed description of the application of the ESD models will be given in the following lines, the reader is referred to Olmo and Bettolli (2021b) for a comprehensive evaluation of the calibrated models selected for the present work.
Simulations from 12 GCMs from the CMIP5 and CMIP6 experiments (Taylor et al. 2012; Eyring et al. 2016) were used in the application of the ESD models. These sets of models were selected since they have full availability for the different predictor variables and levels, including the Andes mountain range. Historical outputs of the predictor variables mentioned above during the period 1979-2005 (1979-2014) were considered for the CMIP5 (CMIP6) models. Furthermore, in order to obtain downscaled climate projections for the 21st century, future model outputs covering the period 2006-2100 (2015-2100) from the RCP8.5 (SSP585) simulations were employed. Although these two future scenarios are not strictly comparable, both imply the worst-case scenario in each CMIP experiment. GCMs predictor variables were interpolated to a common grid of 2°. Additionally, raw precipitation outputs from the GCMs in their native grid resolution - selecting the closest model grid cell to each station point - were included in the different analyses for comparative purposes. A complete description of the GCMs used in this study is presented in Table 2.
In order to produce a large multi-model ensemble of statistically downscaled GCMs, different families of ESD models were considered in the present work:
The analog method (AN, Zorita and von Storch 1999) consists of finding, for each daily record, the most similar large-scale atmospheric configuration (the nearest neighbour) in the rest of the period to obtain the local prediction according to a similarity metric (in this case, the Euclidean distance). The ANs have been widely applied in different regions due to its simplicity and ability to account for the non-linearity of the predictor-predictand relationships (Timbal et al. 2009; Bettolli and Penalba 2018; Gutiérrez et al. 2019). The main drawback of this technique is that the catalog of analog days is constrained to the observed period and therefore it cannot predict values outside the observed range. This makes the ANs sensible to possible non-stationarities in a climate change scenario, which should be cautiously considered (Benestad, 2010).
- Generalized Linear Models
The generalized linear models (GLMs) are an extension of linear regression in which the predictand variable can belong to a non-linear distribution. Given that in this work special attention is taken in extreme precipitation, stochastic versions of the GLM were adopted, considering a random sample from the predicted distribution at each time step (GLM_STs). As performed by Olmo and Bettolli (2021b), these models consisted in a two-stage implementation, with a Bernoulli distribution and a logit link for precipitation occurrence, and a Gamma distribution and log link for precipitation amounts (San Martin et al. 2017; Bedia et al. 2020). Hence, two simulated time-series were produced at each station point so that the final rainfall time-series was calculated by multiplying these two series (see Bedia et al. (2020) for further details).
- Generalized Linear Models + Weather Types
This family of ESD models is made up of circulation-conditioned GLMs (GLM_WTs) (San Martin et al. 2017; Olmo and Bettolli 2021a). The classification of 16 weather types (WT), found in Olmo and Bettolli (2021a) based on Z500 anomalies over southern South America, was used to calibrate the GLM individually for each WT. Note that these are deterministic GLMs, thus, the precipitation occurrence is obtained by transforming the simulated probability values into binary ones considering the climatological probability of rain as threshold, whereas precipitation amounts are the expected values modelled by optimizing the conditional mean. For this method, the Z500 daily anomalies of the GCMs were projected into the ERA-Interim weather types according to the Euclidean distance.
The artificial neural networks (NNs) are non-linear regression-based models made of a chain of neurons organized in layers of feed-forward networks. The model structure relies on these neurons being connected between consecutive layers from the input layer - daily predictors fields - through a set of hidden layers in the network, to the output layer. These connections are characterized by different weights - that the network optimize from the input data - and activation functions, which are non-linear functions applied in each neuron (Cavazos, 1999; Haylock et al., 2006; Vu et al. 2015; Baño-Medina et al., 2020; Baño-Medina et al., 2021). Following the model construction employed in Olmo and Bettolli (2021b), a two-layers NN configuration with 25 and 15 neurons, respectively, was selected. The learning rate parameter was set to 0.01 and the activation function in the hidden layers was the sigmoid function. Similarly to the GLMs, a stochastic set-up of NNs (NN_ST) was included in this study, which was implemented in two stages: one for precipitation occurrence - for which the output layer performed a linear function - and one for precipitation amounts, for which the output layer performed the sigmoid function. Following the recommendations by Baño-Medina et al. (2020), the NN optimizes the negative log-likelihood of a Bernoulli-Gamma distribution. Specifically, it estimates the rain probability from the Bernoulli distribution and the shape and scale parameters of the Gamma distribution to obtain the rain occurrences and amounts, respectively.
Note that deterministic GLM and NN simulations of daily precipitation were previously analyzed, resulting in systematic underestimations of precipitation amounts - particularly in extremes - and generally presented poorer performances than their stochastic versions or other ESD families. Hence, this work will be centered on the application of the ESD techniques that best represented the main features of daily precipitation as shown by Olmo and Bettolli (2021b). For each ESD model family (AN, GLM_ST, GLM_WT and NN_ST), two models with different predictors configuration were chosen based on their performance in the cross-validation procedure (Olmo and Bettolli 2021b). In this way, in the present study 8 ESD models - outlined in Table 1 - were applied to the GCMs outputs for their historical and future periods. The ESD methods were implemented using the R-based climate4R open framework (Bedia et al. 2020).
c) Evaluation framework
The perfect prognosis approach assumes that GCMs are able to adequately reproduce the predictors as depicted by the reanalysis dataset. For this reason, the distributional similarity between the representation of predictor variables by the GCMs and ERA-Interim was addressed here by performing a two-sample Kolmogorov-Smirnov (KS) test, which is a non-parametric test for checking the null hypothesis that two datasets come from the same distribution. The KS statistic presents values from 0 to 1, where the lowest values indicate greater distributional similarity. Following the methodology suggested by Bedia et al. (2020), standardized anomalies of the predictor variables were compared using an effective sample size due to the strong serial correlation on the daily time-series, at the 99% level of confidence.
Regarding the precipitation simulations, statistically downscaled GCMs were evaluated separately according to the CMIP version of each GCM, that is, model analysis was carried out for their respective historic and future periods individually (Table 2). The meteorological station points over SESA were used as reference. On one hand, multiple validation indices adapted from the VALUE intercomparison experiment (Maraun et al. 2015) were estimated during the historical period, encompassing the main features of daily precipitation. The indices are the relative bias (Bias Rel.), the frequency of rainy days (R01), the precipitation intensity (as depicted by the simple precipitation intensity index, SDII), the intensity and frequency of days with precipitation above 20 mm (R20 and R20p, respectively) and the 98th percentile of the daily precipitation distribution in wet days (P98). In all cases (observations and models), the threshold for a rainy day was set at 1 mm per day. On the other hand, the spatio-temporal variability of precipitation extremes was assessed by analyzing their spatial behaviour, intensity and intra-annual variability. To this end, extreme rainfall was defined at each station point as the precipitation values in those days when the accumulated precipitation exceeded the 95th percentile (P95) of the empirical distribution of rainy days. This daily percentile was calculated in the base period 1986-2005, considering a 29-days moving window centered on each calendar day. In addition, Taylor diagrams (Taylor, 2001) were illustrated for summarizing the representation of the spatial patterns of precipitation over SESA. These diagrams quantify the degree of statistical similarity between the stations reference dataset and the models, reporting the Pearson correlation coefficient, the standard deviation and the centered root mean squared error.
For the analysis of the downscaled projections, anomalies from the 1986-2005 base period were estimated in each model output and time-series of specific precipitation indices were constructed using both historical and future simulations. The indices used here were chosen to study changes in both the intensity and frequency of rainy days (PPmean and PPfrequency, respectively) and precipitation extremes based on the estimation of the P95 as described above (R95 and R95p for extremes intensity and frequency, respectively). Moreover, agreement maps showing the coincidence among the multi-model ensemble in the projected changes were plotted over SESA.
In order to inter-compare the statistically downscaled projections obtained here with other dynamical downscaling models, the regional climate models RegCM4v7 and REMO2015 from the CORDEX-CORE initiative (Gutowski et al. 2016) were considered. Particularly, daily RCMs precipitation from the historical (1979-2005) and RCP8.5 (2006-2100) experiments driven by three CMIP5 GCMs were employed (MPI-ESM-LR, MPI-ESM-MR and NorESM1-M). Note that these three GCMs were also used for the ESD simulations in the present work (Table 2). In this way, 4 simulations were compared with the different statistical downscaling outputs of the common GCMs: MPI-ESM-MR_RegCM4v7, NorESM1-M_RegCM4v7, MPI-ESM-LR_REMO2015 and NorESM1-M_REMO2015. A more detailed description of the RCMs simulations is displayed in Table 3. Even though the validation of these RCMs is beyond the scope of this study, they have already been evaluated and used over the region in the literature and their skills and shortcomings were identified (Olmo and Bettolli 2021a; Teichmann et al. 2020). Thereby, the focus of this analysis will be on the comparison of the variety of downscaled projections in terms of model agreement and spread in the intensity of changes over SESA.