Performance Evaluation of Bias Correction Methods for Climate Change Monthly Precipitation Projections over Costa Rica

Abstract

Six bias correction (BC) methods; delta-method (DT), linear scaling (LS), power transformation of precipitation (PTR), empirical quantile mapping (EQM), gamma quantile mapping (GQM) and gamma-pareto quantile mapping (GPQM) were applied to adjust the biases of historical monthly precipitation outputs from five General Circulation Models (GCMs) dynamically downscaled by two Regional Climate Models (RCMs) for a total of seven different GCM-RCM pairs over Costa Rica. High-resolution gridded precipitation observations were used for the control period 1951–1980 and validated over the period 1981–1995. Results show that considerable biases exist between uncorrected GCM-RCM outputs and observations, which largely depend on GCM-RCM pair, seasonality, climatic region and spatial resolution. After the application of bias correction, substantial biases reductions and comparable performances among most BC methods were observed for most GCM-RCM pairs; with EQM and DT marginally outperforming the remaining methods. Consequently, EQM and DT were selectively applied to correct the biases of precipitation projections from each individual GCM-RCM pair for a near-future (2011–2040), mid-future (2041–2070) and far-future (2071–2100) period under Representative Concentration Pathways (RCPs) 2.6, 4.5 and 8.5 using the control period 1961–1990. Results from the bias-corrected future ensemble-mean anticipate a marked decreasing trend in precipitation from near to far-future periods during the dry season (December, January, February (DJF) and March, April, May (MAM)) for RCP4.5 and 8.5; with pronounced drier conditions for those climatic regions draining towards the Pacific Ocean. In contrast, mostly wetter conditions are expected during the dry season under RCP2.6, particularly for the Caribbean region. In most of the country, the greatest decrease in precipitation is projected at the beginning of the rainy season (June, July, August (JJA)) for the far-future period under RCP8.5, except for the Caribbean region where mostly wetter conditions are anticipated. Regardless of future period, slight increases in precipitation with higher radiative forcing are expected for SON excluding the Caribbean region, where precipitation is likely to increase with increasing radiative forcing and future period. This study demonstrates that bias correction should be considered before direct application of GCM-RCM precipitation projections over complex territories such as Costa Rica.

1. Introduction

Global warming caused by anthropogenic increases of greenhouse gases concentrations in the atmosphere is associated with climate change around the world [1]. Higher temperatures will increase atmospheric moisture holding capacity, which can lead to increases in precipitation [2]. Changes in precipitation patterns will alter the world’s terrestrial hydrological systems, strongly influencing water resources availability, which will directly impact natural ecosystems along with critical socioeconomic sectors such as water supply, energy, crop production, animal raising and fishery [3]. Furthermore, rises in temperature along with future changes in precipitation patterns will be highly variable at regional and local scales, with some regions becoming drier while other regions are becoming wetter [4]. Managing water resources at local and regional scales under changing climate is a major challenge faced by societies today, especially due to considerable uncertainties regarding future atmospheric greenhouse gas emissions scenarios [5,6]. Central America is a region of the world that could increase its vulnerability to extreme climatic events in the future [7]. The region has been identified as a climate change “hot-spot,” where decreasing trends in precipitation have already been highlighted [8,9]. Central America is well-known for its climatic fluctuations including major droughts [10], impact of hurricanes [11] and increasing temperature trends [12] and it is expected to be the most responsive tropical region to future climate change [13]. Moreover, future changes in precipitation extremes and characteristics are among the most critical factors determining the overall impact of climate change on the hydrology and the water resources of the Central American region [14].

Relevant sectors such as hydropower generation, food production and cattle ranching are especially vulnerable to future climate extremes due to the strong dependence on the availability and distribution of water resources in Central America, especially for countries like Costa Rica and Panama [15,16]. Therefore, there is a need to produce reliable high-quality projections of future changes in precipitation due to climate change in the region.

The General Circulation Models (GCMs) developed in the fifth phase of the Coupled Model Intercomparison Project (CMIP5) of the Intergovernmental Panel on Climate Change (IPCC) are the most advanced tools currently available for understanding the mechanisms of past climate, investigating the relative contributions of anthropogenic and natural drivers and projecting possible future climate changes under idealized emission or concentration assumptions [17,18]. General Circulation Models (GCMs) are mathematical representations of physical processes in the atmosphere, ocean and land surface, which generate meteorological variables such as precipitation, temperature, wind speed, relative humidity and solar radiation among others [19]. GCMs meteorological variables are generated by solving the primitive equations of thermodynamics, mass and momentum [20]. However, GCMs outputs remain relatively coarse in terms of their spatial and temporal resolutions [21]. Systematic biases caused by imperfect conceptualization, spatial averaging and non-inclusion of local features make the direct use of GCMs outputs unreliable to force fine-scale impact assessment models [22]. Furthermore, GCMs can only provide limited spatial details in regions that feature narrow continental extensions, high mountains and steep slopes because the altitude-dependence of precipitation is not properly represented [23]. Specifically, most of the CMIP5 GCM models show less skill in precipitation simulation over regions with complex topography [24]. The Central America region for instance, has considerable topographic relief, which implies the existence of large gradients in climate variables such as precipitation, temperature, humidity and wind [25,26]. The interactions between the complex topography, narrow land configuration and climate patterns determined by the Pacific and Atlantic oceans give rise to numerous climate and ecological zones both horizontally and vertically that no single coarse-resolution GCM is capable of resolving, due in part to locally intense convective activity [27,28]. Moreover, GCMs do not provide reliable information for most hydrological relevant variables (e.g., precipitation and temperature) on resolutions below 200 km at monthly or daily scales, which is notably crude compared to that of hydrologic processes [29].

In light of the abovementioned constraints, higher temporal and spatial resolution precipitation projections that better represent key aspects of weather systems as demanded by local climate change studies must be produced to overcome the limitations of GCMs outputs, especially for geographically complex regions such as Central America. One approach to bridge this gap is the use of dynamical downscaling, which preserves the physical coherence of GCMs outputs [30,31].

Dynamical downscaling techniques transfer information from GCMs to finer scales by applying higher resolution Regional Climate Models (RCMs) over a limited area with initial boundary conditions taken from a driving GCM [32]. The RCMs added value is related to their ability to providing additional climate change signals that are not resolved in the coarser resolution GCMs [33,34]. Compared to GCMs, RCMs are more reliable at reproducing relevant patterns of local precipitation, since they explicitly solve mesoscale atmospheric processes so that the topographic effects on precipitation can be much better represented at a regional scale [35].

RCMs future projections of precipitation have extensively been used to assess the impact of climate change on various water resources systems including aquifer recharge [36,37], hydrological modelling [38,39,40], extreme river flows and flood analysis [30], water infrastructure systems design [41,42], land use and forestry [43,44] and tracking of tropical cyclones [45,46]. The improved performance of RCMs over GCMs is mainly due to (1) better representation of geographical features such as orography related to higher spatial resolutions (12.5–50 km) and (2) better description of physical processes by means of sub-grid scale parameterization and more detailed land surface schemes, providing spatially and physically consistent outputs for operational use on climate timescales [47,48]. The use of RCMs is preferable since they have proven to provide more reliable results for impact study of climate change on regional and local scales than GCMs [49]. However, RCMs model outputs may still produce considerable systematic biases; hence their direct use as input for impact models is inadvisable since they might lead to inaccurate conclusions [50,51,52]. Biases are defined as long term average deviations between simulations and observations [53].

RCMs biases originate from various sources including (1) boundary errors transferred from GCMs to RCMs, (2) insufficiently resolved surface properties related to spatial resolution (3) errors due to numeric resolutions and (4) internal model parameterization [54,55]. RCMs biases also depend on GCMs simulated variables, which could be significant for precipitation due to its highly nonlinear nature and large spatial variability, making them highly dependent on model resolution [56,57]. RCMs accuracy is seasonally and regionally dependent, showing diminished accuracy at representing summer precipitation than autumn and winter precipitation due to difficulties in modelling convective rainfall [58]. Furthermore, RCMs perform better in regions with temperate climate conditions than tropical regions, since tropical rainfall is often predominantly convective in nature and occurs at a sub daily time scale [59]. On the other hand, RCMs local feedbacks and forcings can lead to significantly different climate change responses in different parts of the world [39]. For instance, to address the uncertainty associated with model outputs, it is advisable to consider the use of multiple GCM-RCMs and different emission scenarios [60,61]. For scenario uncertainty, the spread across the Representative Concentration Pathway (RCP) emission scenarios, as proposed in the Fifth Assessment Report (AR5) by Intergovernmental Panel on Climate change (IPCC) [7], is the dominant source of uncertainty for projections to the end of the century [62]. Therefore, even though RCMs considerably improve the temporal and spatial resolutions of GCMs outputs, there is still a scale gap between the simulated grid cell values of precipitation and point-scale data [63].

