High-Resolution Precipitation Gridded Dataset on the South-Central Zone (34° S–41° S) of Chile

1. Introduction

Meteorological spatial data consistent with observational information are critical for several scientific fields: environmental, hydrology, agriculture, application of renewable energy, biology, economy, and sociology, among others (Parra et al., 2004; Hijmans et al., 2005; Abatzoglou, 2013; Cannon et al., 2015; Liu et al., 2017; Sun et al., 2018). Within weather variables, precipitation is the driver of the hydrological cycle and the most difficult to estimate (Michaelides et al., 2009; Kidd and Levizzani, 2011; Tapiador et al., 2012; Beck et al., 2017b). Precipitation datasets provide information for hydrological models such as SWAT or TOPMODEL (Berezowski et al., 2016) and environmental verification (Ji et al., 2015; Berezowski et al., 2016; Brinckmann et al., 2016; Fick and Hijmans, 2017).

Interpolation methods use information from local weather networks, atmospheric reanalysis products, weather radar, satellite data, or a combination of these products to construct global or local gridded rainfall datasets because each one has its negative sides, in particular precipitation. For instance, datasets based only on instrumental observation networks depend on the spatiotemporal distribution, quality, and length of records (Sun et al., 2014). Atmospheric reanalyses, on the other hand, are at low resolution and sometimes inadequately parameterize sub-grid processes, leading to misrepresentation of synoptic-scale and convective dynamics (Roads, 2003; Ebert et al., 2007; Kidd et al., 2013). Another example is weather radar, since it provides data at high temporal and spatial resolution, but of limited spatial coverage (Koistinen, 1991; Kitchen and Blackall, 1992; Chen et al., 2008; Beck et al., 2017a). Satellite data allow for direct coverage of large areas $\left(\sim {0.25}^{\circ }\right)$ (Sorooshian et al., 2014; Ashouri et al., 2015) but are partly insensitive to light rainfall or drizzle events, sensitive to systematic error, and also unfitted for snow- and ice-covered surfaces on mountain areas (Ferraro et al., 1998; Ebert et al., 2007; Kidd and Levizzani, 2011; Kidd et al., 2012; Laviola et al., 2013; Beck et al., 2017a).

Some of these sources for obtaining precipitation grids (satellite, reanalysis, atmospheric modeling, or statistical downscaling techniques) may present problems in areas of complex topography where high spatial-temporal heterogeneity is difficult to estimate, as for example Chile (Daly et al., 1994; Chen et al., 2014; Herold et al., 2016; Beck et al., 2017a; Zambrano M. et al., 2017). Several studies have assessed the insufficiency of these datasets (e.g., Nastos et al., 2016; Beck et al., 2017b; Camera et al., 2017; Liu et al., 2017; Hu et al., 2018; Sun et al., 2018; Timmermans et al., 2019) including some zones in Chile (Muñoz et al., 2011; Ward et al., 2011). Just recently, this problem has been addressed at the national level in Chile (Zambrano F. et al., 2017; Zambrano M. et al., 2017).

High-resolution datasets properly covering Chile are scarce. This is surprising, given that the country is among those most affected by climate change (Kreft et al., 2016), although the complex topography and multiplicity of climates make development of these datasets difficult. Besides, precipitation coverage is inadequate and inhomogeneous, limiting water resources studies (Hosseini-Moghari et al., 2018). To date, datasets are available only in specific regions by local weather stations (Román and Andrés, 2010; Zambrano, 2011; Jacquin and Soto-Sandoval, 2013; Reyes, 2013; Castro et al., 2014; Sijinaldo, 2015; Cifuentes, 2017) or from low-resolution gridded datasets (Morales-Salinas et al., 2012).

Despite these difficulties, there are known methods that allow for building high-resolution databases in an efficient way. Among statistical downscaling techniques for spatial interpolation, regression-based methods are computationally cheap to run, easy to understand, and statistically efficient and straightforward (Semenov et al., 1998; Dibike and Coulibaly, 2005; Hashmi et al., 2011; Pahlavan et al., 2018). Thus, interpolation methods based on functions describing spatial changes of the target variable (Brinckmann et al., 2016) generate continuous spatial data (e.g., Hay et al., 1998; Marquínez et al., 2003; Ceccherini et al., 2015; Beck et al., 2017b).

The main goal of this study was to improve, quantify, and offer a new basis for average monthly precipitation at high spatial resolution $\left(\sim 800\text{\hspace{0.17em}m}\right)$ using a dynamic-statistical method. We focus this analysis on the central-south zone of Chile (Figure 1) in the period $2000\text{–}2011$. We structured this article as follows: in Section 2, we describe the regional climate and topography. Section 3 includes description of the dynamic topo-climatic model and the datasets used. This is followed by Section 4, which reports the major findings. The article ends with a discussion and conclusion on the robustness of our model and main findings.

FIGURE 1. (A) Study area in central-southern Chile (O’Higgins–Los Ríos). (B) River basins used for spatial yield analysis of the weather research and forecasting model, the gridded precipitation sets (Table 2), and the topo-climatic model. The red dots indicate the 136 stations used as input for the topo-climatic model and the 60 green dots indicate the stations used for performance analysis.

2. Study Area

