Abstract

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.

Key words
Biogeographic Choco; hotspots; La Plata Basin; rainfall; South America and spatial interpolation

INTRODUCTION

Climate is one of the most important environmental factors in terrestrial ecosystems (Tuhkanen 1980TUHKANEN S. 1980. Climatic Parameters and Indices in Plant Geography. Acta Ohytogeogr Suec. Stockholm: Almqvist Wiksell International, 115 p., Chiu et al. 2009CHIU CA, LIN PH LU KC. 2009. GIS-based tests for quality control of meteorological data and spatial interpolation of climate data. Mt Res Dev 29: 339-349. doi: 10.1659/mrd.00030.). 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 2015SÁNCHEZ Ó VÉLEZ V. 2015. Complejidad de la estructura espacio-temporal de la precipitación. Rev Acad Colom Cienc 39: 304-320. doi: https://doi.org/10.18257/raccefyn.196.). 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 2000GUISAN A ZIMMERMANN NE. 2000. Predictive habitat distribution models in ecology. Ecol Model 135: 147-186. doi: https://doi.org/10.1016/S0304-3800(00)00354-9., Anderson & Martínez-Mayer 2004ANDERSON RP MARTÍNEZ-MAYER E. 2004. Modeling species’ geographic distributions for preliminary conservation assessments: an implementation with the spiny pocket mice (Heteromys) of Ecuador. Biol Conserv 16: 167-179. doi: https://doi.org/10.1016/S0006-3207(03)00187-3., Silva 2009SILVA JMC. 2009. Os Desafios Científicos para a Conservação da Biodiversidade no Brasil. Megadiversidade 5: 1-114.). In this respect, spatial interpolation techniques are a trustworthy approach to estimating climatic information from close measurements for locations without observations (Berndt & Haberlandt 2018BERNDT C HABERLANDT U. 2018. Spatial interpolation of climate variables in Northern Germany—Influence of temporal resolution and network density. J Hydrol: Regional Studies 15: 184-202. doi: 10.1016/j.ejrh.2018.02.002.).

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 2013ANDRADE L MOREANO R. 2013. Sistema de información para la interpolación de datos de temperatura y precipitación del Ecuador. Revista Politécnica 32: 70-75., Basconcillo et al. 2017BASCONCILLO JQ, DURAN GAW, FRANCISCO AA, ABASTILLAS RG, HILARIO FD, JUANILLO EL MARATAS SLA. 2017. Evaluation of Spatial Interpolation Techniques for Operational Climate Monitoring in the Philippines. Sci Online Lett Atmos 13: 114-119. doi:10.2151/sola.2017-021., Berndt & Haberlandt 2018BERNDT C HABERLANDT U. 2018. Spatial interpolation of climate variables in Northern Germany—Influence of temporal resolution and network density. J Hydrol: Regional Studies 15: 184-202. doi: 10.1016/j.ejrh.2018.02.002.). 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. 2007ATORRE F, ALFO M, DE SANCTIS M, FRANCESCONI F BRUNO F. 2007. Comparison of interpolation methods for mapping climatic and bioclimatic variables at regional scale. Int J Climatol 27: 1825-1843. doi: https://doi.org/10.1002/joc.1495.). 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. 2017BASCONCILLO JQ, DURAN GAW, FRANCISCO AA, ABASTILLAS RG, HILARIO FD, JUANILLO EL MARATAS SLA. 2017. Evaluation of Spatial Interpolation Techniques for Operational Climate Monitoring in the Philippines. Sci Online Lett Atmos 13: 114-119. doi:10.2151/sola.2017-021.). 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 1998CARUSO C QUARTA F. 1998. Interpolation methods comparison. Comput Math Appl 35: 109-126. doi: https://doi.org/10.1016/S0898-1221(98)00101-1., Ly et al. 2011LY S, CHARLES C DEGRÉ A. 2011. Geostatistical interpolation of daily rainfall at catchment scale: the use of several variogram models in the Ourthe and Ambleve catchments, Belgium. Hydrol Earth Syst Sci 15: 2259-2274. doi: doi:10.5194/hess-15-2259-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 2004POVEDA G. 2004. La hidroclimatología de Colombia: una síntesis desde la escala inter-decadal hasta la escala diurna. Rev Acad Colom Cienc 28: 201-222.) 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. 2000MYERS N, MITTERMEIER RA, MITTERMEIER CG, DA FONSECA GAB KENT J. 2000. Biodiversity hotspots for conservation priorities. Nature 403: 853-858. doi: https://doi.org/10.1038/35002501., Mittermeier et al. 2004MITTERMEIER RA, ROBLES-GIL P, HOFFMANN M, PILGRIM J, BROOKS T, MITTERMEIER C, LAMOREUX, J, DE FONSECA GAB, FORD H SELIGMANN PA. 2004. Hotspots Revisited: Earth’s Biologically Richest and Most Endangered Terrestrial Ecoregions. Mexico City: Conservation International, 392 p., 2011).