As a consequence, postprocessing of RCMs outputs through the application of bias correction (BC) methods is generally a prerequisite step for most climate change impact studies [64,65]. Bias correction is the adjustment of biased simulated data to observations [66]. BC methods aim to add value to RCMs outputs by removing systematic biases of simulated data so that they can be used in climate change impact modelling [67]. Several bias correction methods have been developed to adjust meteorological variables from the RCMs ranging from simple scaling approaches to sophisticated distribution mapping [58,68,69]. BC methods are not expected to modify the original climate change signals [34] but most of them assume that the statistical relationship between the historic model simulation and observations will remain stable in the future period. Nevertheless, it has been found that this assumption may not be true because of the nonstationarity of the climate change processes [15]. Even though BC methods are capable of reducing biases from RCM outputs [55], their performances are most likely regionally dependent and should be evaluated and validated over a recent-past period prior to any climate change application [21]. Bias correction of precipitation is equally sensitive to the selection of a particular time period, so a good temporal structure of original observations is imperative [38,70]. Located in the Central America region, this study focuses on Costa Rica where most of the territory has experienced some degree of precipitation change during the period 1961–1990, showing increases on the north-western Caribbean side and decreases on the Pacific side [12,71]. As well as its vulnerability to perturbations and fluctuations in precipitation associated with climate change and its importance to ecology and relevance for the hydrological cycle, the performance of six bias correction methods, used to adjust the biases of two RCMs driven by five GCMs is evaluated over the recent past period 1951–1995.

Since the applicability of bias correction methods has not been sufficiently investigated in Costa Rica, this study aims at answering the following research questions regarding precipitation: (1) Can significant biases be identified between uncorrected GCM-RCMs outputs and observations? (2) Is the performance of the various bias correction methods spatially, regionally and seasonally dependent? (3) What are the best bias correction methods for the selected GCM-RCMs? (4) What are the projected bias-corrected impacts of climate change over each climatic region of Costa Rica considering the various RCPs emission scenarios and future periods considered? Projecting the impact of future climate change on precipitation and its effects on water resources is of paramount importance to the sustainable development of the country.

2. Materials and Methods

2.1. Study Area

Costa Rica is located across the Central American isthmus between Panama and Nicaragua and occupies an area of 51,060 km² (Figure 1b). The country is bordered by the Caribbean Sea to the east and the Pacific Ocean to the west, which favors oceanic and climatological influences from both oceans. Precipitation variability in Costa Rica is driven by interactions between the local topography and a combination of the seasonal migration within the Intertropical Convergence Zone (ITCZ), which includes sea breeze effects, monsoonal circulations, strong easterly trade winds such as the Caribbean Low-Level Jet (CLLJ) [72], cold air masses from high and mid-latitudes during winter and perturbing influences of hurricanes and tropical cyclones in the Atlantic Ocean, all of which promotes a strong seasonal cycle and regional variability [73]. Moreover, the Central American isthmus is surrounded by one of the warmest body of waters known as the Western Hemisphere warm pool (WHWP) that covers the eastern tropical Pacific, the Gulf of Mexico, the Caribbean Sea and the western tropical North Atlantic [74]. The seasonal cycle of sea surface temperatures (SSTs) within this region affects the annual cycle of precipitation, convective activity and storm development [75]. Additionally, the temporal and spatial variability of precipitation in Costa Rica is heavily influenced by El Niño-Southern Oscillation (ENSO), which complex warm or wet responses vary in terms of their signs, magnitudes, duration and seasonality between those areas draining towards the Pacific and those draining towards the Caribbean [76].

On the other hand, the Costa Rica territory is meridionally divided by northwest-southeast trending cordilleras of high complexity, characterized by mountainous topographies and elevations ranging from sea-level to around 3400 m (Figure 1). Since the climate in Costa Rica is equally dependent on geographical altitude, it has been divided into six separate climatic regions (North-Pacific, Central-Pacific, South-Pacific, Central-Valley, North and Caribbean), which domains are determined by the position and elevation of the abovementioned cordilleras (Figure 1). These climatic regions act as individual domains and reflect the separation that elevation imposes between Caribbean (northeastern) and Pacific (southwestern) sources of moisture, revealing the important role played by elevation in controlling the circulation patterns of local precipitation [77]. The northeastern side of the country covers the Caribbean and North climatic regions.

Runoff from the North region flows towards the San Juan River, which represents the border between Costa Rica and Nicaragua and ultimately discharges in the Caribbean. The Caribbean region as its names suggests, drains directly towards the Caribbean Sea (Figure 1). Precipitation along these two climatic regions is year-round persistent as accumulated precipitation averages over 100 mm/month during December, January, February (DJF), even though the rest of the country experiences the driest period during this season (Figure 2). Peaks for the Caribbean and North regions are present in summer (June, July, August (JJA)) and autumn (September, October, November (SON)), when the northeast trade winds are strongest and easterly wave and tropical cyclone activity is most frequent [78]. Conversely, spring (March, April, May (MAM)) exhibits a transitional precipitation pattern from drier to wetter conditions, with May already approaching the summer peak. Precipitation during DJF is mainly associated with cold air outbreaks and frontal systems that pick up moisture over the Caribbean and advect precipitation on the North and Caribbean regions [79].

The southwestern side of the country on the other hand, includes the North-Pacific, Central-Pacific and South-Pacific regions, which directly drain towards the Pacific Ocean; whereas the Central-Valley region drains into the Tárcoles River, which ultimately discharges into the Pacific Ocean (Figure 1). Precipitation on the southwestern Pacific coast peaks during summer (JJA) and autumn (SON), as sea surface temperature (SST) also peaks and the ITCZ progresses through the region (Figure 2). During JJA, SSTs in the WHWP get considerably higher (above 28.5 C), resulting in intense convection and precipitation [80]. Atlantic tropical cyclones also drive onshore winds on the west side of the cordilleras during late summer (August) and early autumn (SON), which promote heavy precipitation events along the Pacific coast [71]. Precipitation on the southwestern Pacific coast considerably decreases during boreal winter (DJF) when surrounding SSTs are lowest and trade winds are weaker with a transitional dry to wet pattern during MAM, being May the wettest month during spring (Figure 2). There is also evidence of the Mid-Summer Drought (MSD) that essentially impacts the regions on the southwestern Pacific coast (Figure 2) and which manifests as a relative minimum within the overall JJA precipitation peak [81].

2.2. Datasets and Experimental Setup

A 1x1 km gridded monthly precipitation climatology based on rain-gauge historical observations (Figure 1) provided by Instituto Meteorológico of Costa Rica (IMN) [82,83] was selected as the observational dataset for the period 1951–1995. The Global Precipitation Climatology Centre (GPCC) [84] dataset was used on a comparative basis to evaluate the overall spatial and temporal coherence of the locally generated IMN climatology. Monthly precipitation totals from five General Circulation Models (GCMs), dynamically downscaled by the RCA4 Regional Climate Model [85] available in the context of the Coordinated Regional Climate Downscaling Experiment (CORDEX) over Central America (CA) were used in this study (

Table 1. General Circulation Models (GCMs) and Regional Climate Models (RCMs) configuration. Included Representative Concentration Pathways (RCPs) are marked with a X symbol, whereas excluded RCPs are marked with an O symbol.

Driving GCM-Model	GCM Model-Center	RCM Model-Center	RCM Resolution	RCP2.6	RCP 4.5	RCP 8.5
PRECIS	Met Office Hadley Centre	HadRM3P Met Office Hadley Centre	0.22° × 0.22°	X	X	X
CanESM2	Canadian Centre for Climate Modelling and Analysis	RCA4 Swedish (SHMI)	0.44° × 0.44°	O	O	X
EC-EARTH	Irish Centre For High-End Computing	RCA4 Swedish (SHMI)	0.44° × 0.44°	X	X	X
HadGEM2-ES	Met Office Hadley Centre	RCA4 Swedish (SHMI)	0.44° × 0.44°	X	X	X
MIROC5	Tokyo Center for Climate System Research	RCA4 Swedish (SHMI)	0.44° × 0.44°	X	O	X
MPI-ESM-LR	Max Planck Institute for Meteorology	RCA4 Swedish (SHMI)	0.44° × 0.44°	X	X	X

).

