Comparison of spatial interpolation methods for annual and seasonal rainfall in two hotspots of biodiversity in South America

: The Colombian Biogeographic Choco (CBC) and the La Plata Basin (LPB) are regions with high biodiversity. However, these areas are characterized by scarce climatological information, complex orography, and rain-gauge network unevenly distributed. Interpolated data from the ground station might overcome these aspects. For this reason, is necessary to identify the best technique for the spatial interpolation of rainfall. Hence, the spatial interpolation techniques were applied to annual and seasonal rainfall in the CBC and LPB. Geostatistical results and deterministic approaches were compared by cross-validation. Cokriging with spherical (gaussian) model is the best interpolator in the CBC (LPB), as indicated by the lowest root mean square error (RMSE) and a standardized RMSE close to one. The CBC shows three rainfall cores: the northern, 9,000 mm/year; the central-southern, 10,000 mm/year; and the southern, 7,000 mm/ year. The LPB shows a west-east rainfall gradient, with a minimum to the west (450 mm/ year) and a maximum in the mid-west (2,000 mm/year). To the north of the LPB, rainfall reaches 1,500 mm/year, while in the south it reaches only 900 mm/year. The results in our study may be useful for scientists and decision-makers for use in environmental and hydrological models for the CBC and the LPB.


INTRODUCTION
Climate is one of the most important environmental factors in terrestrial ecosystems (Tuhkanen 1980, Chiu et al. 2009). Rainfall plays a key role in the hydrological cycle, a determinant of the global climate system and participates in the dynamics and atmospheric composition (Sánchez & Vélez 2015). The high complexity of the spatial and temporal structure of this climatic element hinders the generalization of typical behaviors, even in small regions. However, climate information is generally recorded as timely data through the weather stations at each location, and many environmental and hydrological models require input information at unobserved locations, requiring spatially continuous climatic data, usually in the form of interpolated grids; furthermore, a common difficulty in the study of events associated with rainfall is the lack of timely and reliable information (Guisan & Zimmermann 2000, Anderson & Martínez-Mayer 2004, Silva 2009). In this respect, spatial interpolation techniques are a trustworthy approach to estimating climatic information from close measurements for locations without observations (Berndt & Haberlandt 2018).
To overcome the low-density network constraint and sparse distribution of the meteorological stations, the spatial interpolation becomes an important tool to estimate climatological variables not geographically covered by the existing observation network (Andrade & Moreano 2013, Basconcillo et al. 2017, Berndt & Haberlandt 2018. However, there is little evidence that a single interpolation method is ideal for several conditions. Thus, it is important to determine the best method for each situation (Atorre et al. 2007). Essentially, this happens because each technique depends on the dataset characteristics. So, a technique may be suitable for some variables or regions but may not work for others (Basconcillo et al. 2017). The generation of continuous surfaces can be performed by a variety of methods, but the difficulty is to choose the one that best reproduces the real surface (Caruso & Quarta 1998, Ly et al. 2011. Concerning the rainfall, it may be influenced by small-scale processes and by orography. Given the orographic diversity and the influence of many atmospheric factors in South America (SA), the local and regional climates show high complexity (Poveda 2004) with diverse climatic dynamics such as the Colombian Biogeographic Choco (CBC) and the La Plata Basin (LPB), regions renowned worldwide as biodiversity hotspots for conservation prioritization (Myers et al. 2000, Mittermeier et al. 2004. Most biodiversity hotspots are in tropical developing countries that face great challenges, such as high demographic pressure, food shortage, poverty, and corruption (Veech 2003, Williams 2011. The CBC and LPB are in extreme northwestern and southeastern SA, respectively. They are regions with climate and biodiversity influenced by the moisture recycling, in particular, evapotranspiration in forests of the Brazilian Amazon the continent's biggest rainmaker which often takes water long distances that contributes substantially to rainfall regionally as well as over other remote regions such as the CBC and LPB (Zemp et al. 2014).
In the CBC and LPB, the spatial distribution of rainfall is an aspect important that requires detailed scientific research since are regions with rainfall measurements scarce and with raingauge network sparse and irregular. Hence, it is necessary the spatial interpolation process where points with known values are used to estimate unknown values at other points. This type of interpolated surface is often called a statistical surface, that describes and explains the spatial trend of the variables (Mitas & Mitasova 2005). There is a wide range of rainfall interpolation methods proposed in the literature, divided into three main categories: deterministic (geometric), statistical and geostatistical. Deterministic interpolation methods create surfaces from the measured points, based on similarity, such as polygons of Thiessen, also known as nearest neighbor or Voronoi diagrams (Thiessen 1911) and Inverse Distance Weighting (IDW) (Shepard 1968) or degree of smoothing, such as Radial Basis Functions (RBF), a generalized version of the multi-quadratic method (MQ) developed by Hardy (1971); besides the irregular triangular network (TIN) or linear regression and neural networks (Sluiter 2008). Statistical methods use linear or multiple regressions to correlate rainfall with predictive variables, such as elevation, longitude, distance to the sea, among others (Castro et al. 2014). However, the greater the number of predictors considered in the regression, the greater the number of stations required in the estimation (Chen & Liu 2012). The kriging method (Matheron 1971), creates surfaces that incorporate the statistical properties of the measured data, considering spatial autocorrelation between known points (Firdaus & Talib 2016, Angulo et al. 2009, Hao & Chang 2013. These techniques produce not only predictive surfaces, but also surfaces of error or uncertainty, thus give an indication of how good the predictions are (Firdaus & Talib 2016).
Univariate geostatistical methods (simple kriging) generally smooth the interpolated variable, and therefore have difficulty in reproducing the spatial variability accurately. To improve the interpolation performance (Wagner et al. 2012), multivariate methods use complementary spatial information from covariates such as elevation or dynamic radar data (Foehn et al. 2018), since the inclusion of geographic and topographical information increases the estimation capacity of the interpolation schemes (Portalés et al. 2010, Borges et al. 2016. Studies have incorporated radar information into interpolation methods (Schiemann et al. 2011, Verworn & Haberlandt 2011, Wagner et al. 2012). However, these methods require a large amount of data that is often not available to users (Borges et al. 2016). Goovaerts (2000) using the elevation as secondary data incorporated into the multivariate geostatistics for monthly and annual rainfall, compared these results with the deterministic methods, and found that the IDW and the Thiessen Polygon provide greater prediction errors. More recently, Pellicone et al. (2018) compared different rainfall interpolation algorithms in southern Italy to identify the method that best reproduces the surface of the rainfall field; the results clearly indicate that geostatistical methods outweigh the IDW, and that kriging with external drift (KED) shows the smallest prediction error.
In the case of Colombia, Álzate et al. (2018) used the regionalized rainfall model Regionalisierte Niederschlage (Regnie) to interpolate rainfall in the Andean, Caribbean and Pacific regions, integrating slope and elevation as secondary variables. Multiple linear regression models and geoprocessing tools were used to generate the interpolated surfaces. The Regnie surface tests were like those obtained with Spline and IDW interpolations for rainfall. The Institute of Hydrology, Meteorology and Environmental Studies (IDEAM 2005) has applied the IDW method to represent the main climatic variables for the Colombian territory. On the other hand, Correa et al. (2014) found that circular, spherical and gaussian models with the ordinary kriging interpolation method can be used with satisfactory performance in interpolating the annual rainfall data in the state of Mato Grosso in southern Brazil (high basin of LPB).
In this context, the present study compares different spatial interpolation techniques commonly used, aiming to identify the technique with the best rainfall estimation in the CBC and LPB. Furthermore, the performance of the spatial interpolation techniques is explored through its capability to represent the annual and seasonal rainfall. The techniques evaluated here include non-geostatistical (IDW) and geostatistical approaches (Ordinary Kriging and Ordinary Cokriging with altitude as a secondary variable) with different theoretical models (exponential, spherical, and gaussian). The segment below presents the study domain and data, followed by the description of the methods used. The results and discussions are presented in Colombian biogeographic Choco and in La Plata Basin.

Study domain and data
Analyses are carried out in two regions in SA: the Colombian Biogeographic Choco (CBC) and the La Plata Basin (LPB). The Colombian Pacific basin located in extreme northwestern SA is composed of a coastal tropical forest and a variety of ecosystems proper among mangroves, marshes, flood forests and moorland that structures an enclave of the specific diversity known as the CBC (Figure 1) (Myers 1988, 1990, Myers et al. 2000, Mittermeier et al. 2004, Marchese 2015. This region hosts approximately 3% of the global plant species, and is characterized by high rainfall, up to 8,000 mm/year (Poveda & Mesa 2000). Despite being a region of great importance for Colombia, it is an area with scarce climatological information due to the orographic complexity, the difficult access and the little instrumentation and investment of the Colombian state in the region. On the other hand, LPB in southeastern SA is the fifth largest basin in the world and the second most extensive in SA after the Amazon ( Figure 1) (García & Vargas 1998, Mechoso et al. 2001. It covers ecosystems of special interest such as the Cerrado and the Atlantic Forest; however, does not realize its full potential due to the lack of coordinated efforts of the five countries sharing the basin (Saurral & Barros 2009). These two regions are part of the biodiversity hotspots, which cover only 17.3% of the terrestrial surface but retain 77% of all endemic plant species, 43% of vertebrates and 80% of all endangered amphibians (Mittermeier et al. 2011).
Data used in the study include monthly rainfall averages of 193 meteorological stations provided by the meteorological and hydrology institutes of LPB countries and 166 for the CBC ( Table I). The stations were selected according to criteria such as: extent of the historical records reported by each institute, location of stations within the basin and in the nearby basins, and those with less than 10% of missing data. Based on these criteria, of the 193 LPB stations, 170 were selected with historical series of rainfall spanning over 31 years with a common period between 1987-2017. The stations in the LPB exhibit an irregular spatial distribution, with some information deficits in the north, mainly in the State of Mato Grosso (Brazil), and in the southeast of the basin ( Figure 1). On the other hand, the 156 stations in the CBC present data during a common period from 1983 to 2016 and show a more homogeneous distribution; however, the northern coast offers a lower density of information ( Figure 1).
To obtain coherent historical series, the data were evaluated for absence of errors, completeness and consistency. The data consistency was tested with missing data analysis (MDA) by applying the linear regression with data from the nearest stations. MDA, also called data imputation, is a procedure that uses the information contained in the sample to assign a value to those variables that have records with missing values, either because there is no information or because it is detected some unexpected behavior. Then, the exploratory analysis was performed, and the frequency charts, multiannual average monthly rainfall and quantil-quantil graphs were obtained.
Elevation data were obtained from the Shuttle Radar Topography Mission (SRTM; National Aeronautics & Space Administration-NASA 2017). Data was extracted for each station.

Interpolation methods
To find the best interpolation technique for each region and differences in seasonal performance, we compared the IDW, Kriging and Cokriging techniques with different semivariogram models (exponential, spherical and gaussian), and altitude as a secondary variable, for a total of seven interpolation tests.

Inverse distance weighting (IDW)
The IDW, the method most widely used by GIS analysts (Firdaus & Talib 2016), employs the first Law of Geography (Tobler 1970), and estimates unknown measures as weighted averages of measures known in the near points, which has the greatest weight in the procedure (Longley et al. 2011). That is, each climatic value at a nonsampled point z(x) is a scoring average of the distance from the values at the sampling points in the neighborhood z(x 1 ), z(x 2 ), ..., z(x n ). Climatic values are more like nearest distances, so the inverse distance (1/d i ) between z(x i ) and z(x) is used as a weighting factor (Eq. 1).
where z(x) is the predicted value, z(x i ) is the climatic value at a neighboring meteorological station, d ij is the distance between z(x) and z(x i ) and r is an empirical parameter. The ideal weighting value can be obtained by minimizing the root mean square error, which is a statistic calculated during cross-validation (Hao & Chang 2013).

Kriging and cokriging
The Kriging method considers both the degree and distance of variation between observed points. If the spatial variation of an attribute is not totally random (stochastic) or deterministic (Ly et al. 2011), and the spatial variation of a continuous climate variable is irregular to be modeled by a continuous mathematical function, the spatial variation can best be predicted by a probabilistic surface (Angulo et al. 2009). This continuous variable is called a regionalized variable, which consists of a derivation component and a spatially correlated random component (Burrough et al. 2015). Thus, the spatially localized climatic variable z(x) is expressed by Eq. 2.
the structural variation of the climatic variable, ɛ'(x) are the spatially correlated residuals, or the difference between the derivation component and the sampling data values, and ɛ'' is the spatially independent residue.
A function that relates the spatial variance of the variable is determined using the semivariogram that indicates the semi variance ( γ ) that describes how similar the observed values are grouped in space (Tobler 1970). The calculated variance is a measure that determines the similarity between observations, where the highest similarity, the lower semivariance (Lozano et al. 2004). Kriging focuses on the correlated spatial component and uses the adjusted semivariogram (Atorre et al. 2007). On the other hand, the Cokriging allows to consider the influence of external variables (co-varied, in this case height) when analyzing the crosscorrelation between the errors of the different variables ɛ' 1 (x), ɛ' 2 (x), etc. (Angulo et al. 2009).
Since the geostatistical estimates are based on the spatial structure and the spatial variability of the data, the semivariogram is the most appropriate means to represent the spatial dependence of the data (Aziz et al. 2019). This spatial statistic tool allows to express the dissimilarity as a function of distance (Sabzipour et al. 2019). Essentially, a semivariogram is a description of the spatial data continuity (Pellicone 2018). The experimental semivariogram is computed as half the average squared difference between the components of data pairs (Eq. 3) to quantify the spatial dependence on the data.  (2000) observed that the use of multiple secondary variables could lead to unstable cokriging systems. Therefore, only the elevation in this study was considered to improve the spatial interpolation. Thus, cokriging requires the semivariogram models of each variable γ (Z i , Z j ) and γ (Y i , Y j ), equal to precipitation and elevation, and cross semi-variograms of primary and secondary variables respectively for spaced distances γ (Z i , Y i ). The cross semivariogram was computed as Eq. 4.
However, the difficulty lies in the fact that the three models cannot be constructed independently from one another. Thus, the easiest approach to estimating rainfall consists of modeling the three semi-variograms as linear combinations of the same set of basic semivariogram models (Ly et al. 2011, Goovaerts 1997. The coefficients of the adjusted models are used to determine the weight of the stations in the interpolation. Besides, the estimated variance depends on the semivariogram model, the number N of rain gauges, and its spatial location. In this way, an optimal choice of the semivariogram model is crucial for the evaluation of the data (Bohling 2007). Nevertheless, there are several different semivariogram models, and it is a difficult task to determine among the models which one produces the closest results to the observed conditions (Mazzella & Mazzella 2013, Nadiah et al. 2016, Aziz et al. 2019. The method's advantages and its disadvantages hence depend strongly on the characteristics of the dataset used to define their suitability. Also, Castro et al. (2010) point out that many studies do not define the theoretical model that best fits the data studied. Based on the above, the study aimed to compare different semivariogram models to select the one that best suits to obtain the spatial distribution of annual and seasonal rainfall in the regions under study. Here, the three more used models for fitting semivariogram, the spherical, gaussian and exponential models (Mcbratney & Webster 1986, Goovaters 1997, Deutsch & Journel 1998, Ly et al. 2011, Mazzella & Mazzella 2013, Aziz et al. 2019, were used. Therefore, the performances of the spherical, gaussian, and exponential semivariogram models were evaluated in this study. The equations of the semivariogram models can be seen in Ly et al. (2011). Then, the best semivariogram model was selected based on cross-validation statistics, which are described in the next section.

Evaluation of error assessment
The IDW, ordinary Kriging and Cokriging interpolation methods were explored and compared. The performance of each method was evaluated through cross-validation. This procedure compares interpolation methods by repeating the following procedure for each interpolation (Chang 2006): (1) eliminates a known point from the dataset, (2) the remaining points are used to estimate the value at the previously removed point and 3) the error of the estimate is calculated by comparing the estimated value with the known one. After completing the procedure for each known point, two common diagnostic statistics, the Root Mean Square Error (RMSE, Eq. 5) and the standardized RMSE (RMSSE, Eq. 6), are obtained to evaluate the accuracy of the interpolation method as is shown in the following equations: where Zi and Z are the values measured and estimated at the sampling point i (i=1, 2, ... n); n is the number of values used for the estimate; and S is the standard error.
The RMSE statistic is available for all exact local methods, but RMSSE is only available for the Kriging and Cokriging because variance is required for the computation. The best interpolation method should produce a smaller RMSE and a RMSSE closer to 1 (Chang 2006). Additionally, the Average Standard Error (ASE, Eq. 7) was calculated to assess the variability of predictions: ( ) If the ASE is close to the RMSE, the variability in the prediction is evaluated correctly. In addition, if the ASE is greater (less) than the

RESULTS AND DISCUSSIONS
Colombian Biogeographic Choco to the RMSE value, exponential and spherical cokriging show the best prediction values of the variability in the seasonal scale. However at the annual scale, only the spherical model has good predictive capability. This shows that spherical semivariogram for cokriging is the best fitted experimental semivariogram for the study of annual and seasonal CBC rainfall. Besides, since the ASE is greater (less) than the RMSE, the cokriging variance is larger (smaller) than the true variance and indicates that the variogram model overestimates (underestimates) the prediction variability mainly in March-May (MAM), July-August (JJA) and September-November (SON) (DJF and Annual), specifically for Cokriging spherical. These results are confirmed in the Quantil-Quantil graphs (Figure 2), which indicate that the values predicted by the spherical cokriging represent with more precision the distribution of the observed data. This occurs because the CBC rainfall is strongly associated with the topography of the Western Cordillera region, the Atrato and San Juan valleys, the Baudó Mountain Range and the coastal plain (Guarín & Poveda 2013). The surface winds of the Pacific Ocean interact with the eastern trade winds on the western cordillera, which combined with the surface warming effect and the orographic rise favor the deep convection, elevation of moist air, high amounts of condensation, and therefore, high rainfall (Poveda & Mesa 1999, Trenberth 1999. Thus, the prediction could be improved considering elevation as a correlated secondary variable. Figure 3 presents the results of the interpolation performed with cokriging for the spherical model. Also, the interpolation results show three main cores of strong rainfall, which occur on the western flank of the western ridge, over the middle part of the region in the westeast direction. The first one is located in the central-north of the CBC, in the sub-region of the Northern Pacific, where the smallest rainfall occur in the first semester of the year and the largest between July and November (variation between 700 mm and 800 mm). The second core located in the Patía Basin shows the highest amounts of rainfall in the region, with seasonal values over 800-900 mm during MAM and SON. The last nucleus of higher rainfall is in the South Pacific, where there is greater rainfall during the first semester of the year (600 to 700 mm seasonal), and the lowest in the JJA, with rainfall between 400 and 500 mm.
On the other hand, the eastern flank of the western mountain range exhibits the lowest values of precipitation, mainly in the southern and northern extremes in the Pastos Knot and the Middle Cauca River sub-regions respectively, however, in different seasons of the year. During the boreal summer-JJA (winter-DJF), the smallest rainfall occur in the southern (northern and central) part of the region, when the Intertropical Convergence Zone (ITCZ) is in its northern (southern) position (Figure 3). The analysis of the directional trends of annual precipitations indicates a marked difference between the western flank (CBC) and east (geographical valley of the Cauca river) of the western mountain range. On the western flank, the highest rainfall is concentrated at the mountain slopes and fell to the northern and southern, with lower values in the northern and higher parts of the mountains. In the case of the eastern flank, the rainfall increases from west to east and from south to north (Figure 3).
The results presented here suggest that the CBC rains are unevenly distributed, with three intense precipitation cores throughout the region, reaching average values between 3,000 and 11,000 mm annually (Figure 3). The lower intensity occurs in the subregion of Urabá (less than 3,000 mm/year) and in the South Pacific subregion (precipitation from 3,000 to 7,000 mm/year). High annual rainfall occurs in the subregion of the Northern Pacific, with mean values between 8,000 and 10,000 mm/year, and in the Patía Basin, between the states of Valle del Cauca and Cauca with annual totals between 7,000 and 10,000 mm. This rainfall distribution agrees with previous findings for the Colombian Pacific (Eslava 1994, Poveda 2004, Guzmán et al. 2014. These results contrast with the rainfall regime on the eastern side of the western cordillera, where rainfall in the High Cauca River and Pastos Knot reaches mean values between 1,000 and 2,000 mm/year. These results are consistent with Poveda et al. (2004), who established ecogeographic subregions for the CBC with three major climatic belts: the South Pacific subregion, the Northern Pacific subregion (mainly in the upper parts of the Baudo and El Carmen de Atrato basins), and the Patía Basin. The first subregion, in the southern part, presents low humidity characteristics and annual rainfall between 730 and 3,318 mm/ year, the second in the central northern part, wetter conditions and rainfall between 8,494 and 13,670 mm/year, and the last subregion in the southern-central part, rainfall between 5,909 and 8,494 mm. On the other hand, in the subregional scale Hurtado (2009) identified a strong variability of the precipitation in small distances, which is determined by gradients induced by the topography, mainly in the Pacific plain and the mountain slopes, results that correspond to those found in this study. Figure 4 shows the relationship between annual rainfall and altitude, indicating a large difference between the stations of the western margin (Pacific slope) and those of the eastern (Andean slope). In general, the western slope stations have an average altitude of 250 m.a.s.l., and average annual rainfall higher than 5,000 mm (e.g., Bellavista, Bocas de Patía, Tutunendo, Junín and Pto. López); some stations such as Carmen de Atrato and La Cumbre, located between the limit of the two slopes (1,500 m.a.s.l.) present values below 3,000 mm/year. While the eastern slopes have an average height of 1,400 m.a.s.l and rainfall less than 1,700 mm/year (e.g., Tormento, Clarita, Toscana and, Univalle). Also, some stations in the western margin of the central mountain range of the Colombian Andes -Nariño, Potreros and, Paraiso-present a mean rainfall between 2,000 and 3,000 mm/ year. Conversely, less than 1,000 mm/year has been recorded for stations in the eastern margin of the western mountain range such as Vijes and Guasca.
The difference between the two slopes is mainly due to the uniform rainfall regime throughout the year in the western region, without a defined dry season, due to the persistence of the ITCZ, the convergence of trade winds in the equator and the Choco Low-Level Jet (CJ), which interacts with the mesoscale convective systems (Poveda & Mesa 2000, Zea et al. 2000. Meanwhile, the Andean region experiences a bimodal cycle with two wet stations, one in MAM and the other in SON, and two dry stations during JJA and DJF, due to the migration of the ITCZ. This influence is combined with the orographic effects of the Western Cordillera of the Colombian Andes, which acts as a barrier to the moisture transport from the Pacific Ocean that discharges more moisture on the western slope through orographic rainfall (Poveda et al. 2002, 2014, Trojer 2018. Besides, the differences in the Andean Region have been previously documented by Sedano (2017), which observed greater rainfall on the central mountain range, associated with a moisture transport from the Pacific Ocean that crosses the western mountain range and enhances moisture convergence on the western margin of the central mountain range between 1,600 and 2,000 m. a.s.l. which tend to overestimate the rainfall in SON (ASE > RMSE) and underestimate in the other stations. However, on the annual scale, the results of gaussian cokriging indicate a better ability to predict the rainfall variability. The results suggest that the gaussian semivariogram for cokriging is the best fitted experimental semivariogram for the study of LPB rainfall. Following the previous results, Figure 5 shows that the values predicted by gaussian cokriging are better representing the distribution of observed data. So, the annual and seasonal rainfalls of the LPB noted in Figure 6 represent the results of cokriging with gaussian model. The annual rainfall for the 1987-2017 period shows a west-east gradient with a minimum of 450 mm in the highest part of the LPB in Bolivia and highest values toward the center-west, with maximum values of 2,000 mm on the triple border between southern Brazil, Paraguay and northeastern Argentina, near the Atlantic Ocean ( Figure 6). On the other hand, rainfall is also abundant in the north of the basin, where values reach 1,500-1,650 mm in the Cerrado and Pantanal regions, while in the humid Argentine Pampas rainfall reaches only 900 mm per year. This distribution is consistent with previous results (Boulanger et al. 2005, Barros et al. 2006, Bidegain et al. 2017. However, rainfall maps provide more information about the Chaco and Andean region. The annual cycle clearly shows the lowest precipitation in the high Andes, especially in the west and southwest; highlighting the relationship between the annual precipitation and the altitude of 170 stations located in the LPB. Only a limited number of stations are more than 1,500 meters above sea level (m.a.s.l.) and all are less than 1,000 mm/year (Padilla, Azurduy, Sucre, Higueras and Comarapa stations; Figure 7), those located at more than 3,000 m.a.s.l. have values lower than 500 mm/year (Potosí, Pabellón and Villazón stations; Figure 7). Similar results were found for the stations of the Andean region in the Amazon Basin (Espinoza et al. 2009) and in Bolivia (Ronchail & Gallaire 2006). At low altitudes, stations that register more than 2,000 mm/year are lower than 1,000 m.a.s.l. and close to the Atlantic Ocean such as Paranaguá and Bernardo in Brazil, and Iguaçu and Campo Grande in Argentina (Figure 7).

La Plata Basin
The seasonal rainfall patterns in Figure 6 show considerable differences between seasons. During DJF, the highest values are over the north (200mm) and the lowest over the south of the basin (< 120mm). During JJA, the highest precipitations are concentrated on the eastern flank of the basin (60 to 140 mm), close to the Atlantic Ocean, while the north and western regions present values less than 60 mm quarterly.
During MAM and SON, the climatological patterns show a better distribution of the rains on the LPB, however, a core of precipitations is perceptible on the Center of the LPB, that diminishes towards the west of the basin; this core of precipitations presents higher values during SON (140 to 200mm) in comparison with MAM (120 to 180mm). The climatological results coincide with the seasonal precipitation regimes found by several authors (Berbery & Barros 2002, Vera et al. 2002, Gan et al. 2004, Penalba & Vargas 2008. As previously pointed out, the seasonal characteristics can vary significantly from one region to another (Grimm et al. 1998, Boulanger et al. 2005. These authors identified three types of precipitation regimes. The first is characterized by a minimum marked during JJA and abundant maximum during DJF, when the superficial heating together with the advection of steam in the north favors the convection associated with the South American Monsoon System that can be observed up to 20° S. The annual monsoon cycle begins in September when precipitation patterns slowly extend from the equatorial region to the south and connect with the South Atlantic Convergence Zone along the ocean. The development and connection of the South Atlantic Convergence Zone with South American Monsoon System generates a precipitation pattern along a northwesternsoutheastern axis during DJF with the highest precipitations north of 20° S (Grimm et al. 1998, Boulanger et al. 2005, Barros et al. 2006). The Pantanal floodplain plays a key role in storing rainfall in Upper Paraguay, delaying its major contributions to Paraná in almost six months (Comité Intergubernamental Coordinador de los Países de la Cuenca del Plata-CIC 2016).
On the other hand, south of 20° S, the central portion of the basin reaches its maximum at different times of the year, what suggests that  there is more than one mechanism of action, not just the forcing of monsoons. In effect, the second type of regime in the central and southern ranges of the basin tends to be distributed more evenly throughout the year. Thus, a maximum center occurs in the central range near the common border between Argentina, Brazil and Paraguay, during both transition seasons (MAM and SON) (Figure 6), accounts for a large part of the rainfall in the study domain (Rusticucci & Penalba, 2000), and relates to the mesoscale convective complexes (Velasco & Fritsch 1987). The most important contribution for the rainfall with a maximum center in this central region during winter is due to the activity of midlatitude synoptic scale systems (Vera et al. 2002). In some parts of this central region, a third peak can also be observed during the summer monsoon (third type of rainfall regime).

CONCLUSIONS
Seven interpolation models for annual and seasonal average rainfall of 34 years in the CBC and 31 years in the LPB were tested and compared. The accuracy of the interpolation was determined by the cross-validation method. The RMSE, RMSSE and ASE were chosen validation metrics and were presented in Tables II and III. Using the cross-validation method, it was observed that of the seven tested interpolation surfaces, the spherical model Cokriging (with the gaussian model) was the best precipitation interpolator in the CBC region (in LPB region). These models provided the lowest root mean square error (RMSE) and a standardized RMSE close to one for almost all cases. The estimated values from the Cokriging represent more accurately the distribution of the observed data than those from the Kriging. On the other hand, the Kriging and Cokriging methods with the gaussian model presented the worst performance for the CBC, with higher RMSE and ASE; although IDW shows a lower RMSE, the validation showed larger errors on the interpolated surface for the two regions.
It was also noted a positive impact of elevation as a predictive variable to characterize precipitation in the study domain, besides corroborating that the Geographic Information Systems are powerful tools, that allow to surpass the subjectivity of the traditional empirical method of interpolating the climate data (Andrade & Moreano 2013). In geostatistical methods, there are several possibilities for incorporating secondary data to improve primary data, such as radar, satellite, among others. However, in areas with information deficits, elevation is a widely available and accessible data, which contributes to the multivariate analysis of the rainfall (Goovaerts 2000, Vicente et al. 2003, Lloyd 2005, Atorre et al. 2007. Some authors selected the IDW and Kriging methods to describe the precipitations in Colombia, but they did not present the criteria to choose the method (Guzmán 2015, Estupiñan 2016). The results show the advantage of using multivariate geostatistical methods compared to the IDW deterministic method for annual and seasonal precipitation, when elevation data is used as a secondary variable (Goovaerts 2000, Lloyd 2005, Ly et al. 2013, providing a great benefit in improving the use of these methods to interpolate the precipitation of the two regions under study. The results in our study may be useful for scientists, engineers, hydrologists, and decision makers as a prerequisite for their use in future environmental and hydrological models for the CBC and the LPB.