Most biodiversity hotspots are in tropical developing countries that face great challenges, such as high demographic pressure, food shortage, poverty, and corruption (Veech 2003VEECH JA. 2003. Incorporating socioeconomic factors into the analysis of biodiversity hotspots. Appl Geogr 23: 73-88. doi: https://doi.org/10.1016/S0143-6228(02)00071-1., Williams 2011WILLIAMS JN. 2011. Human population and the hotspots revisited: a 2010 assessment. In: In . Zachos FE Habel JC (Eds), Biodiversity hotspots 61-81. doi: https://doi.org/10.1007/978-3-642-20992-5_4.). 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. 2014ZEMP DC, SCHLEUSSNER CF, BARBOSA HMJ, ENT VDRJ, DONGES JF, HEINKE J, SAMAPIO E RAMMIG A. 2014. On the importance of cascading moisture recycling in South America. Atmos Chem Phys 14: 13337-13359. doi: https://doi.org/10.5194/acp-14-13337-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 rain-gauge 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 2005MITAS L MITASOVA H. 2005. Spatial Interpolation. In: Longley PA, Goodchild MF, Maguire DJ Rhind DW (Eds), Geographic Information Systems: Principles, Techniques, Management and Applications. Hoboken: J Wiley Sons, p. 481-492.). 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 1911THIESSEN AH. 1911. Precipitation averages for large areas. Mon Weather Rev 39: 1082-1089. doi: http://dx.doi.org/10.1175/1520-0493(1911)391082b:PAFLA2.0.CO;2.) and Inverse Distance Weighting (IDW) (Shepard 1968SHEPARD D. 1968. A two-dimensional interpolation function for irregularly-spaced data. In: Proceedings of the 1968 23rd Association for Computing Machinery (ACM) national conference, New York: ACM Press, p. 517-524.) or degree of smoothing, such as Radial Basis Functions (RBF), a generalized version of the multi-quadratic method (MQ) developed by Hardy (1971)HARDY RL. 1971. Multiquadric equations of topography and other irregular surfaces. J Geophys Res 76: 1905-1915. doi: https://doi.org/10.1029/JB076i008p01905.; besides the irregular triangular network (TIN) or linear regression and neural networks (Sluiter 2008SLUITER R. 2008. Interpolation methods for climate data. Literature review. Netherlands: De Bilt, 24 p.). 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. 2014CASTRO LM, GIRONÁS J FERNÁNDEZ B. 2014. Spatial estimation of daily precipitation in regions with complex relief and scarce data using terrain orientation. J Hydrol 517: 481-492. doi: https://doi.org/10.1016/j.jhydrol.2014.05.064.). However, the greater the number of predictors considered in the regression, the greater the number of stations required in the estimation (Chen & Liu 2012CHEN FW LIU CW. 2012. Estimation of the spatial rainfall distribution using inverse distance weighting (IDW) in the middle of Taiwan. Paddy Water Environ 10: 209-222. doi: 10.1007/s10333-012-0319-1.). The kriging method (Matheron 1971MATHERON GF. 1971. The theory of regionalized variables and its applications (Les cahiers du Centre de Morhologie Mathématique de Fontainebleau). Paris: École national supérieure des mines, 211 p.), creates surfaces that incorporate the statistical properties of the measured data, considering spatial autocorrelation between known points (Firdaus & Talib 2016FIRDAUS N TALIB S. 2016. Spatial Interpolation of Monthly Precipitation in Selangor, Malaysia – Comparison and Evaluation of Methods. J Appl Phys Sci 2: 1-9. doi: 10.20474/japs-2.1.1., Angulo et al. 2009ANGULO M, LOPEZ M, VICENTE S BEGUERÍA S. 2009. Mapping rainfall erosivity at a regional scale: a comparison of interpolation methods in the Ebro Basin (NE Spain). Earth Syst Sci 13: 1907-1920. doi: https://doi.org/10.5194/hess-13-1907-2009., Hao & Chang 2013HAO W CHANG X. 2013. Comparison of Spatial Interpolation Methods for Precipitation in Ningxia, China. Int J Sci Res 2: 181-184.). 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 2016FIRDAUS N TALIB S. 2016. Spatial Interpolation of Monthly Precipitation in Selangor, Malaysia – Comparison and Evaluation of Methods. J Appl Phys Sci 2: 1-9. doi: 10.20474/japs-2.1.1.).

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. 2012WAGNER PD, FIENER P, WILKEN F, KUMAR S SCHNEIDER K. 2012. Comparison and evaluation of spatial interpolation schemes for daily rainfall in data scarce regions. J Hydrol 464: 388-400. doi: https://doi.org/10.1016/j.jhydrol.2012.07.026.), multivariate methods use complementary spatial information from covariates such as elevation or dynamic radar data (Foehn et al. 2018FOEHN A, HERNÁNDEZ JG, SCHAEFLI B DE CESARE G. 2018. Spatial interpolation of precipitation from multiple rain gauge networks and weather radar data for operational applications in Alpine catchments. J Hydrol 563: 1092-1110. doi: 10.1016/j.jhydrol.2018.05.027.), since the inclusion of geographic and topographical information increases the estimation capacity of the interpolation schemes (Portalés et al. 2010PORTALÉS C, BORONAT N, PARDO-PASCUAL JE BALAGUER-BESER A. 2010. Seasonal precipitation interpolation at the Valencia region with multivariate methods using geographic and topographic information. Int J Climatol 30: 1547-1563. doi: https://doi.org/10.1002/joc.1988., Borges et al. 2016BORGES PA, FRANKE J, DA ANUNCIAÇÃO YMT, WEISS H BERNHOFER C. 2016. Comparison of spatial interpolation methods for the estimation of precipitation distribution in Distrito Federal, Brazil. Theor Appl Climatol 123: 335-348. doi: 10.1007/s00704-014-1359-9.). Studies have incorporated radar information into interpolation methods (Schiemann et al. 2011SCHIEMANN R, ERDIN R, WILLI M, FREI C, BERENGUER M SEMPERE-TORRES D. 2011. Geostatistical radar-raingauge combination with nonparametric correlograms: methodological considerations and application in Switzerland. Hydrol Earth Syst Sci 15: 1515-1536. doi: doi:10.5194/hess-15-1515-2011., Verworn & Haberlandt 2011VERWORN A HABERLANDT U. 2011. Spatial interpolation of hourly rainfall-effect of additional information, variogram inference and storm properties. Hydrol Earth Syst Sci 15: 569-584. doi: doi:10.5194/hess-15-569-2011., Wagner et al. 2012WAGNER PD, FIENER P, WILKEN F, KUMAR S SCHNEIDER K. 2012. Comparison and evaluation of spatial interpolation schemes for daily rainfall in data scarce regions. J Hydrol 464: 388-400. doi: https://doi.org/10.1016/j.jhydrol.2012.07.026.). However, these methods require a large amount of data that is often not available to users (Borges et al. 2016BORGES PA, FRANKE J, DA ANUNCIAÇÃO YMT, WEISS H BERNHOFER C. 2016. Comparison of spatial interpolation methods for the estimation of precipitation distribution in Distrito Federal, Brazil. Theor Appl Climatol 123: 335-348. doi: 10.1007/s00704-014-1359-9.).

Goovaerts (2000)GOOVAERTS P. 2000. Geostatistical approaches for incorporating elevation into the spatial interpolation of rainfall. J Hydrol 228: 113-129. doi: https://doi.org/10.1016/S0022-1694(00)00144-X. 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)PELLICONE G, CALOIERO T, MODICA G GUAGLIARDI I. 2018. Application of several spatial interpolation techniques to monthly rainfall data in the Calabria region (southern Italy). Int J Climatol 38: 3651-3666. doi: 10.1002/joc.5525. 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)ÁLZATE VDF, ARAUJO CGA, ROJAS BEO, GÓMEZ LDA MARTÍNEZ MFE. 2018. Interpolación Regnie para lluvia y temperatura en las regiones andina, caribe y pacífica de Colombia. Colombia Forestal 21: 102-118. doi: http://dx.doi.org/10.14483/2256201x.11601. 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 2005IDEAM - INSTITUTO DE HIDROLOGÍA, METEOROLOGÍA Y ESTUDIOS AMBIENTALES. 2005. Atlas Climatológico de Colombia. Bogotá: IDEAM: http://documentacion.ideam.gov.co/openbiblio/bvirtual/023777/023777.html (last access: May 7, 2019).
http://documentacion.ideam.gov.co/openbi... ) has applied the IDW method to represent the main climatic variables for the Colombian territory. On the other hand, Correa et al. (2014)CORREA CCG, TEODORO PE, CUNHA ER, OLIVEIRA-JÚNIOR JF, GOIS G, RIBEIRO LP TORRES FE. 2014. Spatial interpolation of annual rainfall in the State Mato Grosso do Sul (Brazil) using different transitive theoretical mathematical models. Int J Innov Res Sci Engineering and Technol 3: 16618-16625. doi: 10.15680/ijirset.2014.0310006. 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 1988MYERS N. 1988. Threatened biotas: ‘‘hot spots’’ in tropical forests. Environmentalist 8: 187-208. doi: https://doi.org/10.1007/BF02240252., 1990MYERS N. 1990. The biodiversity challenge: expanded hot-spots analysis. Environmentalist 10: 243-256. doi: https://doi.org/10.1007/BF02239720., Myers et al. 2000MYERS N, MITTERMEIER RA, MITTERMEIER CG, DA FONSECA GAB KENT J. 2000. Biodiversity hotspots for conservation priorities. Nature 403: 853-858. doi: https://doi.org/10.1038/35002501., Mittermeier et al. 2004MITTERMEIER RA, ROBLES-GIL P, HOFFMANN M, PILGRIM J, BROOKS T, MITTERMEIER C, LAMOREUX, J, DE FONSECA GAB, FORD H SELIGMANN PA. 2004. Hotspots Revisited: Earth’s Biologically Richest and Most Endangered Terrestrial Ecoregions. Mexico City: Conservation International, 392 p., Marchese 2015MARCHESE C. 2015. Biodiversity hotspots: A shortcut for a more complicated concept. Global Ecol Conservation 3: 297-309. doi: 10.1016/j.gecco.2014.12.008.). This region hosts approximately 3% of the global plant species, and is characterized by high rainfall, up to 8,000 mm/year (Poveda & Mesa 2000POVEDA G MESA O. 2000. On the Existence of Lloró (the Rainiest Locality on Earth): Enhanced Ocean-Land-Atmosphere Interaction by a Low-Level Jet. Geophys Res Lett 27: 1675-1678.). 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 1998GARCÍA NO VARGAS WM. 1998. The temporal climatic variability in the Rıo de La Plata Basin displayed by the river discharges. Clim Change 38: 359-379. doi: https://doi.org/10.1023/A:1005386530866., Mechoso et al. 2001MECHOSO CR ET AL. 2001. Climatology and hydrology of the Plata Basin, 2001: http://www.clivar.org/sites/default/files/documents/vamos/laplata.pdf (last access: November 26, 2019).
http://www.clivar.org/sites/default/file... ). 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 2009SAURRAL R BARROS V. 2009. Estudio de la climatología y la hidrología de la Cuenca del Plata en un conjunto de modelos climáticos globales. Meteorológica 34: 05-15.). 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. 2011MITTERMEIER RA, TURNER WR, LARSEN FW, BROOKS TM GASCON C. 2011. Global Biodiversity Conservation: The Critical Role of Hotspots. In: Zachos F Habel J (Eds), Biodiversity Hotspots. Springer, Berlin, Heidelberg, p. 3-22. doi: 10.1007/978-3-642-20992-5_1.).

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 2017NASA. 2017. Shuttle Rada Topography Mission (SRTM): http://www2.jpl.nasa.gov/srtm/ (last access: May 12, 2017).
http://www2.jpl.nasa.gov/srtm/... ). 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 2016FIRDAUS N TALIB S. 2016. Spatial Interpolation of Monthly Precipitation in Selangor, Malaysia – Comparison and Evaluation of Methods. J Appl Phys Sci 2: 1-9. doi: 10.20474/japs-2.1.1.), employs the first Law of Geography (Tobler 1970TOBLER WR. 1970. A computer movie simulating urban growth in Detroit region. Econ Geogr 46: 234-240.), 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. 2011LONGLEY P, GOODCHILD M, MAGUIRE D RHIND D. 2011. Geographic Information Systems and Science. New York: J Wiley Sons, 560 p.). That is, each climatic value at a non-sampled point z(x) is a scoring average of the distance from the values at the sampling points in the neighborhood z(x₁), z(x₂), ..., 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 2013HAO W CHANG X. 2013. Comparison of Spatial Interpolation Methods for Precipitation in Ningxia, China. Int J Sci Res 2: 181-184.).

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. 2011LY S, CHARLES C DEGRÉ A. 2011. Geostatistical interpolation of daily rainfall at catchment scale: the use of several variogram models in the Ourthe and Ambleve catchments, Belgium. Hydrol Earth Syst Sci 15: 2259-2274. doi: doi:10.5194/hess-15-2259-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. 2009ANGULO M, LOPEZ M, VICENTE S BEGUERÍA S. 2009. Mapping rainfall erosivity at a regional scale: a comparison of interpolation methods in the Ebro Basin (NE Spain). Earth Syst Sci 13: 1907-1920. doi: https://doi.org/10.5194/hess-13-1907-2009.). This continuous variable is called a regionalized variable, which consists of a derivation component and a spatially correlated random component (Burrough et al. 2015BURROUGH PA, MCDONNELL R, MCDONNELL RA LLOYD CD. 2015. Principles of geographical information systems Oxford: Oxford University Press, 371 p.). Thus, the spatially localized climatic variable z(x) is expressed by Eq. 2.

where m(x) is the derivation component, i.e. 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 1970TOBLER WR. 1970. A computer movie simulating urban growth in Detroit region. Econ Geogr 46: 234-240.). The calculated variance is a measure that determines the similarity between observations, where the highest similarity, the lower semi-variance (Lozano et al. 2004LOZANO Z, BRAVO C, OVALLES F, HERNÁNDEZ RM, MORENO B, PIÑANGO L VILLANUEVA JG. 2004. Selección de un diseño de muestreo en parcelas experimentales a partir del estudio de la variabilidad espacial de los suelos. Bioagro 16: 1-17.). Kriging focuses on the correlated spatial component and uses the adjusted semivariogram (Atorre et al. 2007ATORRE F, ALFO M, DE SANCTIS M, FRANCESCONI F BRUNO F. 2007. Comparison of interpolation methods for mapping climatic and bioclimatic variables at regional scale. Int J Climatol 27: 1825-1843. doi: https://doi.org/10.1002/joc.1495.). On the other hand, the Cokriging allows to consider the influence of external variables (co-varied, in this case height) when analyzing the cross-correlation between the errors of the different variables ɛ’₁(x), ɛ’₂(x), etc. (Angulo et al. 2009ANGULO M, LOPEZ M, VICENTE S BEGUERÍA S. 2009. Mapping rainfall erosivity at a regional scale: a comparison of interpolation methods in the Ebro Basin (NE Spain). Earth Syst Sci 13: 1907-1920. doi: https://doi.org/10.5194/hess-13-1907-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. 2019AZIZ MKBM, YUSOF F, DAUD ZM, YUSOP Z KASNO MA. 2019. Comparison of Semivariogram Models in Rain Gauge Network Design. Matematika 35: 157-170. https://doi.org/10.11113/matematika.v35.n2.1155.). This spatial statistic tool allows to express the dissimilarity as a function of distance (Sabzipour et al. 2019SABZIPOUR B, ASGHARI O SARANG A. 2019. Evaluation and optimal redesigning of river water-quality monitoring networks (RWQMN) using geostatistics approach (case study: Karun, Iran). Sustain Water Resour Manag 5: 439-455. https://doi.org/10.1007/s40899-017-0176-6.). 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.

where h is the distance, N(h) number of sample pairs within the distance h, x_i and x_i + h sampling locations separated by a distance h and, Z(x_i)and Z (x_i + h) measured values of the variables in the corresponding locations. By changing h, both in distance and direction, a set of the experimental semivariograms for the data is obtained (Aziz et al. 2019AZIZ MKBM, YUSOF F, DAUD ZM, YUSOP Z KASNO MA. 2019. Comparison of Semivariogram Models in Rain Gauge Network Design. Matematika 35: 157-170. https://doi.org/10.11113/matematika.v35.n2.1155.).

For cokriging, Goovaerts (2000)GOOVAERTS P. 2000. Geostatistical approaches for incorporating elevation into the spatial interpolation of rainfall. J Hydrol 228: 113-129. doi: https://doi.org/10.1016/S0022-1694(00)00144-X. 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. 2011LY S, CHARLES C DEGRÉ A. 2011. Geostatistical interpolation of daily rainfall at catchment scale: the use of several variogram models in the Ourthe and Ambleve catchments, Belgium. Hydrol Earth Syst Sci 15: 2259-2274. doi: doi:10.5194/hess-15-2259-2011., Goovaerts 1997GOOVAERTS P. 1997. Geostatistics for natural resources evaluation. New York: Oxford University Press, 483 p.). 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 2007BOHLING GC. 2007. Introduction to geostatistics. Kansas Geological Survey Open File Report no. 26 (2007), 50 p.). 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 2013MAZZELLA A MAZZELLA A. 2013. The importance of the model choice for experimental Semivariogram modeling and its consequence in evaluation process. J Eng (United Kingdom) 2013: 1-11. doi: https://doi.org/10.1155/2013/960105., Nadiah et al. 2016NADIAH N, FIRDAUS M TALIB SA. 2016. Spatial interpolation of monthly precipitation in selangor, Malaysia – comparison and evaluation of methods. J Appl Phys Sci 2: 346-357. doi: https://doi.org/10.20474/-japs2.1.1., Aziz et al. 2019AZIZ MKBM, YUSOF F, DAUD ZM, YUSOP Z KASNO MA. 2019. Comparison of Semivariogram Models in Rain Gauge Network Design. Matematika 35: 157-170. https://doi.org/10.11113/matematika.v35.n2.1155.). 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)CASTRO FS, PEZZOPANE JEM, CECÍLIO RA, PEZZOPANE JRM XAVIER AC. 2010. Avaliação do desempenho dos diferentes métodos de interpoladores para parâmetros do balanço hídrico climatológico. Rev Bras Enge Agrí Amb 14: 871-880. doi: https://doi.org/10.1590/s1415-43662010000800012. 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 1986MCBRATNEY A WEBSTER R. 1986. Chosing functions for semivariograms of soil properties and fitting them to sampling estimates. J Soil Sci 37: 617-639. doi: https://doi.org/10.1111/j.1365-2389.1986.tb00392.x., Goovaters 1997, 2000, Deutsch & Journel 1998DEUTSCH C JOURNEL A. 1998. Gslib: Geostatistical software library and user’s guide. New York: Oxford University Press, 147 p., Ly et al. 2011LY S, CHARLES C DEGRÉ A. 2011. Geostatistical interpolation of daily rainfall at catchment scale: the use of several variogram models in the Ourthe and Ambleve catchments, Belgium. Hydrol Earth Syst Sci 15: 2259-2274. doi: doi:10.5194/hess-15-2259-2011., Mazzella & Mazzella 2013MAZZELLA A MAZZELLA A. 2013. The importance of the model choice for experimental Semivariogram modeling and its consequence in evaluation process. J Eng (United Kingdom) 2013: 1-11. doi: https://doi.org/10.1155/2013/960105., Aziz et al. 2019AZIZ MKBM, YUSOF F, DAUD ZM, YUSOP Z KASNO MA. 2019. Comparison of Semivariogram Models in Rain Gauge Network Design. Matematika 35: 157-170. https://doi.org/10.11113/matematika.v35.n2.1155.), 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)LY S, CHARLES C DEGRÉ A. 2011. Geostatistical interpolation of daily rainfall at catchment scale: the use of several variogram models in the Ourthe and Ambleve catchments, Belgium. Hydrol Earth Syst Sci 15: 2259-2274. doi: doi:10.5194/hess-15-2259-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 2006CHANG K. 2006. Introduction to Geographic Information Systems. Boston: McGraw-Hill Higher Education, 419 p.): (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 2006CHANG K. 2006. Introduction to Geographic Information Systems. Boston: McGraw-Hill Higher Education, 419 p.). Additionally, the Average Standard Error (ASE, Eq. 7) was calculated to assess the variability of predictions:

where N is the number of values in the dataset and $${\sigma }^2$$ is the interpolation variance for the location x_i (Robinson & Metternicht 2006ROBINSON TP METTERNICHT G. 2006. Testing the performance of spatial interpolation techniques for mapping soil properties. Computers and electronics in agriculture 50: 97-108. doi: https://doi.org/10.1016/j.compag.2005.07.003.).

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 RMSE, the variability of the predictions is overestimating (underestimating) (Kalivas et al. 2013KALIVAS DP, KOLLIAS VJ APOSTOLIDIS EH. 2013. Evaluation of three spatial interpolation methods to estimate forest volume in the municipal forest of the Greek island Skyros. Geo-Spatial Information Science 16: 100-112. doi: 10.1080/10095020.2013.766398.).

RESULTS AND DISCUSSIONS

Colombian Biogeographic Choco

Table II shows the RMSE, RMSSE and ASE values ​​for the seven interpolation surfaces tested. In the case of the Cokriging model, the elevation was used as a second variable because of its influence on rainfall. Although the IDW shows lower values ​​of RMSE during December-January-February (DJF), in space considerable failures were observed in the interpolated surface, which leads to discarding this method. Based on Table II, for exponential and spherical cokriging models, RMSSE values are closer to 1, while kriging results show lower RMSSE for spherical and gaussian models, and gaussian cokriging and exponential kriging models show values well above 1, which means that the exponential and spherical cokriging models have a better ability to represent the variability of the dataset. In addition, since ASE values must be around 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 2013GUARÍN GW POVEDA G. 2013. Variabilidad Espacial y Temporal del Almacenamiento de Agua en el Suelo en Colombia. Rev Acad Colomb Cienc 37: 89-113.). 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 1999POVEDA G MESA O. 1999. La corriente de chorro superficial del Oeste (“del Chocó”) y otras dos corrientes de chorro en Colombia: climatología y variabilidad durante las fases del ENSO. Rev Acad Colom Cienc 23: 517-528., Trenberth 1999TRENBERTH KE. 1999. Atmospheric moisture recycling: Role of advection and local evaporation. J Climate 12: 1368-1381. doi: https://doi.org/10.1175/1520-0442(1999)0121368:AMRROA2.0.CO;2.). 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 west-east 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 1994ESLAVA J. 1994. Climatología del Pacífico Colombiano. Academia Colombiana de Ciencias Geofísicas. Bogotá: Colección Eratóstenes., Poveda 2004POVEDA G. 2004. La hidroclimatología de Colombia: una síntesis desde la escala inter-decadal hasta la escala diurna. Rev Acad Colom Cienc 28: 201-222., Guzmán et al. 2014GUZMÁN D, RUÍZ JF CADENA M. 2014. Regionalización de colombia según la estacionalidad de la precipitación media mensual, a través Análisis de Componentes Principales (ACP). Grupo de Modelamiento de Tiempo, Clima y Escenarios de Cambio Climático. Bogotá: http://www.ideam.gov.co/documents/21021/21141/Regionalizacion+de+la+Precipitacion+Media+Mensual/1239c8b3-299d-4099-bf52-55a414557119 (last access: November 26, 2019).
http://www.ideam.gov.co/documents/21021/... ). 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)POVEDA MI, ROJAS C, RUDAS A RANGEL JO. 2004. El Chocó Biogeográfico: Ambiente físico. In: Rangel JO (Ed), Diversidad Biótica IV. El chocó Biogeográfico/Costa Pacífica, Bogotá: Universidad Nacional de Colombia, Instituto de Ciencias Naturales and Conservación Internacional, p. 1-21., 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 sub-regional scale Hurtado (2009)HURTADO AF. 2009. Estimación de los campos mensuales históricos de precipitación en el territorio colombiano. MSc thesis. Universidad Nacional de Colombia Sede Medellín: http://bdigital.unal.edu.co/2305/1/71364267.2009.pdf (last access: November 26, 2019). (Unpublished).
http://bdigital.unal.edu.co/2305/1/71364... 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 2000POVEDA G MESA O. 2000. On the Existence of Lloró (the Rainiest Locality on Earth): Enhanced Ocean-Land-Atmosphere Interaction by a Low-Level Jet. Geophys Res Lett 27: 1675-1678., Zea et al. 2000ZEA JA, LEÓN GE ESLAVA JA. 2000. Influencia de sistemas extratropicales en Colombia. Meteorol Colomb 1: 31-40.). 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. 2002POVEDA G, VÉLEZ JI, MESA O, HOYOS C, MEJÍA J, BARCO OJ CORREA PL. 2002. Influencia de fenómenos macroclimáticos sobre el ciclo anual de la hidrología Colombiana: Cuantificación lineal, no lineal y percentiles probabilísticos. Meteorol Colomb 121-130., 2011POVEDA G, ÁLVAREZ DM RUEDA ÓA. 2011. Hydro-climatic variability over the Andes of Colombia associated with ENSO: A review of climatic processes and their impact on one of the Earth’s most important biodiversity hotspots. Clim Dyn 36: 2233-2249. doi: https://doi.org/10.1007/s00382-010-0931-y., 2014POVEDA G, JARAMILLO L VALLEJO LF. 2014. Seasonal precipitation patterns along pathways of South American low-level jets and aerial rivers. Water Resour Res 50: 98-118. https://doi.org/10.1002/2013WR014087., Trojer 2018TROJER H. 2018. Meteorología y climatología de la vertiente del Pacífico colombiano. Rev Acad Colomb Ciencias Exactas, Físicas y Nat 41: 467. https://doi.org/10.18257/raccefyn.583.). Besides, the differences in the Andean Region have been previously documented by Sedano (2017)SEDANO CRK. 2017. Influencia de la variabilidad climática y factores antrópicos en los extremos hidrológicos en el Valle Alto del río Cauca, Colombia. PhD dissertation. Universistat Polotècnica de València: https://riunet.upv.es/handle/10251/90579# (last access: December 11, 2019).
https://riunet.upv.es/handle/10251/90579... , 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.

La Plata Basin

Table III presents the RMSE, RMSSE and ASE values for the seven tested interpolation surfaces in the LPB. The results of the seasonal and annual RMSSE indicate that gaussian cokriging and, gaussian and exponential kriging have the best results with values close to 1 in most cases. Exponential (spherical) cokriging presented very high RMSSE values for DJF and MAM (SON). Low RMSSE values are observed in MAM and JJA for spherical kriging. Based on the differences between ASE and RMSE, the best prediction of variability are presented for gaussian cokriging and exponential kriging, 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. 2005BOULANGER JP, LELOUP J, PENALBA O, RUSTICUCCI M, LAFON F VARGAS W. 2005. Observed precipitation in the Paraná-Plata hydrological basin: long-term trends, extreme conditions and ENSO teleconnections. Clim Dynam 24: 393-413. doi: 10.1007/s00382-004-0514-x., Barros et al. 2006BARROS V, CLARKE R DIAS PS. 2006. Climate change in the La Plata basin. São José dos Campos: Publication of the Inter-American Institute for Global Change Research (IAI): http://www-atmo.at.fcen.uba.ar/cordex/Barrosetal2006ChapterII.pdf (last access: November 26, 2019).
http://www-atmo.at.fcen.uba.ar/cordex/Ba... , Bidegain et al. 2017BIDEGAIN M, PASTEN M, NAGY G, CORONEL G, FERRER J ARANA I. 2017. Tendencias recientes de las precipitaciones e impactos asociados con ENSO en la cuenca del Río de la Plata. Paraquaria Nat 5: 8-18. https://paraquaria.org.py/product/tendencias-recientes-de-las-precipitaciones-e-impactos-asociados-con-enso-en-la-cuenca-del-rio-de-la-plata/ (last access: November 26, 2019)., among others). 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. 2009ESPINOZA VJC, RONCHAIL J, GUYOT JL, COCHONNEAU G, NAZIANO F, LAVADO W, DE OLIVEIRA E, POMBOSA R VAUCHEL P. 2009. Spatio-temporal rainfall variability in the Amazon basin countries (Brazil, Peru, Bolivia, Colombia, and Ecuador). Int J Climatol 29: 1574-1594. doi: 10.1002/joc.1791.) and in Bolivia (Ronchail & Gallaire 2006RONCHAIL J GALLAIRE R. 2006. ENSO and rainfall along the Zongo valley (Bolivia) from the Altiplano to the Amazon basin. Int J Climatol 26: 1223-1236. doi: 10.1002/joc.1296.). 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).

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 2002BERBERY EH BARROS VR. 2002. The Hydrologic Cycle of the La Plata Basin in South America. J Hydrometeorol 3: 630-645. doi: https://doi.org/10.1175/1525-7541(2002)0030630:thcotl2.0.co;2., Vera et al. 2002VERA CS, VIGLIAROLO PK BERBERY EH. 2002. Cold season synoptic-scale waves over subtropical South America. Mon Weather Rev 130: 684-699. doi: https://doi.org/10.1175/1520-0493(2002)1300684:CSSSWO2.0.CO;2., Gan et al. 2004GAN MA, KOUSKY VE ROPELEWSKI CF. 2004. The South America Monsoon circulation and its relationship to rainfall over west-central Brazil. J Clim 17: 47-66. doi: https://doi.org/10.1175/1520-0442(2004)0170047:TSAMCA2.0.CO;2., Penalba & Vargas 2008PENALBA OC VARGAS WM. 2008. Variability of low monthly rainfall in La Plata Basin. Meteorol Appl 323: 313-323. https://doi.org/10.1002/met.).

As previously pointed out, the seasonal characteristics can vary significantly from one region to another (Grimm et al. 1998GRIMM AM, FERRAZ SE GOMES J. 1998. Precipitation anomalies in southern Brazil associated with El Niño and La Niña events. J Climate 11: 2863-2880. doi: https://doi.org/10.1175/1520-0442(1998)011%3C2863:PAISBA%3E2.0.CO;2., Boulanger et al. 2005BOULANGER JP, LELOUP J, PENALBA O, RUSTICUCCI M, LAFON F VARGAS W. 2005. Observed precipitation in the Paraná-Plata hydrological basin: long-term trends, extreme conditions and ENSO teleconnections. Clim Dynam 24: 393-413. doi: 10.1007/s00382-004-0514-x.). 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 northwestern-southeastern axis during DJF with the highest precipitations north of 20° S (Grimm et al. 1998GRIMM AM, FERRAZ SE GOMES J. 1998. Precipitation anomalies in southern Brazil associated with El Niño and La Niña events. J Climate 11: 2863-2880. doi: https://doi.org/10.1175/1520-0442(1998)011%3C2863:PAISBA%3E2.0.CO;2., Boulanger et al. 2005BOULANGER JP, LELOUP J, PENALBA O, RUSTICUCCI M, LAFON F VARGAS W. 2005. Observed precipitation in the Paraná-Plata hydrological basin: long-term trends, extreme conditions and ENSO teleconnections. Clim Dynam 24: 393-413. doi: 10.1007/s00382-004-0514-x., Barros et al. 2006BARROS V, CLARKE R DIAS PS. 2006. Climate change in the La Plata basin. São José dos Campos: Publication of the Inter-American Institute for Global Change Research (IAI): http://www-atmo.at.fcen.uba.ar/cordex/Barrosetal2006ChapterII.pdf (last access: November 26, 2019).
http://www-atmo.at.fcen.uba.ar/cordex/Ba... ). 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 2016COMITÉ INTERGUBERNAMENTAL COORDINADOR DE LOS PAÍSES DE LA CUENCA DEL PLATA - CIC. 2016. CIC Cuenca del Plata. http://www.cicplata.org/ (last access: November 26, 2019).
http://www.cicplata.org/... ).

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, 2000RUSTICUCCI M PENALBA O. 2000. Interdecadal changes in the precipitation seasonal cycle over Southern South America and their relationship with surface temperature. Clim Res 16: 1-15. doi: 10.3354/cr016001.), and relates to the mesoscale convective complexes (Velasco & Fritsch 1987VELASCO I FRITSCH JM. 1987. Mesoscale convective complexes in the Americas. J Geophys Res-Atmos 92: 9591-9613. doi: https://doi.org/10.1029/JD092iD08p09591.). 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. 2002VERA CS, VIGLIAROLO PK BERBERY EH. 2002. Cold season synoptic-scale waves over subtropical South America. Mon Weather Rev 130: 684-699. doi: https://doi.org/10.1175/1520-0493(2002)1300684:CSSSWO2.0.CO;2.). 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 2000GOOVAERTS P. 2000. Geostatistical approaches for incorporating elevation into the spatial interpolation of rainfall. J Hydrol 228: 113-129. doi: https://doi.org/10.1016/S0022-1694(00)00144-X., Vicente et al. 2003VICENTE SM, SAZ MA CUADRAT JM. 2003. Comparative analysis of interpolation methods in the middle Ebro Valley (Spain): application to annual precipitation and temperature. Clim Res 24: 161-180. doi: doi:10.3354/cr024161., Lloyd 2005LLOYD C. 2005. Assessing the effect of integrating elevation data into the estimation of monthly precipitation in Great Britain. J Hydrol 308: 128-150. doi: doi:10.1016/j.jhydrol.2004.10.026., Atorre et al. 2007ATORRE F, ALFO M, DE SANCTIS M, FRANCESCONI F BRUNO F. 2007. Comparison of interpolation methods for mapping climatic and bioclimatic variables at regional scale. Int J Climatol 27: 1825-1843. doi: https://doi.org/10.1002/joc.1495.). 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 2015GUZMÁN A. 2015. Análisis del comportamiento espacio temporal de la precipitación en la región Pacífica Colombiana en presencia de los fenómenos El Niño y La Niña, en el periodo de 1983-2009. Dissertation. Universidad de Nariño: http://sired.udenar.edu.co/1584/ (last access: May 5, 2018). (Unpublished)., Estupiñan 2016ESTUPIÑAN A. 2016. Estudio de la variabilidad espacio temporal de la precipitación en Colombia. PhD dissertation. Universidad Nacional de Colombia: http://bdigital.unal.edu.co/54014/1/1110490004.2016.pdf (last access: November 26, 2019).
http://bdigital.unal.edu.co/54014/1/1110... ). 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 2000GOOVAERTS P. 2000. Geostatistical approaches for incorporating elevation into the spatial interpolation of rainfall. J Hydrol 228: 113-129. doi: https://doi.org/10.1016/S0022-1694(00)00144-X., Lloyd 2005LLOYD C. 2005. Assessing the effect of integrating elevation data into the estimation of monthly precipitation in Great Britain. J Hydrol 308: 128-150. doi: doi:10.1016/j.jhydrol.2004.10.026., Ly et al. 2013LY S, CHARLES C DEGRÉ A. 2013. Different methods for spatial interpolation of rainfall data for operational hydrology and hydrological modeling at watershed scale: a review. Biotechnol Agron Soc 17: 392-406.), 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.

ACKNOWLEDGMENTS

The first author was supported by the Doctoral Scholarship of the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brazil (CAPES) - Finance Code 001 and Universidad del Valle (Cali-Colombia). The second and third authors were supported by the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) of Brazil under grants 305611/2019-4 and 302322/2017-5, respectively. The fifth author was supported by the CEIBA Foundation for your doctoral studies. We thank the Water Resources and Soil Engineering Research Group (IREHISA), Universidad del Valle, Cali, Colombia. The present research article is the result of the doctoral project: “Study of the space-time variability of precipitation associated with the transport of moisture from low-level jets in South America”. The authors thank the anonymous reviewers for their useful suggestions.

REFERENCES

Publication Dates

History