Data from CORDEX-CA models, which covers the period 1950 to 2100 have a 0.44° × 0.44° (~50 km) spatial resolution and its domain aims to adequately simulate key large-scale circulation features that affect regional climate in Central America including the Intertropical Convergence Zone (ITCZ) displacements, the Caribbean Low-Level Jet (CLLJ), the occurrence of northerly cold fronts and the seasonal changes of the Western Hemisphere warm pool (WHWP). Additionally, monthly precipitation totals from the HadGEM2-ES GCM were dynamically downscaled by the Met Office Hadley Centre HadRM3P Regional Climate Model [42] (also referred to as PRECIS) were run at a 0.22° × 0.22° (~25 km) spatial resolution over the CA domain for the period 1950 to 2100. The Representative Concentration Pathway (RCP) emission scenarios, as proposed in the Fifth Assessment Report (AR5) by Intergovernmental Panel on Climate change (IPCC) were adopted in this study. RCPs refer to the possible range of radiative forcing values in the year 2100 relative to pre-industrial values (+2.6, +4.5, +6.0 and +8.5 W/m², respectively), depending on how much greenhouse gases are eventually emitted in the years to come.

Since emissions trends can vary significantly from gradual development to rapid industrialized growth, RCPs 2.6, 4.5 and 8.5 were considered in this study to provide a wider picture of possible climate change impacts over the coming century. The RCP2.6 is a low emission scenario representing an equivalent concentration of CO₂ of 440 ppm by mid-century. The RCP4.5 is an intermediate stabilization emission scenario with an equivalent concentration of CO2 ranging from 580 to 720 ppm in 2100, whereas the RCP8.5 is a rising pathway scenario leading to an equivalent concentration of CO₂ larger than 1000 ppm in 2100. It is worth noticing that CanESM2-RCA4 was only run under RCP8.5, whereas MIROC5-RCA4 excluded RCP4.5 (Table 1). For instance, the multimodel ensemble-mean (EM) was calculated based on the total number of GCM-RCMs pairs available for specific periods. Because of the various spatial resolutions adopted by different GCM-RCM pairs, all model outputs, including the local climatology were resampled to a common 25 × 25 km grid using the bilinear interpolation method to match the spatial resolution of the gridded PRECIS and GPCC outputs. All data processing and spatial interpolation, including data-converting, masking, resampling, dating, reprojecting and raster-stacking were executed using the R programming language (v3.5.2) [86] by combining functionalities of the gstat (v2.0.3), sp (v1.3.1), raster (v3.0.7), ncdf4 (v1.16.1) rgdal (v1.4.6) and rgeos packages. Reprojection from WGS84 followed the official Costa-Rica Transverse-Mercator (CRTM05) projected-coordinate-system (Figure 3).

2.3. Bias Correction Methods

Six bias correction (BC) methods covering the major types of existing methods were compared in this study to adjust the biases of precipitation GCM-RCM outputs over Costa Rica (Table 1). The BC methods consist of delta-method (DT) linear scaling (LS), power transformation of precipitation (PTR), empirical quantile mapping (EQM), gamma quantile mapping (GQM) and gamma-pareto quantile mapping (GPQM). All of them are conducted on a monthly basis for each calendar month. A detailed explanation of these methods can be found in literature [2,6,22,29,32,43,69,87]. A brief description is provided below.

2.3.1. Delta-Change (DT)

The delta-change method [88,89,90] generates climate scenarios by adding the future change signal (anomalies) from GCM-RCM simulations for a perturbation of the observational datasets rather than to use the GCM-RCM simulations of future conditions directly according to:

$$P_{contr}^{BC}\left(t\right)=P_{obs}\left(t\right)$$

(1)

$$P_{frc}^{BC}\left(t\right)=P_{obs}\left(t\right)·\left[\frac{{\mu }_m{P}_{frc}\left(t\right)}{{\mu }_m{P}_{contr}\left(t\right)}\right],$$

(2)

where P is precipitation, contr is GCM-RCM simulated time series during the control period, obs is the observational time series during the control period, frc is the future forecast time series to be corrected, BC is the final bias-corrected time series, t is the time step and µ_m is the long-term monthly mean.

For the future forecast scenarios, the GCM-RCM-simulated anomalies between control and scenario runs are superimposed upon the observational time series, which is commonly done on a monthly basis using a multiplicative correction precipitation. In the delta-change method, observational time series are used as baseline which makes it a robust method capable of correcting the mean values. Nevertheless, standard deviation (variance), wet frequencies and intensities are not corrected. By definition, the DT perfectly adjusts GCM-RCM simulated time series during the calibration but potential future changes in climate dynamics and variability are not accounted for due to its simple transfer nature.

2.3.2. Linear Scaling (LS)

The linear scaling (LS) method aims to perfectly match the long-term monthly mean of corrected values with those of observed values [59]. Precipitation is typically adjusted with a multiplier factor. By definition, corrected GCM-RCM simulations will perfectly agree in their monthly mean values with the observations. The LS method operates with monthly correction values based on the differences between observed control and uncorrected data:

$$P_{contr}^{BC}\left(t\right)=P_{contr}\left(t\right)·\left[\frac{{\mu }_m{P}_{obs}\left(t\right)}{{\mu }_m{P}_{contr}\left(t\right)}\right]$$

(3)

$$P_{frc}^{BC}\left(t\right)=P_{frc}\left(t\right)·\left[\frac{{\mu }_m{P}_{obs}\left(t\right)}{{\mu }_m{P}_{contr}\left(t\right)}\right].$$

(4)

2.3.3. Power Transformation of Precipitation (PTR)

Since transformation methods like DT and LS account for a bias in the mean, they do not correct standard deviation (variance), wet frequencies or intensities. Consequently, a non-linear correction can be used to adjust the variance statistics of precipitation time series. The power transformation of precipitation method (PTR) applies an exponential correction factor in the form _aP^b [30,91]. The Parameter b is estimated with a distribution-free approach on a monthly basis using a 90 time steps window centered on the interval. Initially, b is identified by matching the coefficient of variation (CV) of the corrected GCM-RCM precipitation (P^b) with the CV of observed precipitation (P_obs) for each month (m) using a root-finding algorithm based on Brent’s method [92] finding b_m such that:

$$f\left(b_m\right)=0=C{V}_m\left(P_{obs}\left(t\right)\right)-C{V}_m\left(P_{contr}^{bm}\left(t\right)\right)$$

(5)

$$P_{contr}^{BC,1}\left(t\right)=P_{contr}^{bm}\left(t\right)$$

(6)

$$P_{frc}^{BC,1}\left(t\right)=P_{frc}^{bm}\left(t\right).$$

(7)

Afterwards, the long-term monthly mean of observed precipitation is matched with the monthly mean of the intermediary series precipitation BC,1 using the standard linear scaling parameter abovementioned:

$$P_{contr}^{BC}\left(t\right)=P_{contr}^{BC,1}\left(t\right)·\left[\frac{{\mu }_m{P}_{obs}\left(t\right)}{{\mu }_m{P}_{contr}^{BC,1}\left(t\right)}\right]$$

(8)

$$P_{frc}^{BC}\left(t\right)=P_{frc}^{BC,1}\left(t\right)·\left[\frac{{\mu }_m{P}_{obs}\left(t\right)}{{\mu }_m{P}_{contr}^{BC,1}\left(t\right)}\right].$$

(9)

In PTR, a precipitation threshold can be introduced a priori to minimize drizzle events (very low but non-zero precipitation).

2.3.4. Empirical Quantile Mapping (EQM)

Quantile mapping (QM) attempts to correct the distribution of the GCM-RCM simulated data so that it matches the distribution of the observational dataset [93,94], which is done by applying a transfer function (TF). The TF maps the GCM-RCM simulated data to the observations based statistical moments of their respective cumulative distribution functions (CDF):

$$P_{contr}^{BC}\left(t\right)=CD{F}_{obs,m}^{-1}\left(CD{F}_{conr,m}\left(P_{contr}\left(t\right)\right)\right)$$

(10)

$$P_{frc}^{BC}\left(t\right)=CD{F}_{obs,m}^{-1}\left(CD{F}_{frc,m}\left(P_{frc}\left(t\right)\right)\right).$$

(11)

QM non-linearly corrects the mean, standard deviation (variance), quantiles, wet frequencies and intensities preserving the extremes. Therefore, the variability of corrected data is more consistent with the original GCM-RCM data. Similar to PTR, a precipitation threshold can be introduced to avoid substantial distortion of the distribution caused by drizzle events. However, QM lies on the assumption that both the observed and GCM-RCM datasets follow the same proposed distribution, which may introduce potential new biases. If empirical non-parametric CDFs are used without any assumption on precipitation distribution, the method is referred to as empirical quantile mapping (EQM) [55,64,95].

2.3.5. Gamma Quantile Mapping (GQM)

Quantile mapping (QM) can also be used to correct the biases of GCM-RCM projections by replacing the empirical CDFs with a parametric distribution. Similar to EQM, it is assumed that both the GCM-RCM simulated and observed precipitation follow a specific frequency distribution. The gamma distribution is commonly used for representing the CDF of precipitation and the method is referred as to gamma quantile mapping (GQM) [50]. The gamma distribution depends on a shape parameter α and a scale parameter β:

$$f\gamma \left(x|\alpha ,\beta \right)=x^{\alpha -1}·\frac{1}{{\beta }^{\alpha }·\Gamma \left(\alpha \right)}·e^{\frac{-x}{\beta }};x\ge 0;\alpha ,\beta \succ 0.$$

(12)

The shape parameter α controls the profile of the distribution, whereas the scale parameter β determines the dispersion of the Gamma distribution, which has been proven to be effective in precipitation [64,93,96].

2.3.6. Gamma-Pareto Quantile Mapping (GPQM)

The combination of a gamma distribution and a Generalized Pareto distribution in quantile mapping is referred to as gamma-Pareto quantile mapping (GPQM) [36]. By this combination, a better correction of higher percentiles (extremes) is expected. It is assumed that the lower 95% of the distribution follow the gamma distribution. The upper 5% are assumed to follow the General Pareto distribution. GPQM has successfully been applied for the correction of climate change scenarios [63,97].

2.4. Bias Correction Implementation and Performance Assessment

The evaluation of the various bias correction (BC) methods was based on a two-sample cross-validation approach [21,29], where observed and GCM-RCMs precipitation outputs were split into a 30-year calibration period (1951–1980) and a 15-years validation period (1981–1995), defined mainly on the basis of the observational database temporal extension and data quality and consistency (Figure 3). The R programming language (v3.5.2) [86] package hyfo (v1.4.0) was used to execute the empirical quantile mapping (EQM), gamma quantile mapping (GQM) and linear scaling (LS) methods, whereas the downscaleR (v3.0.5), loadeR (v1.4.12) and transformeR (v1.4.7) packages were used to execute the delta-method (DT), gamma-pareto quantile mapping (GPQM) and power transformation of precipitation (PTR) methods independently for each grid cell within the various climatic regions (Figure 1). Bias correction parameters were generated and applied separately for each month of the year for each of the GCM-RCMs [22].

In the execution of the various BC methods, specialized R functions processed observed and GCM-RCMs simulated gridded time series data of the same duration for the calibration period (1951–1980) to adjust the biases of the validation period (1981–1995) time series. Afterwards, during the performance evaluation stage, the root-mean square error (RMSE), the mean absolute error (MAE), the percentage bias (PBIAS) and the correlation coefficient (R) objective functions were applied to quantify the bias between each individual BC method output and the observations for the validation period (1981–1995) according to:

$$RMSE=\sqrt{\frac{1}{n}{{\displaystyle \sum _{i=1}^n\left(M_i-O_i\right)}}^2}$$

(13)

$$MAE=\frac{1}{n}{\displaystyle \sum _{i=1}^n\left|M_i-O_i\right|}$$

(14)

$$PBIAS=100\left(\frac{{\displaystyle \sum _{i=1}^n\left(M_i-O_i\right)}}{{\displaystyle \sum _{i=1}^n\left(O_i\right)}}\right)$$

(15)

$$R=\frac{\mathrm{cov}\left(M,O\right)}{\sqrt{\mathrm{var}\left(M\right)}\mathrm{var}\left(O\right)},$$

(16)

where: M_i and O_i are the modelled and observed values respectively.

The RMSE ranges from 0 to +∞ and it is used for checking the estimation accuracy between observed and modelled values. A RMSE value close to 0 indicates a higher accuracy in estimation. The MAE is an absolute measure of bias that varies between 0 to +∞. A MAE value close to 0 indicates an unbiased prediction. Both the root mean square error (RMSE) and the mean absolute error (MAE) have extensively been employed in climate change and hydrological model evaluation studies [2,6,29,98,99,100,101,102].

The PBIAS on the other hand, calculates the relative volume difference between modelled and observed volume. A negative value indicates under-prediction whereas a positive value indicates over-prediction. Additionally, the standard deviation (SD) and the correlation coefficient (R) were also used to check the performance of each BC method mainly through the use of Taylor diagrams and heatmaps. Complementarily, quantile- quantile plots (Q-Q plots) were used to compare the uncorrected and bias-corrected CDFs against observations for the multimodel ensemble-mean (EM) during the validation period. Metrics from uncorrected GCM-RCMs and GPCC outputs were also included in the cross-validation process for comparison purposes. Based on the validation (1981–1995) period performance, selected BC methods were applied to correct the biases of precipitation projections from each individual GCM-RCM pair (Table 1) for a near-future (2011–2040), mid-future (2041–2070) and far-future (2071–2100) period under IPCC Representative Concentration Pathways (RCPs) 2.6, 4.5 and 8.5 (Figure 3). As suggested by the World Meteorological Organization [103], the standard 30-year control period (1961–1990) was used for the observed and GCM-RCM simulated gridded time series data, whereas the forecast time series are associated with future periods. After the bias corrected future time series were obtained using the abovementioned R functions, the percentage-change (PC) between the mean values of the observed control period (1961–1990) and the bias corrected ensemble-mean future projections were compared for each RCP scenario and future period. After standardization from the various netCDFs internal classes, future gridded time series products were exported to GeoTIFF format.

3. Results and Discussion

3.1. Historical Simulation

The spatially-averaged annual precipitation cycle over Costa Rica (Figure 4) shows that considerable biases exist between uncorrected GCM-RCM outputs and observations, which ultimately advise against the direct application of GCM-RCM precipitation projections in local climate change studies over the Costa Rica territory. The GCM-RCM pair ability to reproduce the shape of the observed bimodal annual precipitation cycle largely depends on seasonality, climatic region and spatial resolution. During DJF, when the lower surrounding sea surface temperatures (SSTs) and weaker trade winds within the intertropical convergence zone (ITCZ) mark the beginning of the dry season for most of the country; the majority of GCM-RCM pairs exhibit a systematic dry bias except for the Central-Valley region, where most GCM-RCM pairs overestimate precipitation. The same dry bias is observed during the transitional period from dry to wet conditions (MAM), when sea breeze effects and monsoonal circulations are triggered within the ITCZ although less evident for the Caribbean region.

During JJA, when all climatic regions experience a summer precipitation peak as northeast trade winds get stronger and SSTs also increase through the ITCZ, most GCM-RCM pairs overestimate precipitation in all regions but the North, where the majority of simulations substantially underestimate precipitation. Large wet biases can also be observed during the period of extreme precipitation (SON), when Atlantic tropical cyclone activity is most frequent and hurricanes and tropical storms drive onshore winds on the west side of the cordilleras promoting heavy precipitation along the Pacific and Atlantic coasts. Once again, the opposite applies for the North region by late SON, as most models exhibit a tendency to underestimate precipitation. Notwithstanding, CanESM2, MIROC5 and in many instances HadGEM2-ES, exhibit a systematic dry bias throughout the year, with CanESM2 reproducing only one of the two precipitation peaks commonly seen during JJA-SON. At the same time, MIROC5 exhibits an unexpected minimum during September in most of the country. The mid-summer drought of July-August however, is only captured by PRECIS, HadGEM2-ES and MIROC5 but exaggerated in most cases.

In general terms, GCM-RCM pairs downscaled by RCA4 exhibit higher biases when compared against PRECIS (HadGEM2-ES/HadRM3P) regardless of climatic region. Most RCA4 models face difficulties in representing the observed wetter conditions experienced by climatic regions along the Pacific coast, together with specific climate processes operating in those region, including complex land and coastal configurations between the northeastern and southwestern parts of the country, which ultimately renders a very location-dependence performance. In addition, the fact that GCMs dynamically downscaled by RCA4 show such extreme variations suggests that the performance of individual model pairs largely depends on GCMs simulations themselves rather than the conceptualization of RCA4, indicating that local processes and relevant parameterizations are important in determining RCMs responses to boundary forcing. As pointed out by Kotlarski et al. [57], biases related to parameterization schemes and unresolved processes thus remain, in addition to systematic biases related to the large-scale circulation patterns as part of GCMs outputs that later propagate into RCMs. The fact that the two RCMs considered (HadRM3P and RCA4) driven by the same GCM (HadGEM2-ES) exhibit considerably different biases suggests that RCMs spatial resolution and internal parameterization are equally important.

The lower spatial resolution (50 km), convection schemes and surface flux parameterization of the RCA4 could also be associated with such a tendency, particularly over Central America and especially Costa Rica, which narrowest part is barely 125 km wide (Figure 1) between the Pacific and Atlantic oceans. Unwanted marine influence from both oceans could collaterally be incorporated into land cells during the bilinear resampling process to 25 km due to spatial averaging.

A similar tendency was found by Diro et al. [40] using the RegCM4 Regional Climate Model to conduct an analysis over the Central America CORDEX domain at 50 km × 50 km spatial resolution to evaluate the sensitivity of rainfall to land and sea surface schemes. Their results identified a marked underestimation of precipitation mainly related to an erroneous numerical misplacement of the location of the Caribbean Low Level Jet (CLLJ) by the GCM-RCM.