Variability of the South Pacific subtropical anticyclone and high-latitude pressure centers dominates Chile’s synoptic climate (Fuenzalida, 1982; Rojas, 2016). The narrow width (180 km on average) between the Pacific coast and the Andes (https://www.gob.cl/nuestro-pais/), along with a north–south extension of 4,000 km, contributes to a variety of climates (Valdés et al., 2016a). Chile holds a distinct climatic characteristic from low to high moisture from north to south. Because of the orographic effect, the main precipitation pattern shows an annual accumulation moving to the south and the Andes (Pizarro et al., 2012; Quintana and Aceituno, 2012; Valdés-Pineda et al., 2014; Valdés et al., 2016a), with annual precipitations in the north 30°S of about 100 mm increasing up to 3,000 mm in the south 40°S (Pizarro et al., 2012; Barrett and Hameed, 2017). Since elevation increases from the sea level in the west to some thousand meters in a few hundred kilometers to the east, Andean precipitation can double coastal rainfall (Viale and Garreaud, 2015; Barrett and Hameed, 2017).

Within Chile, our study area corresponds to the center-south zone of the country between latitudes $\left({34}^{\circ }\text{S–}{41}^{\circ }\text{S}\right)$, from O’Higgins to Los Ríos regions (Figure 1). This region presents a well-defined annual cycle characterized by a peak of precipitation in austral winter and lower values in austral summer (Valdés et al., 2016a; Valdés et al., 2016b). Besides, this is one area with the highest humidity in the southern hemisphere (Miller, 1976) and one of the extra-tropical areas most affected by El Niño and La Niña events (Grimm et al., 2000; Waylen and Poveda, 2002). Due to its transitional location between subtropical and high-latitude temperatures, climate is characterized by a significant interannual rainfall variability (Aceituno et al., 1993; Quintana, 2004; Le Quesne et al., 2009). The El Niño Southern Oscillation, mostly in the northern section, and the Pacific Decadal Oscillation to the south affect hydrological regimes (Arumí-Ribera and Oyarzún-Lucero, 2006; Araya-Ojeda and Isla, 2016). The south-central zone presents yearly accumulated precipitation between 100 and 3000 mm from the north to south (Quintana and Aceituno, 2012; Valdés et al., 2016a). This is a zone of active recharge for surface and groundwaters as a result of infiltration and transport processes in valleys, channels, and fractured systems in the high zones and the pre-mountain ranges (Arumí et al., 2012; Carling et al., 2012). Along the Andes, particularly toward the south, where the role of oceanic fronts is prominent, the orographic effect increases rainfall.

3. Methodology

The choice and organization of observed precipitation data are presented in Section 3.1. Section 3.2 shows the development and execution of the atmospheric model, and then the adjustment and selection (which considers the database presented in Section 3.1 of the modeled precipitation fields). The collection of global gridded precipitation sets and their later data choice are discussed in Section 3.3. Section 3.4 presents the structure and procedure in the construction of the topo-climatic dynamic model. The statistical analysis strategy required for the process of calibration and validation of the atmospheric model, choice of global gridded data and topo-climatic dynamic model, and final construction of the high-resolution gridded product is addressed in Section 3.5.

3.1. Instrumental Database

The Dirección General de Aguas and the Dirección Meteorológica de Chile (DMC) support Chile’s meteorological station network. Although these data are reliable and cover several decades, many stations are conventional (report only daily accumulation) and at low elevations (Garreaud et al., 2016). Also, in most records (Figure 2), data gaps are common (Zambrano M. et al., 2017b). Within the study zone, we initially considered 301 stations reporting daily observations (Figure 1). We eliminated stations with more than $30\text{\%}$ monthly data gaps between 2000 and 2011 (105 stations). From the remaining stations (196), 136 correspond to input for model calibration and 60 for validation.

FIGURE 2. Number of precipitation stations in Chile used by Global Precipitation Climatology Centre dataset (Schamm et al., 2015). Source: Zambrano M. et al. (2017).

3.2. Atmospheric Modeling

An atmospheric numerical simulation was developed using the Weather Research and Forecasting model (WRF), version 3.6 (Skamarock et al., 2008) for the period $2000\text{–}2011$ to obtain a high-resolution robust precipitation dataset. The simulation was initialized in October 1999 and run through December 2014, with the first 3 months discarded as model spin-up. As recommended by the DMC (DMC, 2015), this model must consider a domain covering some parts of the Pacific Ocean to capture the phenomena approaching the continent (Yáñez-Morroni et al., 2018). We used three nested domains for downscaling (Figure 3). The largest domain, $D_1$, at 36 km resolution, contained two nested domains, $D_2$ and $D_3$, each with a ratio of 1/3 relative to its respective parent (Table 1). Initial and boundary conditions for the $D_1$ domain correspond to the ERA-Interim global reanalysis, at a spatial grid resolution of ${0.75}^{\circ }\times {0.75}^{\circ }$ latitude–longitude and 60 vertical hybrid levels (Dee et al., 2011). Ingestion of ERA-Interim’s atmospheric fields was every 6 hours. The twelve-year simulation took about 22 real computational months in the Troquil cluster at the Center of Excellence of Scientific Modeling and Computing at Universidad de La Frontera. In atmospheric numerical simulations, parameterizations strongly influence the patterns of simulated precipitation (Wang and Seaman, 1997; Gallus and William, 1999; Jankov et al., 2005) and also affect model response to changes in grid spacing (Gallus and William, 1999) or soil moisture (Gallus and William, 1999; Gallus et al., 2000). A summary of the settings used is shown in Table 1, selected following recommendations of the DMC (DMC, 2015).

FIGURE 3. Display of the domains used in weather research and forecasting model. D₁: parent domain with 36-km resolution. D₂: first nesting with 12-km resolution (1/3 parent domain size). D₃: second nesting, with 4-km resolution (1/9 parent domain size).

TABLE 1. Domain configuration and parameterization schemes selected for the weather research and forecasting model.

We implemented a method for systematic error correction (bias) applied to every WRF’s generated precipitation field within each regional domain (Figure 1A). We calculated the relationship between the average monthly precipitation of a station and its nearest grid point of the WRF model (a virtual station). For each regional domain, Eq. 1 determined selection of WRF’s virtual stations. Next, we obtained and applied a fit factor $\mu \left(P_o^s\right)/\mu \left(P_{wrf}^s\right)$ to each grid point in the region (Eq. 1). Finally, each adjusted grid point was reevaluated and discarded as input to the dynamic topo-climatic model if it did not meet the requirement of Eq. 1.

where $P_{min}^o$ and $P_{max}^o$ are the minimum and maximum observed rainfall, respectively; $\mu \left(P_o^s\right)$ is the monthly observed average of stations; $\mu \left(P_{wrf}^s\right)$ is the monthly average of the nearest virtual stations delivered by the WRF model; $P_{wrf}$ is the computed gridbox; and $P_{wrf}^{cs}$ is the bias-corrected precipitation.

3.3. Global Gridded Datasets

Our compilation includes eight precipitation products from diverse sources (observation, reanalysis, satellite, or mix) and present different spatial and temporal resolutions, and coverage (semiglobal or global, Table 2). We did not bias-correct these precipitation products as this was performed by the teams that built them. Therefore, each product’s dynamics depends entirely on its source. Three precipitation sets were based on satellite data ($\text{TRMM}$, $\text{CMOPRH}$, and $\text{SM}2\text{RAIN}$), one set was based exclusively on observed data ($\text{CRUv}4.01$), and four sets were built with satellite data, reanalysis, and observations ($\text{PGFv}3$, $\text{CHIRPSv}2.0$, $\text{PERSIANN}-\text{CDR}$, and $\text{MSWEPv}1.2$).

TABLE 2. Description of the 8 (quasi) global (sub-)precipitation grids used in this study.

In order to minimize the difference in bias from these precipitation products (Zhang et al., 2019; Yeh et al., 2020), each gridded precipitation product used in this study was evaluated according to Eq. 2, where every point was discarded when it was either under or over $5\text{\%}$ of the respective minimum and maximum precipitation observed within the respective regional domain (Figure 1A). This procedure tries to cover areas where atmospheric modeling precipitation grids are discarded (see Section 3.2). Thus, once the best set for each regional domain was obtained, no overlap of data with the WRF model was sought. Otherwise, that is, if overlap happens, the one with the smallest difference in minimum and maximum precipitation bias is selected according to Eq. 2. Note that if neither atmospheric nor the gridded data met the admissibility conditions, that place is not used.

where $P_g$ is the precipitation gridbox of each set in Table 2, $P_g^s$ is the selected gridbox of the grid set, and $P_{min}^o$ and $P_{max}^o$ are the minimum and maximum precipitation of each subregion, respectively. Otherwise, we delete this gridbox.

3.4. Dynamic Topo-Climatic Model

Our model is a multiple linear regression between precipitation as the dependent variable (Swain et al., 2017; Navid and Niloy, 2018; Devi et al., 2020) and the following independent variables: elevation, slope, exposure, continentality, latitude, and longitude (Zambrano, 2011; Camera et al., 2014; Cifuentes, 2017). The model fits the observed values by a least square procedure between the observed and predicted values (Delbari et al., 2019). The linear dependence function is given by

where α is the intercept, and the values ${\beta }_i,i=1,\mathrm{.},p$ are called regression coefficients associated with the variables elevation $\left(\text{elev}\right)$, exposure $\left(\text{exp}\right)$, slope $\left(\text{slp}\right)$, continentality $\left(\text{cont}\right)$, latitude $\left(\text{lat}\right)$, and longitude $\left(\text{lon}\right)$.

To get the dataset gridded at 800 m, we used the Shuttle Radar Topography Mission digital elevation model (Farr et al., 2007) at 90-m spatial resolution. This digital elevation model is aggregated to an 800-m resolution (Figure 4). Note that the model domain covers five administrative regions: O’Higgins, El Maule, Biobío (now Biobío and Ñuble), La Araucanía, and Los Ríos (Figure 1A). However, the domains correspond to regional boundaries (Figure 1A), along with a $25\text{\%}$ increase in the contribution from the neighboring regions. With this, we ensured that basins shared between different regions are covered. Also, mutual covering between domains ensures greater consistency and homogeneity to the global solution. Each month, data used for the construction of the field change with input data (Sections 3.2 and 3.3).

FIGURE 4. Digital terrain model with independent variables: (A) elevation, (B) slope, (C) exposure, and (D) continentality.

For each regional domain, we solve the topo-climatic model monthly, trying different configurations of the independent variables described in Eq. 3. We carry out tests in which either a single topographic variable or a set of them is eliminated. Thus, each month, six $\left(6!\right)$, different tests can be performed. The structure is tested monthly, and at the end of the 12-year period (144 months), the most effective model is used as the basis for the calculation of the monthly topo-climatic model coefficient. Because it plays a significant role in a topo-climatic model, elevation is the only variable used in all calculations. Similarly, the different tests indicated that the elevation and exposure variables are necessary for all regional domains. In the following, only the most representative models are shown.

Each model developed is labeled according to the variables in Eq. 3: $\left(all\right)$, all the variables; $\left(wexp\right)$, without exposure; $\left(wslp\right)$, without slope; $\left(wlat\right)$, without latitude; $\left(wlon\right)$, without longitude; and $\left(wcon\right)$, without continentality variable. Subsequently, each model was evaluated through an analysis of variance (ANOVA) test and a residual analysis. We use different statistics to determine the most appropriate structure for each regional domain, namely, the standardized standard deviation (${\sigma }_e$, Eq. 6), standardized mean square error ($nmrse$, Eq. 7), percentage of the standardized systematic error ($pbias$, Eq. 9), correlation coefficient (${\rho }_{xy}$, Eq. 8), Nash–Sutcliffe coefficient ($nse$, Eq. 10), and the modified Willmott coefficient (d, Eq. 11).

Figure 5 depicts the three phases carried out for the construction of the dynamic topo-climatic model. The first phase $\left(F1\right)$ determines the number of local precipitation stations and their subsequent selection (Section 3.1). The second phase $\left(F2\right)$ corresponds to climatic modeling that produces precipitation fields, which are then bias-corrected (Section 3.2) using the local data (Section 3.1). The third phase $\left(F3\right)$ (Section 3.3) selects the global precipitation gridded sets which are used alternatively for the locally adjusted precipitation modeled fields. Finally, within the fourth phase $\left(F4\right)$ (Section 3.4), on each geographical domain, the tests allow the choice of the coefficients of variables (elevation, slope, exposure, continentality, latitude, and longitude) to develop the high-resolution fields.

FIGURE 5. Schematic overview of the dynamic topo-climatic model development.

3.5. Statistical Analysis

This section presents the statistical analysis strategy required for the atmospheric model correction and adjusting process, global gridded dataset selection, the topo-dynamic model calibration and validation, and the final high-resolution gridded product construction.

Spatial fitting, selection of monthly fields from the WRF model, and spatial selection of the best monthly products of global grid precipitation are analyzed for each regional domain (Figure 1A): The results of these processes are displayed for eight hydrological basins (Figure 1B). Our strategy uses a range of statistics following Cardoso et al. (2013), Ji et al. (2015), and Akhter et al. (2019): standardized standard deviation (${\sigma }_e$, Eq. 6), standardized mean square error ($nrmse$, Eq. 7), correlation coefficient (${\rho }_{xy}$, Eq. 8), percentage of the standardized systematic error ($pbias$, Eq. 9), Nash–Sutcliffe coefficient ($nse$, Eq. 10), and modified Willmott coefficient (d, Eq. 11). We used six indices in the data comparison because each of them quantifies only one aspect of data agreement or disagreement.

Once the selection of the best grid points is completed, a different monthly topo-climatic model is constructed over each regional domain. For this process, a calibration phase is followed by a validation phase. As previously indicated, we used $70\text{\%}$ of the local rainfall data $\left(136\text{\hspace{0.17em}stations}\right)$ for calibration and $30\text{\%}$$\left(60\text{\hspace{0.17em}stations}\right)$ for validation. In both processes, due to the stations’ inhomogeneous distribution and notorious inadequate coverage in mountainous areas, the spatial distribution of the local stations chosen tried to maximize homogeneous coverage of each regional domain.

Calibration of the dynamic topo-climatic model includes application of an ANOVA test, and the goodness-of-fit parameters $R^2$ (Eq. 4) and $R_j^2$ (Eq. 5) with the respective residual analysis. For the validation stage, the following statistics are used: standardized standard deviation (${\sigma }_e$, Eq. 6), standardized mean square error ($nrmse$, Eq. 7), correlation coefficient (${\rho }_{xy}$, Eq. 8), percentage of the standardized systematic error ($pbias$, Eq. 9), Nash–Sutcliffe coefficient ($nse$, Eq. 10), and the modified Willmott coefficient (d, Eq. 11). We also analyze results using Taylor’s diagrams (Taylor, 2001), to visualize results model behavior relative to observations in terms of the spatial structure generated by the Pearson correlation coefficient $\left({\rho }_{xy}\right)$, standardized mean square error ($nrmse$), and absolute standardized standard deviation $\left(\left|{\sigma }_e\right|\right)$. Also, we used boxplot diagrams to visualize standard deviation scores $\left({\sigma }_e\right)$, the mean systematic error percentage $\left(pbias\right)$, coefficient of Nash–Sutcliffe ($nse$), and modified Willmott coefficient (d). This calibration and validation processes are performed monthly, thus generating every month a different topo-climate model. Calibration and validation results presented here correspond to the average of all these months. Also, given that the whole domain includes a climatic transition zone, the results of the eight basins are summarized in three subdivisions: northern, central, and southern areas (Figure 1B).

The statistics used are defined in the following equations:

where X and Y are observed and simulated precipitation, respectively.

4. Results

In Sections 4.1 and 4.2, the average precipitation behavior of the adjusted atmospheric model (see Section 3.2) and the selected global gridded products (see Section 3.3) is presented over three zones (north, center, and south), covering eight basins (Figure 1B). All local stations (196) are used, and these are presented using Taylor and boxplot diagrams (see Section 3.5). Section 4.3 presents the best results associated with the calibration process of the dynamic topographic-climatic model within each regional domain (Section 4.3.1) using 132 precipitation stations (Section 3.1). The model validation stage of the chosen model using 60 stations in Section 4.3.2 is presented for the eight basins (Figure 1B). This part of the article ends by showing the annual accumulation of the constructed monthly climatology.

4.1. Atmospheric Model

In the northern zone (Figures 6A-C), ${\rho }_{xy}$ shows an average of 0.85, but we can appreciate atypical cases with values below 0.65. Meanwhile, $nmrse$ shows agreements close to 0.65, although all basins present cases above 0.5. For ${\sigma }_e$ (Figure 7A), the average agreement is close to 0.1, with Rapel showing some cases above 0.3. Finally, $pbias$ (Figure 7B) gives an average agreement of $14\text{\%}$, with Mataquito presenting the lowest values since $75\text{\%}$ of the analyzed stations there are close to $10\text{\%}$.

FIGURE 6. Taylor’s diagrams corresponding to weather research and forecasting model performance (Table 1) once fitted (Eq. 1) in basins: (A) Rapel, (B) Mataquito, (C) Maule, (D) Itata, (E) Biobío, (F) Imperial, (G) Toltén, and (H) Valdivia.

FIGURE 7. Box diagrams corresponding to the fit of the weather research and forecasting model, once adjusted (Eq. 1) on the river basins of the area under study (Figure 1B). (A) Standardized standard deviation $\left({\sigma }_e\right)$, (B) percentage of the systematic error $\left(pbias\right)$, (C) Nash–Sutcliffe coefficient $\left(nse\right)$, and (D) Willmott coefficient $\left(d\right)$.

Results for the central zone (6D,E) indicate an average for ${\rho }_{xy}$ of 0.83 and 0.55 for $nrmse$, although with several cases well above 0.65. For ${\sigma }_e$ (Figure 7A), average is 0.14, with the Biobío basin showing values above 0.25. For $pbias$ (Figure 7B), the $17\text{\%}$ average is punctuated by certain basins exceeding $25\text{\%}$.

Good estimation for ${\rho }_{xy}$ characterizes the southern region (Figures 6F-G), averaging a value of 0.87. Within this region, Valdivia shows cases close to 0.6. $nrsme$ exhibits an average of 0.56, although all basins have cases above 0.7. In turn, ${\sigma }_e$ (Figure 7A) shows an average of 0.15, with Imperial and Toltén above 0.3. Finally, $pbias$ (Figure 7B) shows an average of $9\text{\%}$. Toltén gives the lowest average estimate with $7\text{\%}$.

We show the $nse$ and Willmott coefficient $\left(d\right)$ results in all regions (and basins) in Figure 7D. In all basins, the average value for $nse$ is greater than 0.7; only Rapel and Maule show values smaller than 0.6. The Willmott coefficient (d) shows values close to 0.8 in the north zone (Rapel, Mataquito, and Maule) and a value of 0.75 for the central zone (Itata and Biobío), whereas for the southern zone (Imperial, Toltén, and Valdivia) $25\text{\%}$ of stations show values below 0.75 with a large spread.

4.2. Global Precipitation Datasets

In the northern zone (Figures 8A-C), ${\rho }_{xy}$ presents an average value of 0.74. However, $\text{CMORPH}$ and $\text{SM}2\text{RAIN}$ show ${\rho }_{xy}$ values under 0.6 in all zones. Variations in $nrmse$ show an average of 0.4, with Mataquito and Maule displaying 0.5 for the $\text{PERSIAN}$ dataset. For ${\sigma }_e$, the average is 0.27, while for $\text{CMORPH}$ and $\text{SM}2\text{RAIN}$ the values are higher than 0.6 for the whole zone (Figure 9A). $pbias$ (Figure 9B) shows an average of $20\text{\%}$, with anomalous estimates (over $30\text{\%}$) in Mataquito and Maule.

FIGURE 8. Taylor’s diagrams corresponding to the gridded set performance (Table 2), once fitted (Eq. 2) on the river basins of (A) Rapel, (B) Maule, (C) Itata, (D) Biobío, (E) Imperial, (F) Toltén, (G) Mataquito, and (H) Valdivia.

FIGURE 9. Box diagrams corresponding to the gridded set performance (Table 2), once adjusted (Eq. 2) on the river basins of the area under study (Figure 1B). (A) Standardized standard deviation $\left({\sigma }_e\right)$, (B) percentage of the systematic error $\left(pbias\right)$, (C) Nash–Sutcliffe coefficient $\left(nse\right)$, and (D) Willmott coefficient $\left(d\right)$.

Over the central zone (Figures 8D,E), ${\rho }_{xy}$ presents an average of 0.74. Again, $\text{CMORPH}$ and $\text{SM}2\text{RAIN}$ show the worst results for ${\rho }_{xy}$, with scores under 0.6. Meanwhile, $nrsme$ gives an average of 0.44, although the entire area shows values above 0.3. Continuing with ${\sigma }_e$ (Figure 9A), the average is 0.34. As above, $\text{CMORPH}$ and $\text{SM}2\text{RAIN}$ give estimates fairly different relative to the average, higher than 0.5. Here, $pbias$ (Figure 9B) shows an average error of $20\text{\%}$, while the highest estimates are below $30\text{\%}$ in the Itata and Biobío basins.

For the southern zone (Figures 8F-H), ${\rho }_{xy}$ shows an average of 0.73. Similar to that described in the previous paragraph, $\text{CMORPH}$ and $\text{SM}2\text{RAIN}$ show the worst results with estimates for ${\rho }_{xy}$ under 0.6. A high $nrmse$ in all basins delivers an average of 0.48. $\text{MSWEP}$ and $\text{PERSIAN}$ show the highest $nrmse$ with values above 0.6. The average ${\sigma }_e$ is 0.34, with $\text{CMORPH}$ and $\text{SM}2\text{RAIN}$ above 0.6 (Figure 9A). Meanwhile, $pbias$ (Figure 9B), exhibits an average of $20\text{\%}$. In this case, it is important to highlight that in all basins, the different datasets present errors more significant than $30\text{\%}$.

The Nash–Sutcliffe ($nse$) and Willmott (d) coefficients are displayed for all basins (Figures 9C,D). Here, $nse$ shows a mean adjustment of about 0.6 from the Rapel to the Imperial basin. There is a larger spread for Toltén and Valdivia, with the minima below 0.55. The Willmott coefficient (d) exhibits a similar distribution from Rapel to Imperial, with almost $90\text{\%}$ of the datasets above 0.7. We again see more heterogeneity in the Toltén and Valdivia basins with values below 0.5. $\text{SM}2\text{RAIN}$ shows the weakest performance for all basins.

4.3. Dynamic Topo-Climatic Precipitation

4.3.1. Regional Domain Analysis

The average fit of the dynamic topo-climatic model for each regional domain is shown in Figure 10. The coefficient of determination $\left(R^2\right)$ (Figure 10A) shows an excellent model performance in O’Higgins and Biobío regions with values above 0.8 and in Araucanía and Los Ríos regions with values below 0.77. The adjusted coefficient of determination $\left(R_j^2\right)$ (Figure 10B) also shows skilled fitting (over 0.8) from O’Higgins to Biobío regions, although the Maule basin shows the lowest estimate for $R_j^2$. The same occurs in Araucanía and Los Ríos regions. The $F_{\text{value}}$ from the ANOVA shows an excellent fit of the model with a $P_{\text{value}}$ less than 0.05 (Figure 10C). In turn, the average model performance for ${\rho }_{xy}$ is higher than 0.9 in all domains (Figure 10H). We also find an average $pbias$ below $22\text{\%}$, with O’Higgins and Biobío below $15\text{\%}$ (Figure 10F). The outliers identified for ${\sigma }_e$ (Figure 10G) exhibit an error below 0.22, supporting that the model is skilled. $nrmse$ (Figure 10E) shows average model values below 0.4, while values above 0.48 are shown in Araucanía and Los Ríos. In these regions, a high error confirms the model is more uncertain in estimating extreme precipitation events. Table 3 shows the results for the calibrated model configuration ($R^2$, $R_j^2$, $F_{\text{value}}$, and $P_{\text{value}}$) and error estimate (${\sigma }_e$, $nrmse$, $pbias$, and ${\rho }_{xy}$) for each regional domain.

FIGURE 10. Box diagrams corresponding to the adjustment of the topo-climatic model on the domains in execution. (A) Multiple determination coefficient $\left(R^2\right)$, (B) adjusted multiple determination coefficient $R_j^2$, (C) value F, (D) value P, (E) standardized mean square error $\left(nrmse\right)$, (F) percentage of systematic error $\left(pbias\right)$, (G) standardized standard deviation $\left({\sigma }_e\right)$, (H) Pearson correlation coefficient $\left({\rho }_{xy}\right)$.

4.3.2. Zonal and Basin Analysis

The northern zone (Figures 11A-C) shows an average ${\rho }_{xy}$ of 0.88, although some stations in Rapel get values close to 0.6. Meanwhile, the $nrmse$ exhibits an average of 0.46, although in Rapel and Maule some cases are above 0.55. Interestingly, ${\sigma }_e$ (Figure 12A) shows an average fit of 0.07, with that in Maule below $-0.4$. $pbias$ (Figure 12B) shows an average adjustment of $18\text{\%}$, but Rapel and Maule show the most considerable difference with an estimate over $50\text{\%}$.

FIGURE 11. Taylor diagrams corresponding to the performance of the topo-climatic model on the river basins (Figure 1B): (A) Rapel, (B) Mataquito, (C) Maule, (D) Itata, (E) Biobío, (F) Imperial, (G) Toltén, and (H) Valdivia.

FIGURE 12. Box diagrams corresponding to the adjustment of the topo-climatic model on the domains in execution. (A) Standardized standard deviation $\left({\sigma }_e\right)$, (B) percentage of the systematic error $\left(pbias\right)$, (C) Nash–Sutcliffe coefficient $\left(nse\right)$, and (D) Willmott coefficient $\left(d\right)$.

The central zone (Figures 11D,E), ${\rho }_{xy}$ shows an average of 0.89 and Biobío exhibits values below 0.8. The $nrmse$ displays an average of 0.52, while Biobío shows an irregular value greater than 1. The average for ${\sigma }_e$ is 0.1 (Figure 12A) with only Biobío being below $-0.5$. For $pbias$ (Figure 12B) there is a $24\text{\%}$ average error with Biobío showing an estimate below $100\text{\%}$.

The southern zone (Figures 11F-H) exhibits an average ${\rho }_{xy}$ over 0.9, while the lowest estimates are above 0.85. $nrsme$ shows an average of 0.52, with all basins having some values above 0.66. On the other hand, while ${\sigma }_e$ (Figure 12A) shows an average of 0.1, only the Imperial basin exhibits a low estimate of $-0.2$. $pbias$ (Figure 12B) shows an average of $24\text{\%}$. $nse$ shows an average above 0.7 in all basins (Figure 12C); only two atypical cases are detected in Maule and Biobío with values of 0.5 and $-1$, respectively. On the other hand, d (Figure 12D) shows an average fitting of 0.75. Only Biobío exhibits an outlier in the Liucura station with a value below 0.5.

On each regional domain, a different monthly model is adjusted to the calibration data and checked using validation points. This allows us to construct a monthly solution for each regional domain, delivering a dataset of monthly fields covering the period $2000\text{–}2011$ that subsequently are used to produce monthly climatologies from January to December. Figure 13 shows the observed accumulated annual rainfall and the one obtained by the topo-climatic dynamic model over the study area $\left({34}^{\circ }\text{S}–{41}^{\circ }\text{S}\right)$. Figure 13B indicates an adequate distribution of precipitation from north to south and from the ocean to the Andes. There is a progressive increase from north to south and drastic increases in mountain areas, particularly along the Andes and the coastal cordillera. Thus, in sectors close to the coast (from ${38}^{\circ }\text{S}$), an abundant annual precipitation is identified, sometimes reaching over 2,000 mm. Throughout the Andean areas, the accumulated annual rainfall can reach and exceed 3,000 mm per year. Finally, between the coast and the Andes below ${36}^{\circ }\text{S}$, it is possible to identify a decrease in annual precipitation.

FIGURE 13. Annual cumulative rainfall. (A) Rainfall observed and (B) dynamic topo-climatic model.

5. Discussion and Conclusion

High-resolution precipitation spatial variability is significant in different fields such as hydrology, environment, forestry, and agriculture. In Chile, several authors worked on gridded precipitation datasets built only with the information available from local stations in some areas of the country (Román and Andrés, 2010; Jacquin and Soto-Sandoval, 2013; Reyes, 2013; Castro et al., 2014; Sijinaldo, 2015). Also, in some cases, datasets were generated from observed data and low-resolution gridded data higher than ${0.25}^{\circ }$ (Morales-Salinas et al., 2012) or even poorer quality datasets. Thus, these high-resolution datasets are inadequate to study large areas and/or the whole country. Moreover, the insufficient spatial distribution of the pluviometric network (Muñoz et al., 2018) precludes real hydro-climatic impact under current environmental conditions.

Therefore, we assert that the dataset will be useful for different communities in Chile. It will be useful for agricultural, forestry, and livestock planning, for example, and for national institutions such as the Dirección General de Aguas (DGA), the Instituto de Investigaciones Agropecuarias (INIA), and the Instituto Forestal de Chile (INFOR), among others. This product can serve as a basis of comparison for studies aimed at investigating future climate or agricultural changes (Orrego et al., 2016), allowing investigation of the role of and Antartic Oscillation and many climatic phenomena on some local precipitation-related aspects (Fustos et al., 2020b), which is beyond the scope of this study.

To get an improved dataset and to quantify precipitation, we built a high-resolution dynamic statistical gridded precipitation product at $\sim 800\text{\hspace{0.17em}m}$. Its development involved time series of local precipitation and different global precipitation grid sets (as shown in Table 2). To detect spatial-temporal variability of precipitation, we performed numerical modeling of the atmosphere over the study area at 4 km spatial resolution for the period $2000\text{–}2011$. These outputs are a powerful tool for statistical scale reduction (Widmann et al., 2003; Schmidli et al., 2006; González-Rojí et al., 2019). Also, within each region, according to Eqs. 1and2 requirements, we performed a month-by-month data selection. We also evaluated precipitation time series belonging to eight basins to estimate the responses of the best products that entered our dynamic topo-climatic model (Figure 11).

We evaluated some precipitation grids and the results are presented in Table 2. To classify the optimal precipitation estimates and reduce uncertainty, we used a selection criterion (Eq. 2). As shown by Zambrano F. et al. (2017) and Zambrano M. et al. (2017), the performance of many is unsatisfactory. The products $\text{CHIRPS}$, $\text{MSWEP}$, and $\text{PGFv}3$ show the most accurate estimations $\left({\rho }_{xy}>0.8\right)$. Meanwhile, $\text{CMORPH}$ and $\text{SM}2\text{RAIN}$ provide the worst performance over the whole study area $\left(<0.6\right)$, and the average of ${\sigma }_e$ displays a similar distribution in all basins. However, the $\text{CMORPH}$ set exhibits the highest mismatches, with estimates close to 1 in Biobío and Imperial. $pbias$ shows a similar distribution along all basins between $10$ and $30\text{\%}$. The Nash–Sutcliffe coefficient ($nse$) shows a decreasing yield from Rapel to Itata $\left(>0.7\right)$ along with an increasing dispersion. For example, at Biobío and Imperial, the adjustment is centered on 0.7 with low dispersion. Also, Toltén and Valdivia have an average of 0.7, but with high dispersion. On the other hand, Willmott (d) does not show an improvement in the performance of the datasets, although it shows a similar pattern as $nse$$\left(>0.7\right)$. For all basins, we notice that $\text{SM}2\text{RAIN}$ shows low values. However, model fitting results for some datasets should be seen with caution, since we eliminated any value above or below $5\text{\%}$ of the respective maximal or minimal observed precipitation. Therefore, cleaned data improve the estimates of these sets.

Atmospheric dynamic models, like WRF, simulate precipitation based on convective and large-scale processes (such as precipitation from cumulonimbus fronts or clouds), generating an improved distribution of precipitation as well as other fields such as temperature (Pope and Stratton, 2002; Jung et al., 2006; Gent et al., 2010). However, there are still errors attributed to feedback processes, which are not necessarily reduced by increasing the model spatial resolution. As these errors are pronounced on rough terrain, and even more so in mountainous areas, to take that possibility into account, it was necessary to add a $5\text{\%}$ correction. Despite this, the limited density of observation along the Andes incorporates a high degree of uncertainty, preventing homogeneous fit, thus making the construction of a correct model at high altitude difficult. Taylor and boxplot diagrams show that the model presents a remarkable performance in the northern zone (Figures 6, 7). As far as the central zone is concerned, adequate fitting is again obtained for ${\rho }_{xy}$, $nse$, and d.

To construct our dynamic geostatistical model, we used global precipitation grids and a local atmospheric model (Table 2) in addition to regional domains (Figure 1B). This is a nontrivial process, particularly in regions with a great variety of space–time patterns (Zhang et al., 2016). During the precipitation reconstruction phase, we verified the elevation coefficients relative to each regional domain with available local information. For the mountains, the little local information available prevents adequate correction and adjustment of the WRF model or the gridded data. Consequently, in the heights, the precipitations should generally be overestimated. Each month, the adjustment provides different elevation coefficients for each regional domain. However, for the entire domain, the solution that imposed a single elevation coefficient, considering all the information from north to south, turned out to be the best.

The results of the model chosen for each region (Table 3) show reasonable skill in the assessed basins (Figures 11, 12). However, the model is unable to satisfactorily represent extreme precipitation events, which gradually increase from north to south. It is essential to note that the choice of the regression model is conditioned by the input data selection and the domain used. Only numerical modeling and satellite information are the ones that provide adequate spatial coverage and persistence. In turn, statistical downscaling can improve the accuracy of discrete meteorological measurements (Diez et al., 2005; Fernández-Ferrero et al., 2009; Le Roux et al., 2018; González-Rojí et al., 2019). To our knowledge, this kind of study is scarce in Chile, not only in the construction of grilled climatic data but also in the analysis and comparison of different methodologies for an optimal statistical scale reduction model or the domain type needed to carry this out.

TABLE 3. Statistical summary of the type of topo-climatic model chosen by the political domain of each region. (wslp), (wcost), (all), (wcon).

The database presented in this study results from atmospheric modeling, corrected, adjusted, and ported to higher resolution by a dynamic topo-climatic method using data from local meteorological stations and satellite fields. Note that to improve the accuracy of database selection by different statistical methodologies, we visually inspected each dataset in detail and month by month. However, we still depend on the quality of the atmospheric modeling, which is not reality, and that of the satellite fields and local data. Therefore, we would like to point out that although we are confident that the method used and the resolution will allow a better understanding of local hydrological variability, this database does not replace reality. That is why, contrary to several global gridded databases, which do not stipulate it, we recommend within each basin to check against local data and correct any existing bias.

Of course, there are still many improvements to be made. Note that we will always be dependent on the distribution of the local weather network, which is the only one that allows the calibration and validation process. It is fundamental to improve the quality in many places, notably in the slopes of the Coastal and Andean mountain ranges, where weather coverage is almost nonexistent. A nonhomogeneous distribution of altitude data, used to build a solution over the mountains, may not represent the different mechanisms in both cordilleras, sometimes separated by less than 100 km. Thus, specific measurement campaigns are needed to support, for example, model calibration during brief periods and on a local scale. However, many phenomena throughout Chile develop at high frequencies, where the monthly value does not allow for a correct description (Fustos et al., 2020a). Thus, while it is necessary to expand time–space coverage, north–south and up to 2020, in addition to build other high-resolution grids of essential hydro-climatic parameters (such as temperature), it is also imperative to focus on daily grid products.

Data Availability Statement

The datasets generated for this study are available on request to the corresponding author.

Author Contributions

RA assisted the first author (FA) at different levels of research, such as conception, design, analysis, and editing of the manuscript. Although RA is responsible for the overall concept and design of the research, data collection, development of atmospheric modeling, and statistical analysis of this manuscript was performed by FA; RA and AA contributed to the technical implementation required for atmospheric modeling; and AA improved the mathematical interpretation of the results and participated in the review of the submitted manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We appreciate and are indebted to both the reviewers’ outstanding work, whose rigor and scientific culture helped us improve both form and substance. We are very grateful and indebted to our colleague A. Fernandez, for careful proofreading of the manuscript. We are also grateful for advices by B. Mark. We are grateful for funding support received from the PhD in Physical Sciences program from the Faculty of Physical and Mathematical Sciences, the Research Directorate (DRINV), and Postgraduate Directorate of the University of Concepción. This research was partially supported by the supercomputing infrastructure of the NLHPC (ECM-02) at Centro de Modelación y Computacion Científica, Universidad de La Frontera, CMCC-UFRO.

References