In their study of the CORDEX Central America domain using the RegCM4 Regional Climate Model, Fuentes-Franco et al. [104] concluded that the representation of topography over the narrow Central America isthmus critically depends on GCM-RCM model spatial resolution. The authors also highlighted the relevance of resolving coastal topography for a better simulation of tropical cyclone activity in the region. The complex physiographical features of Costa Rica, which include complicated mountainous terrain in a relatively small area, highlight the need for higher resolution GCM-RCM models.

Olsson et al. [105] also encountered limitations on the 50 km × 50 km CORDEX spatial resolution as it seems too course for small scale hydrological modelling. Conversely, the higher spatial resolution (25 km) and scheme parameterization of PRECIS, which more clearly describes the boundary between land and ocean cells, seems to capture more of the spatiotemporal variability of precipitation over most climatic regions of Costa Rica, especially during the rainy season. This indicates that as individual GCM-RCM pairs; precipitation variability is better reproduced by PRECIS, as higher resolution outcomes add improvement over the coarser RCA4 50 km simulations in most cases. Furthermore, the year-round persistent and complex precipitation regime of the North and Caribbean regions, which are associated with cold air outbreaks and frontal systems, is better described by PRECIS than any RCA-RCM pair.

Similarly, Karmalkar et al. [25] analyzed two sets of experiments at 25 km × 25 km spatial resolution over Central America using PRECIS and found a reasonable representation of the key regional climate features for recent climate conditions. Nevertheless, the authors underline that a better representation of the Mid-Summer Drought (MSD), the Caribbean Low Level Jet (CLLJ), cold surges and tropical cyclones dynamics are needed to improved reliability on GCM-RCM projections. It is worth noticing that the uncorrected model ensemble-mean generally performed better than individual RCA4 pairs and in many instances even better than PRECIS. Such is the case of the North-Pacific, Central-Pacific and South-Pacific regions, as the ensemble-mean simulations get closer to observations across most seasons.

For the remaining climatic regions however, the uncorrected model ensemble-mean seems to overestimate precipitation during the rainy season and underestimates precipitation during the dry season. This is possibly related to the severe dry biases incorporated by CanESM2 and MIROC5 within the spatially-averaged mean, especially over the Central-Valley and North regions. In spite of clear inabilities of individual GCM-RCM pairs to properly reproduce the magnitude of historical precipitation, all models except CanESM2 and MIROC5 generally simulate well the annual cycle of precipitation over most climatic regions of Costa Rica. For instance, the outcomes of these models could be used for climate change impact studies after the application of bias correction methods. In addition, observations from the GPCC dataset are generally in good agreement with the high resolution local precipitation climatology except for the Central-Valley region, where GPCC slightly overestimates precipitation during most of the rainy season. However, on a comparative basis, the overall spatial and temporal coherence between the local climatology and the GPCC products is considered satisfactory, which gives sufficient assurance that the local climatology can be used for bias correction purposes over the Costa Rica territory.

3.2. Evaluation of Bias Correction Methods

Taylor diagrams (Figure 5) summarize the results of the various BC methods applied to all GCM-RCM pairs during the validation period (1981–1995) for each climatic region of Costa Rica in terms of root mean square error (RMSE), correlation coefficient (R) and standard deviation (SD). Results show that most BC methods considerably improve GCM-RCM simulations as RMSE values clearly decrease, SD values get closer to observations and R values get higher after the application of BC methods. Nonetheless, in most cases the extent of improvements is fairly comparable.

Uncorrected GCM-RCM R values are generally lower than their bias corrected counterparts, suggesting than BC values better represent the temporal and spatial patterns of precipitation over Costa Rica. Moreover, SD values of uncorrected GCM-RCM pairs vary significantly across climatic regions, showing marked spatial and temporal heterogeneity. The uncorrected ensemble-mean usually exhibits better performance compared to individual uncorrected GCM-RCM pairs. The ensemble-mean also benefit from bias correction application but there is a tendency to over-adjust SD values by most BC methods in the South-Pacific region. Reasonably higher R values (between 0.70 and 0.85) and lower RMSE values (100 to 125 mm/month) are achieved after bias corrected for most climatic regions except for the Caribbean, where R values in fact decrease from a range of 0.60–0.70 to values below 0.50, even though RMSE values are substantially reduced. This indicates that uncorrected GCM-RCM pairs encounter limitations in properly capturing the spatiotemporal variability of precipitation over the Caribbean climatic region and BC methods are for instance even less effective. In addition, uncorrected RMSE and MAE values are also usually higher for RCA4 models, particularly along the Caribbean coast where the year-round persistent precipitation pattern, associated with cold air outbreaks and frontal systems advect precipitation is clearly not properly described by RCA4 models. Since BC methods results vary substantially across temporal and scales, seasonal heatmaps of performance metrics are presented in terms of root mean square error (RMSE), mean absolute deviation (MAE), correlation coefficient (R) and percentage bias (PBIAS) for the entire dataset (Figure 6, Figure 7, Figure 8 and Figure 9), along with Q-Q plots of uncorrected and bias corrected ensemble-mean data against observations (Figure 10).

3.2.1. DJF

At the beginning of the dry season (DJF), the uncorrected model ensemble-mean (EM) generally outperforms individual GCM-RCM pairs, exhibiting lower RMSE, MAE and PBIAS values and higher R values regardless of climatic region or spatial resolution, which suggests that the ensemble-mean is more effective in capturing the seasonal pattern of precipitation over Costa Rica during DJF. Improvements for the ensemble-mean and individual GCM-RCM pairs are clearly achieved for the North-Pacific, Central-Pacific, Central-Valley and North regions after the application of bias correction methods, with EQM and DT showing marginally lower RMSE, MAE and PBIAS values than the remaining methods. For the North-Pacific, the ensemble-mean PBIAS was reduced from -32% to 3% using EQM, 73% to −5% for the Central-Valley using DT, −62% to 14% for the North region using DT and 16% to 4% for the Central-Pacific using DT as well. Nonetheless, PTR, GQM and GPQM provide similar performances as well. Reductions in terms of RMSE range from 8% to 78% after bias correction; in terms of MAE such reductions range from 24% to 66% depending on the GCM-RCM pair, climatic region and BC method.

Similar to uncorrected values, the bias corrected model ensemble-mean is generally better than individual bias corrected GCM-RCM pairs, only slightly exceeded by PRECIS in some cases. Q-Q plots for the ensemble-mean also support the application of BC methods, since bias corrected quantiles are closer to observed values than their uncorrected counterparts, regardless of the BC method. It can also be noted that uncorrected quantiles show larger variations across climatic regions, especially for the North and Central-Valley regions.

Notwithstanding, for the South-Pacific and Caribbean regions improvements are mainly achieved on RCA4 models in terms of RMSE MAE and R when bias corrected by EQM and DT respectively but biases from PRECIS and the ensemble-mean were generally amplified after bias correction. PBIAS values in both regions also increase when compared to uncorrected values for PRECIS and the ensemble-mean. It can be seen for example that PRECIS PBIAS increases from −1% when uncorrected to 38% after bias corrected by DT in the Caribbean region and from −19% to 49% for the South-Pacific using EQM. A similarly tendency is shown by the ensemble-mean in the Caribbean region, probably driven by PRECIS influence.

Concerning R, the application of most BC methods ultimately decreases correlation for PRECIS and the ensemble-mean. In addition to that, the ensemble-mean bias corrected PBIAS generally shows a systematic wet bias for these two regions, which can be observed in the Q-Q plots.

It is clear that the added value of bias correction on these two regions is considerably limited and even detrimental for PRECIS during DJF, which indicates that all bias correction methods fail to correct the spatial and temporal pattern of precipitation. In this case in particular, bias correction may introduce unnecessary variances compared to those already delivered by HadRM3P itself during the dynamical downscaling process of HadGEM2, which evidently discourages the use of BC methods and promotes the use of uncorrected PRECIS data instead. Consequently, the uncorrected PRECIS outcomes are preferred over their bias corrected counterparts over the South-Pacific and Caribbean climatic regions during DJF. Contrastingly, during DJF the GPCC climatology generally benefited from bias correction for all climatic regions except for the South-Pacific regions, an issue that has been previously identified and analyzed Mendez et al. [82]. In this case, EQM scores slightly higher in terms of RMSE, MAE and R but PTR yields equally satisfactory performances. In terms of PBIAS however, EQM and GQM score comparably well, with corrected PBIAS below 5 mm/month for all climatic regions except the Caribbean, which exhibits a corrected BIAS value in the order of 12 mm/month. Positive PBIAS values indicate over-estimation, whereas negative values indicate under-estimation. Overall, a significant value is added to the GPCC climatology after bias correction for all regions except the South-Pacific.

3.2.2. MAM

The transitional pattern from dry to wet conditions experienced during MAM, when sea breeze effects and monsoonal circulations are triggered within the intertropical convergence zone (ITCZ) promoting precipitation along the Atlantic and Pacific coasts, results in lower variations of performance and ranking among the various BC methods within each climatic region when compared to DJF. Similar to DJF, uncorrected ensemble-mean RMSE, MAE values are generally lower than individual uncorrected GCM-RCM pairs. After bias correction, the ensemble-mean performs as good as or even better than individual RCA4 models or PRECIS in terms of RMSE, MAE and PBIAS, with EQM and DT slightly outperforming the remaining BC methods. In some instances, the bias corrected ensemble-mean ranks moderately lower than its uncorrected counterpart but many individual GCM-RCM pairs perform better than the ensemble-mean as a whole. On the other hand, Q-Q plots prove the added value of bias correction over the ensemble-mean, as the distance between model and observation is reduced and the adjusted distributions better match the observations, especially for the North-Pacific, South-Pacific and North regions. RCA4 models perform as good as or even better than PRECIS before or after bias correction, implying than the higher spatial resolution of PRECIS is of little advantage during MAM. In terms of PBIAS, EQM usually yields lower values when compared to the remaining BC methods for most climatic regions; for the North and Caribbean regions however, DT scores better, passing from −47% to 1% and from 1% to 7% respectively.

Similar to DFJ, the application of most BC methods results in lower R values for the ensemble-mean but generally improved individual GCM-RCM pairs including PRECIS. Concerning GPCC, the Central-Valley and Caribbean regions do exhibit improvements in terms of RMSE, MAE and PBIAS after bias correction, mainly EQM and PTR. For the Central-Valley, GPCC PBIAS is reduced from 21% to 3% by EQM, in the Caribbean region, from 14% to −2% by PTR. Nevertheless, for all other regions bias correction resulted in an adverse effect over most metrics considered.

3.2.3. JJA

During JJA, all climatic regions in Costa Rica experience a precipitation peak, as northeast trade winds get stronger and sea surface temperature (SST) also increase through the ITCZ, promoting heavy convective activity and storm development. As a result, uncorrected and bias corrected RMSE and MAE values get remarkably higher than those shown by either the ensemble-mean or individual GCM-RCM pairs during DJF or MAM. Percentage reductions on these two metrics also increase notably, in particular for the Central-Pacific, Central-Valley and Caribbean regions, with values well above 50% in many cases. Significant reductions and less variability in terms of PBIAS are accomplished by most BC methods for all climatic regions, which indicate that BC methods are substantially more beneficial for GCM-RCM pairs during the rainy season than during the dry season.

For most of the country, EQM, PTR and GPQM yield comparable results, with EQM marginally outperforming the remaining methods. The fact that EQM moderately ranks higher than the remaining BC in so many instances, is most likely related to the fact that the nonparametric transformation technique implemented by this method adjusts the quantiles variabilities with those of the observed data; matching the quantiles of regional empirical cumulative distribution functions (CDFs) of GCM-RCM pairs with those of the local climatology. The superior performance of EQM could also be attributable to the adopted monthly approach used in the BC methods application, which takes different precipitation levels into account on individual basis. Various bias correction evaluation studies have also concluded that EQM is typically superior when compared to other BC methods [21,22,32,50,52,53,55,87,106,107,108]. Contrarily, the distribution-based parametric methods (GQM and GPQM) show slightly larger bias when compared to empirical distributions (EQM), as these methods are based on the assumption that both GCM-RCM and observed data approximate the corresponding theoretical distributions functions [29,35], which could not always be the case for the temporal and spatial precipitation distribution within most climatic regions, especially during the dry season.

Gudmundsson et al. [87] argues that non-parametric methods usually produce a better performance as compared to the parametric methods in the context of QM bias adjustment. For the Caribbean region however, DT is superior to the remaining methods in terms of RMSE and MAE. Once more, after bias correction, the ensemble-mean mean exhibits RMSE, MAE and PBIAS values as low as or lower than individual GCM-RCM pairs. Q-Q plots reveal large variations for the uncorrected ensemble-mean across all climatic regions, particularly for the Caribbean and South Pacific regions, which also exhibit the largest RMSE and MAE uncorrected or bias corrected values. Substantial improvements in terms of CDFs can be observed for all climatic regions after bias correction, nevertheless there are still some discrepancies in the spatial-temporal distribution of biases between uncorrected and bias corrected quantiles relative to the observations for the South-Pacific and North regions, since bias corrected model curves line far from the theoretical dashed line, indicating that GCM-RCM data have different distributions when compared to observed data.

Contrary to DJF and MAM, bias corrected ensemble-mean R values effectively increased for most BC methods the Pacific regions and Central-Valley but R values are substantially lower (below 30%) when compared to DFJ or JJA. For the North and Caribbean region however, bias corrected PBIAS values are generally higher than those of the remaining regions, with a noticeable dry-bias which ranges from −5% to −25%. Furthermore, R exhibits predominantly low (lower than 0.3) close to zero or even negative values, suggesting a reversal of the observed annual precipitation cycle regardless of BC method or GCM-RCM pairs.

One of the factors contributing to such a tendency is the locally intense convective activity of JJA precipitation across most of the Central-American Caribbean coast [27,109], which in this case most GCM-RCM models fail to capture and for instance, BC methods render less effective. In general terms, the diminished R and PBIAS performances for these two regions are most likely related to GCM-RCMs difficulties at representing predominantly convective precipitation during JJA; a situation that has largely been identified and investigated [21,58,70,80]. On the other hand, GCM-RCM pairs can also suffer from considerable systematic biases of the precipitation field in orographic terrain such as Costa Rica due to the windward and leeward effects imposed by the northwest-southeast trending cordilleras, producing overestimation on the windward side and underestimation on the lee side [57,59].

The full amount of these complex factors could also explain why DT, being simpler than the remaining BC methods, normally exhibits better performances than more complex methods such as EQM or GPQM for the North and Caribbean regions, which could also be limited to monthly temporal resolution data as suggested by Yang et al. [47], Shrestha et al. [110] and Hanel et al. [111]. Nevertheless, the diminished GCM-RCM model performances in terms of R and PBIAS does not necessarily indicate severe model biases but rather shows that GCM-RCM pairs cannot properly reproduce JJA monthly deviations within the annual precipitation distribution cycle. In relation to the GPCC dataset, only the Central-Valley region shows clear improvements in terms of RMSE, MAE and PBIAS when corrected by EQM, PTR or GPQM, all of which yield comparable results. For the remaining climatic regions however, GPCC does not benefit from bias correction irrespectively of the BC method applied, this indicates that bias correction is more effective for the GPCC during the dry season and little added value is gained during the rainy season.

3.2.4. SON

During SON uncorrected RMSE and MAE values for most RCA4 models are generally higher than those experienced during the first precipitation peak of JJA. As previously mentioned, SON represents the rainiest period for all climatic regions, evidencing that biases increase as precipitation also increases in the case of RCA4. For PRECIS however, uncorrected RMSE and MAE values during SON are not only lower than any RCA4 model but also lower compared to JJA their counterpart, suggesting that PRECIS is more effective in describing cyclonic-related precipitation than convective precipitation. After BC nonetheless, comparable performances between RCA4 pairs and PRECIS can be observed, with EQM and PTR slightly outperforming the remaining BC methods for most climatic regions except for the North and Caribbean regions, where once more DT seems to better describe the precipitation dynamics. Regardless of RCA4 or PRECIS, RMSE and MAE percentage reductions values accomplished by these BC methods (well above 50% in some cases) are usually higher than those obtained during JJA, indicating once again that the wetter the season the more effective bias correction is. Much like JJA, RCA4 models benefit more from bias correction than PRECIS. For most climatic regions, the bias corrected ensemble-mean exhibit RMSE and MAE values as low as or even lower than individual GCM-RCM pairs. Marginally higher RMSE and MAE values can be observed for PRECIS and the ensemble-mean after bias correction when compared to uncorrected values for all Pacific climatic regions. Nevertheless, such RMSE and MAE performances are still as good as individual bias corrected GCM-RCM pairs.

Contrary to JJA, bias corrected PBIAS values during SON exhibit a predominantly wet-shift (between 5% and 20%) for all climatic regions regardless of GCM-RCM pair or BC method, which could also be related to the primarily cyclonic nature of precipitation during SON. Still, such wet biases are clearly lower for EQM, PTR and DT. In terms of R, significant added value can be observed for North-Pacific, Central-Pacific, South-Pacific and Central-Valley regions, as increasingly positive values are obtained after bias correction irrespectively of GCM-RCMs pairs. Meanwhile, predominantly negative values still remain after the application of bias correction for the North and Caribbean regions. The ensemble-mean Q-Q plots show considerable added value as bias corrected quantile get much closer to observations than uncorrected values, especially for EQM, PTR and DT.

On the other hand, LS generally proved less effective in terms of RMSE, MAE, PBIAS and R when compared to the remaining BC methods during SON irrespectively of climatic region. Similar results were obtained by Teutschbein and Seibert [69] in their evaluation of several bias correction methods over five Swedish catchments using 15 different RCMs driven by ERA40, where LS resulted in slightly larger model error growth and deviations compared to more advanced correction methods such as QM, which makes LS the least reliable method under changed conditions. In contrast, Shrestha et al. [110] found no significant differences between the performances of LS and EQM bias corrected outputs for the Kaligandaki River Basin in Nepal. Such discrepancies however, are most likely related to substantial climatic differences between the Himalayas, Sweden and the Central-American tropics. It also must be noted that LS is essentially used to scale the GCM-RMCs model values with the multiplicative differences between the observed and model means during the validation period.

Under the present configuration, LS also uses unique scaling factors during specific months where no extreme precipitation is specifically considered, often leading to heavy precipitation being greatly underestimated [52,112]. Similar to JJA, GPCC clearly benefits from bias correction mainly for the Central-Valley region, with most BC rendering similar performances. For all other climatic regions however, GPCC gains little from bias correction and for the South-Pacific, the performance in fact decreases in terms of all metrics considered.

3.3. Future Projections

Based on the validation period (1981–1995) results, it can be stated that in the case of Costa Rica bias correction methods considerably decrease the uncertainty of monthly precipitation projections compared to uncorrected GCM-RCM outputs. It can also be stated that the ensemble-mean outperforms individual GCM-RCM pairs in most cases. Therefore, EQM and DT were selectively applied to correct the biases of precipitation projections from each GCM-RCM pair (Table 1) for three 30-year future periods namely, near-future (2011–2040), mid-future (2041–2070) and far-future (2071–2100) under IPCC Representative Concentration Pathways (RCPs) 2.6, 4.5 and 8.5 using the 30-year control period 1961–1990 (

Table 2. Selection of bias correction methods applied to individual GCM-RCM pairs for future periods over all climatic regions of Costa Rica with respect to baseline (1961–1990).

BC-Method/Season	North-Pacific	Central-Pacific	South-Pacific	Central-Valley	North	Caribbean
DJF	EQM	EQM	EQM	DELTA	DELTA	DELTA
MAM	EQM	EQM	EQM	EQM	DELTA	DELTA
JJA	EQM	EQM	EQM	EQM	EQM	DELTA
SON	EQM	EQM	EQM	EQM	EQM	DELTA

).

The abovementioned methods were selected on the basis of marginally higher performances, which should not discourage the use of the remaining BC methods such as GQM, GPQM and PTR. In the case of PRECIS however, DJF uncorrected values were used instead of bias corrected projections for the South-Pacific and Caribbean regions. After individual bias correction of future GCM-RCM pairs, the bias corrected ensemble-mean was selected to provide future projections of precipitation over Costa Rica. This approach has repeatedly been implemented in literature [5,57,64,93,113,114]. For future projections, the monthly averaged precipitation anomalies between the mean values of the observed control period (1961–1990) and the bias corrected multimodel ensemble-mean were calculated and expressed as percentage-change (PC). Results, expressed as spatially-averaged (Figure 11) and spatially-distributed (Figure 12) percentage-changes (PC) for each RCP scenario and future period contain mixed features, with both increasing and decreasing precipitation trends, with mostly well-defined common features for each climatic region and season. Similar to BC evaluation, future projections are analyzed seasonally.

3.3.1. DJF

During DJF, results from the bias-corrected future ensemble-mean anticipate a marked decreasing trend in precipitation with increasing radiative forcing and future period for most of the country under RCP4.5 and 8.5. Pronounced drier conditions are expected for those climatic regions draining towards the Pacific Ocean which also coincides with the driest season, exhibiting percentage-change (PC) values between −10% and −20% for the far-future period. Fairly comparable tendencies between RCP4.5 and 8.5 are appreciated for the North-Pacific, Central-Pacific and South-Pacific regions but a considerably higher decrease and more spatial variation is shown by RCP4.5 during the mid-future period for the Central-Valley region, reaching PC values as low as −25%. The magnitudes of changes in precipitation however, are less significant for those climatic regions draining towards the Caribbean Sea. For the North region, PC changes range between 5% and −5% regardless of radiative forcing or future period. For the Caribbean region, changes are more evident for RCP4.5 and 8.5 and far-future period, with spatially-averaged PC mean values in the order of −10% (+/− 5%). It must be noted that precipitation in these two climatic regions is year-round persistent and evidently humid during DJF when the rest of the country experiences the driest period (Figure 2). A similar decreasing trend in precipitation but of different magnitude can be observed under RCP2.6 when compared to RCP4.5 and 8.5 in all climatic regions as future period also increases. Nevertheless, for the near and mid-future periods changes in fact anticipate wetter conditions, with PC values as high as 10% at the beginning of the century and as low as −5% at the end of the century along the Pacific Coast regions. Once again, the decreases in precipitation are less significant for the North and Caribbean regions under RCP2.6.

3.3.2. MAM

During the transitional period from dry to wet conditions (MAM), decreasing precipitation patterns of higher magnitude can be observed for the North-Pacific, Central-Pacific and South-Pacific regions with increasing radiative forcing (RCP4.5 and 8.5) and future period, with spatially-averaged PC mean values ranging between −15% and −20% under RCP8.5 and far-future period. For the Central-Valley region, the projected decrease in precipitation is more extreme for RCP4.5 and mid-future period, reaching PC mean values of around −15%.

The North and Caribbean regions decreasing precipitation trends with increasing radiative forcing are less extreme when compared to those climatic regions draining towards the Pacific Ocean. RCP4.5 however, tends to anticipate slightly drier conditions than RCP8.5 for these two regions, with PC mean values between −5% and −10% for the far-future period. RCP8.5 on the other hand, anticipates slightly wetter conditions (PC mean values between 0% and 5%) for these two regions for the near and mid-future periods but dryer conditions at the end of the century. Furthermore, considerable spatial variations in precipitation changes are clear for the far-future period between drier mountainous areas and wetter coastal areas of the Caribbean region under RCP8.5, with PC mean values ranging from 5% to −10% respectively (Figure 12). As for RCP2.6, a slightly less significant decreasing trend is expected for those climatic regions draining towards the Pacific ocean as radiative forcing increases with future period, with PC mean values ranging between −3% and 3% in most cases.

On the other hand, mostly wetter conditions are expected for the North and Caribbean regions with increasing radiative force and future periods under RCP2.6, with PC values well above 10% at the end of the century. Historically, drier conditions in the Pacific coast during DJF and MAM are associated with a southernmost displacement of the ITCZ, a decrease in SSTs in the eastern Pacific and the Inter-American Seas (values lower than 28 °C) and an increase in trade winds strength, which results in substantially lower convective activity in most of the WHWP region [80].

Future climate change projections executed by Karmalkar et al. [25] using PRECIS under IPCC Scenario SRES-A2 (quasi-equivalent to RCP8.5) for a far-future period (2070–2100), also suggest large reductions in precipitation (around −15%) along the Pacific Coast of Costa Rica during the dry season (DJF and MAM), which the authors associate to a decrease in low mean sea level pressure (MSLP) in the Caribbean Sea and a more drastic southward shift of the Intertropical Convergence Zone (ITCZ) band, which ultimately alters trade winds magnitudes and directions. The authors also suggest that the general decrease in precipitation along the Pacific Coast of Central America and the subsequent wetter conditions along the Caribbean coast during the dry season is most likely a consequence of orographic uplift of modified north-easterly/easterly trade encountering the Cordilleras, resulting in orographic precipitation along the Caribbean slopes of the region.

Fuentes-Franco et al. [115] used RegCM4 forced by HadGEM-ES and MPI to investigate Inter-annual variability of precipitation over southern Mexico and Central America and its relationship to sea surface temperature for a far-future period (1970–2100) under RCP8.5. The authors found that future drier conditions during the dry season are most likely linked to a greater warming of the of the Tropical Northeastern Pacific (TNP) compared to the Tropical North Atlantic (TNA), which causes stronger wind fluxes from the TNA to the TNP through the Caribbean Low Level Jet and is identified as the main process responsible for these drier conditions.

3.3.3. JJA

At the beginning of the rainy season (JJA), when convective mechanisms dominates precipitation in most of the country, fairly comparable trends can be observed for RCP2.6 and 8.5 during the near and mid-future periods for the Central-Pacific, South-Pacific, Central-Valley and North climatic regions, with noticeably small changes with respect to the 1961–1990 baseline period (PC mean values between −3% and 3%) and RCP2.6 showing no significant deviation towards the end of the century. Nevertheless, these four regions experience an extreme decrease in precipitation at the end of the century with RCP8.5 (PC mean values around −20%), with higher spatial variation for the North region. RCP4.5 on the other hand, suggests drier conditions for these four climatic regions for the mid-future period (PC mean values around −10%) when compared to RCP2.6 and 8.5, which decreasing trend seems to be less extreme for the far-future period when compared to RCP8.5 (PC mean values slightly above −10%). This shows marked heterogeneity in the projected precipitation changes as per the emission scenarios and in different future periods during JJA, which could also be related to the fact that fewer GCM-RCM pairs were run under RCP4.5 (Table 1). The most extreme decrease in precipitation can be observed for the North-Pacific region as radiative forcing increases with future period, with PC mean values in the order of −35% for the far-future period and RCP8.5. Moreover, in this case RCP 4.5 and 8.5 changes in precipitation for near and mid-future periods are fairly comparable. Similar to the remaining Pacific regions, RCP2.6 remains relatively stable regardless of future period.

Remarkably complex trends among the various RCPs and future periods can be observed for the Caribbean region, as precipitation changes do not respond uniformly to radiative forcing and there is considerably less sensitivity to the various RCP scenarios. Despite the initial chaotic appearance of these patterns, RCP2.6 predicts mostly wetter conditions regardless of future period (PC mean values around 10%). RCP4.5 suggests mean drier conditions only for the mid-future period but also presents high spatial variations, with PC mean values as low as −15% and as high as 5%. For the near and far-future periods, RCP4.5 suggests mostly wetter conditions for the Caribbean region. RCP8.5 on the contrary, follows a very similar wet pattern when compared to RCP2.6 up to mid-future period, nevertheless for the far-future period, there are extreme spatial changes between mountainous regions and coastal areas, with PC mean values raging between −20% and 15%. Therefore, for the Caribbean region it is difficult to find any readily apparent long-term trend in the evolution of the various RCPs considered.

Generally, precipitation is projected to decrease during JJA with increasing radiative forcing and future period for most of the country except the Caribbean region, which dry or wet trend tend to significantly vary with future period and RCP scenarios.

3.3.4. SON

During SON, changes for all climatic regions except the Caribbean are fairly similar under RCP2.6, 4.5 and 8.5 regardless of future period, with slightly drier conditions for the near-future period (below −5%) characterized by a minor increase of precipitation mostly towards the end of the century (below 5% as PC mean values in most cases). In contrast, precipitation is expected to considerably increase for the coastal and mountainous areas of the Caribbean with radiative forcing and future period, with PC mean values well above 20% for RCP8.5 towards the end of the century. This trend is consistent with findings by Campbell et al. [19] in their analysis of 2071–2100 future climate change over the Caribbean region using PRECIS forced with HadAM3 under the A2 and B2 SRES scenarios, where the authors suggest an intensification of a gradient pattern in November–January in which the Costa Rica Caribbean region gets wetter (between 10% and 30%) and the Pacific shows a robust drying signal during June-October (between −25% and −30%).

Similarly Karmalkar et al. [25] suggest large reductions in precipitation (between −15% and −30%) along the Pacific Coast of Mexico and the Central American region during the rainy season, with a slight increase in precipitation in the Caribbean regions of Costa Rica and Panama also during the rainy season. They associate these changes to a southward shift of the Intertropical Convergence Zone (ITCZ) band, which is also related with an increase in the low mean sea level pressure (MSLP) in the Caribbean Sea and dominated by the North Atlantic subtropical high (NASH) pressure system variability. Rauscher et al. [116] suggest that the precipitation decrease in Central America is associated with an intensification and early onset of the Mid-Summer Drought (MSD) over the region driven by regional variations in SST over the surrounding equatorial ocean areas.

In their analysis of high spatial-resolution future climate change scenarios for Central America, comparable results over the Costa Rica territory were obtained by Imbach et al. [16], where projections from the HadGEM2-ES/ ETA-RCM under RCP4.5 and the mid-future period (2021–2050) suggest a general decrease in precipitation in all seasons but particularly during the rainy season (JJA and SON) for those climatic regions draining towards the Pacific Ocean. Additionally, they projected strengthening of the Mid-Summer Drought and rather wetter conditions for the Caribbean and North climatic regions. They also found a slight increase in precipitation across the Caribbean coast of Central America, especially Costa Rica and Nicaragua. Hidalgo et al. [14] found significant reductions in precipitation (between −5% and −10%) in Central America through their future climate change projections for the period 2050–2099 using a 30-member General Circulation Model (GCM) model-ensemble for the A1B emission scenario in a 0.5° × 0.5° grid over Central America. The authors also relate the projected dry pattern with a southward displacement of the Intertropical Convergence Zone (ITCZ), suggesting that the prevalence of severe drought may increase significantly in the future with higher emission scenarios.

4. Conclusions

In this study, monthly precipitation outputs from five General Circulation Models (GCMs) dynamically downscaled by two Regional Climate Models (RCMs) were evaluated over all climatic regions of Costa Rica for the period 1951–1995. Results show that the GCM-RCMs ability to reproduce the magnitude and shape of the observed bimodal annual precipitation cycle over Costa Rica largely depends on seasonality, climatic region, spatial resolution and model parameterization but pronounced spatial and temporal heterogeneities are evident between RCA4 models and PRECIS. The complex physiographical features of Costa Rica, which includes a complicated mountainous terrain in a relatively small geographical area highlights the need for a higher resolution GCM-RCM models such as PRECIS. In general, significant biases can be identified between uncorrected GCM-RCM outputs and observations, which highlight the need for the evaluation and application of BC methods over complex territories such as Costa Rica before forcing climatic data into specific climate change impact assessments. Consequently, six bias correction (BC) methods were applied to adjust the biases of historical monthly precipitation outputs from all GCM-RCM pairs using for the control period 1951–1980 and validated over the period 1981–1995. Results show that most BC methods considerably improve GCM-RCM model simulations even though the performance of the various methods clearly varies across spatial and temporal scales.

In many instances, EQM, DT and PTR performed marginally better than parametric quantile mapping (QM) techniques such as GQM and GPQM. LS however, proved less effective than the remaining methods irrespectively of season and climatic region, which clearly prevents its use over the Costa Rica territory. Bias corrections methods are noticeably more beneficial for GCM-RCMs during the rainy season than during the dry season and at the same time more beneficial for lower spatial-resolution models such as RCA4. In addition, both the uncorrected and bias-corrected multimodel ensemble-mean generally exhibit better performances when compared to individual GCM-RCM pairs, which suggests that multimodel ensembles using a range of RCMs driven by different GCMs provides an solid approach in the generation of climate change scenarios over Costa Rica. Since the application of BC methods considerably decreased the uncertainty of monthly precipitation projections over Costa Rica compared to uncorrected GCM-RCM outputs, EQM and DT were selectively applied to correct the biases of precipitation projections from each individual GCM-RCM pair for a near-future (2011–2040), mid-future (2041–2070) and far-future (2071–2100) period under Representative Concentration Pathways (RCPs) 2.6, 4.5 and 8.5 using the control period 1961–1990.

Overall, results from the bias-corrected future ensemble-mean anticipate a decreasing trend in precipitation from near to far-future periods during the dry season (DJF and MAM) for most of the country under RCP4.5 and 8.5, with pronounced drier conditions for those climatic regions draining towards the Pacific Ocean. In contrast, mostly wetter conditions are expected during the dry season under RCP2.6, particularly for the Caribbean and North regions. The greatest decrease in precipitation is projected at the beginning of the rainy season (JJA) for the far-future period under RCP8.5 for all climatic regions except the Caribbean, where mostly wetter conditions are anticipated under RCP8.5, 4.5 and 2.6 for the near and mid-future periods with large spatial variations under RCP8.5 between coastal and mountainous areas.

Regardless of future period, slight increases in precipitation with higher radiative forcing are expected during SON for most of the country excluding the Caribbean region, where precipitation is expected to substantially increase with increasing radiative forcing and future period. The results of this study, in terms of future precipitation changes are mostly in agreement with previous climate change studies over the Central America region, which most authors relate to an extreme southward shift of the Intertropical Convergence Zone (ITCZ) band and an increase in the low mean sea level pressure (MSLP) in the Caribbean Sea. This gives sufficient assurance that the future climate change precipitation scenarios produced by this study can be used in local climate change studies over the Costa Rica territory. Future variability of precipitation due to climate change will most likely impose water management challenges in Costa Rica under global warming in terms of both magnitude and spatial distribution. For instance, uncertainty associated with such projections should permeate across institutional and political structures in such a way that potential implications can be included on national climate change mitigation and adaptation plans. In the future, the authors will focus on the dynamical downscaling of a larger set of GCMs using locally-parameterized higher-resolution RCMs that better represent key aspect of the Central America climate such as the Mid-Summer Drought (MSD), the Caribbean Low Level Jet (CLLJ), northern cold surges and tropical cyclones as to provide more robust statements and higher confidence associated with future climate change precipitation projections over Costa Rica.