This work shows that local-scale climate projections obtained by means of statistical downscaling are sensitive to the choice of reanalysis used for calibration. To this aim, a generalized linear model (GLM) approach is applied to downscale daily precipitation in the Philippines. First, the GLMs are trained and tested separately with two distinct reanalyses (ERA-Interim and JRA-25) using a cross-validation scheme over the period 1981–2000. When the observed and downscaled time series are compared, the attained performance is found to be sensitive to the reanalysis considered if climate change signal–bearing variables (temperature and/or specific humidity) are included in the predictor field. Moreover, performance differences are shown to be in correspondence with the disagreement found between the raw predictors from the two reanalyses. Second, the regression coefficients calibrated either with ERA-Interim or JRA-25 are subsequently applied to the output of a global climate model (MPI-ECHAM5) in order to assess the sensitivity of local-scale climate change projections (up to 2100) to reanalysis choice. In this case, the differences detected in present climate conditions are considerably amplified, leading to “delta-change” estimates differing by up to 35% (on average for the entire country) depending on the reanalysis used for calibration. Therefore, reanalysis choice is an important contributor to the uncertainty of local-scale climate change projections and, consequently, should be treated with as much care as other better-known sources of uncertainty (e.g., the choice of the GCM and/or downscaling method). Implications of the results for the entire tropics, as well as for the model output statistics downscaling approach are also briefly discussed.
Statistical downscaling (SD) techniques are nowadays routinely applied to translate coarse-resolution output from global climate models (GCMs) to local-scale climate information required by impact and adaptation studies. These techniques, however, have been developed and applied almost exclusively for extratropical regions (Hewitson and Crane 1996; Wilby and Wigley 1997; Trigo and Palutikof 2001; Hanssen-Bauer et al. 2005; Fowler et al. 2007; Maraun et al. 2010; Gutiérrez et al. 2013). Nevertheless, for low-latitude regions (e.g., tropical Africa or southeast Asia), where the demand for reliable local-scale climate information is of paramount importance due to a large vulnerability to changing environmental factors (Wilby et al. 2009), studies are rare or even nonexistent to date since manifold problems still hinder the successful application of SD in these regions (Hewitson et al. 2014).
In the so-called perfect-prog approach, local-scale climate variability, typically represented by (gridded) weather station records, is statistically linked to quasi-observations from reanalysis datasets (Marzban et al. 2006). The success of perfect-prog schemes in the extratropics relies on the fact that a large fraction of local-scale climate variability can be described by atmospheric phenomena operating on spatial scales on the order of thousands of kilometers, typically having a lifetime of several days. At this scale, reanalysis datasets are known to be skillful, in the sense that their spatiotemporal resolution captures the relevant processes such as extratropical cyclones and the associated fronts (Grotch and MacCracken 1991; Widmann et al. 2003). At lower latitudes, however, the atmospheric drivers of local-scale climate variability operate on much finer scales (both spatial and temporally) and are generally poorly captured by reanalysis datasets (Manzanas et al. 2014). Moreover, observational coverage is generally sparse in the tropics, leading to considerable differences between distinct reanalyses (Trenberth et al. 2001; Sterl 2004; Brands et al. 2012, 2013) and to errors with respect to observational records (Manzanas et al. 2014), which in turn can complicate the detection of a relationship with the local-scale climate variability.
Therefore, the present work tests whether reanalysis choice is relevant for the application of SD in the tropics. A generalized linear model (GLM) approach is separately calibrated for two distinct reanalyses (ERA-Interim and JRA-25) in order to downscale daily precipitation over the Philippines, using a long-term, quality-controlled precipitation dataset that essentially eliminates predictand-induced uncertainty (Hewitson et al. 2014). Because of its geographical location between the monsoonal and inner tropics, the Philippines provides an ideal testbed for SD studies.
First, following a cross-validation scheme for the period 1981–2000, the downscaling results are shown to be sensitive to reanalysis choice if climate change signal–bearing variables such as temperature and/or specific humidity are used as predictors. Second, when the reanalyses-calibrated coefficients are applied to predictor data from a GCM (MPI-ECHAM5)—in which case signal-bearing predictor variables should be applied in order to capture the “correct” climate change signal (Goodess and Palutikof 1998; Wilby et al. 1998)—the sensitivity to reanalysis choice is largely amplified, leading to differences in the projected “deltas” of up to 35% (on average for the entire country) for both reanalyses.
The paper is outlined as follows: In section 2, the considered datasets are described and a brief introduction to the precipitation climate of the Philippines is provided. The applied downscaling technique is described in section 3 and the results are presented through section 4. Section 5 provides the conclusions and a brief discussion on the implications for the entire tropics as well as the model output statistics downscaling approach.
Daily precipitation amounts from 42 gauges maintained by the Philippine Atmospheric, Geophysical and Astronomical Services Administration (PAGASA) were considered as predictand data for the period 1981–2000 (see Fig. 1b). These station time series, which in the following are classified into the four precipitation climate types (CTs) defined in Coronas (1920), were selected after a rigorous quality control, thus minimizing the predictand-induced uncertainty (Hewitson et al. 2014). As can be seen in Fig. 1c, precipitation along the coastlines of the northern part of the archipelago (CT1 and CT2) exhibits a strong seasonal cycle, which is driven by alternating monsoonal winds. In particular, during the southwest monsoon (June–September), precipitation peaks at the stations pertaining to CT1 while CT2 is affected by relative dryness. However, the opposite is the case during the northeast monsoon (October–February). During the dry months (March–May), easterly winds prevail, leading to orographic precipitation along the mountain ranges in the east of the archipelago (see Fig. 1a) and to relatively high precipitation amounts for the stations pertaining to CT2. At the stations pertaining to CT3 and CT4 (mainly situated in the center and south of the archipelago), precipitation is bounded to mesoscale dynamics and is not directly driven by the monsoons, leading to a weak seasonal cycle. Additionally, interannual variability is larger for CT1 and CT2 than for CT3 and CT4 (Fig. 1d). For a comprehensive description of the climate in the Philippines, the interested reader is referred to Coronas (1920), Flores and Balagot (1969), and Kintanar (1984) as well as to the PAGASA website (http://www.pagasa.dost.gov.ph/).
Atmospheric variables describing circulation, moisture, and convection are generally considered to be among the most informative predictors for perfect-prog SD of precipitation (Charles et al. 1999; Timbal et al. 2003; Bürger and Chen 2005; Cavazos and Hewitson 2005; Dibike and Coulibaly 2005; Haylock et al. 2006; Hewitson and Crane 2006; Fowler et al. 2007; Hertig and Jacobeit 2008; Timbal and Jones 2008; Sauter and Venema 2011). If SD is applied in climate change conditions (i.e., to predictor data obtained from a given GCM), the GCM is assumed to perfectly reproduce the same climatological properties provided by the reanalysis used for calibration (Hewitson and Crane 1996; Wilby et al. 2004). In other words, the “performance” of the GCM (Giorgi and Mearns 2002) must be evaluated for the relevant predictor variables. An important dilemma of perfect-prog SD is that GCMs generally perform better for circulation and temperature variables than for moisture ones (Räisänen 2007; Brands et al. 2011, 2013). Yet, moisture information should be included into the predictor field in order to 1) improve the statistical link-function (i.e., the predictive potential of the SD method) and 2) capture the “correct” climate change signal (Goodess and Palutikof 1998; Wilby et al. 1998).
With these precepts in mind, and after consulting the expertise of local meteorologists as well as the results from previously published studies (Kang et al. 2007; Chu et al. 2008; Paul et al. 2008; Chu and Yu 2010), a set of different predictor combinations was chosen (see Table 1). These combinations consist of circulation variables alone (zonal wind component at 850 and 300 hPa; P1: U850, U300), circulation and specific humidity (P2: U850, U300, Q850), circulation and temperature (P3: U850, U300, T850) and circulation, specific humidity, and temperature (P4: U850, U300, Q850, T850). In addition, Q850 and T850 were considered as single predictor variables. Note that Q850 is used instead of column integrated water vapor or precipitable water since the latter variables are usually not provided by the common GCM databases.
The predictor variables listed in Table 1 were obtained from two distinct reanalyses and one GCM: The European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-Interim reanalysis (Dee et al. 2011), the Japanese 25-year Reanalysis (JRA-25) (Onogi et al. 2007), and the Max Planck Institute (MPI) ECHAM5 GCM (Giorgetta et al. 2006); see the acknowledgments for data sources. For the case of ECHAM5, control and A1B scenario data from the third transient run developed within the ENSEMBLES project were retrieved. To keep consistency between the time steps available for both reanalyses and the GCM, daily instantaneous values at 0000 UTC were chosen in all cases. Because of distinct native resolutions, the predictor data from all sources were regridded onto a common regular 2° grid using bilinear interpolation.
3. Downscaling technique
The downscaling technique used here to build transfer functions from the predictors to the predictand (y) is based on generalized linear models (Nelder and Wedderburn 1972), which allow for nonnormal error distributions. The conditional expected value of the predictand given the predictors is linked via a monotonic function to a linear combination of the predictors , where are the regression coefficients. These models have been used in numerous previous downscaling studies dealing with precipitation (e.g., Brandsma and Buishand 1997; Chandler and Wheater 2002; Abaurrea and Asín 2005; Fealy and Sweeney 2007; Hertig et al. 2013).
In this work, the two-stage implementation commonly used for precipitation downscaling is applied (see, e.g., Chandler and Wheater 2002). First, a GLM with Bernoulli error distribution and logit link-function (also known as logistic regression) is used to downscale daily precipitation occurrence (a threshold of 0.5 mm was used to define occurrence). Probabilities equal or greater (smaller) than 0.5 are considered as rainfall occurrences (absences). Second, a GLM with gamma error distribution and log link-function is applied to downscale daily precipitation amount. Unlike in other studies, the stochastic component of the GLM is excluded from each of the two models (occurrence and amount); that is, expected values are predicted in any case. This is done to isolate the effect of reanalysis uncertainty on the downscaling results.
For each gauge, predictor data at the four nearest grid points are considered for both the occurrence and amount models. For the case of the reanalyses and the GCM in the control period, each predictor variable is standardized grid box by grid box to have zero mean and unit variance. Standardization brings the first- and second-order moments of the reanalysis and GCM data into agreement and thereby provides a better approximation for the assumption of “perfect” GCM performance than using untransformed data. The GCM scenario data are standardized by removing the mean of the control period from the mean of the corresponding scenario period and dividing by the standard deviation of the control period.
To avoid overfitting, a k-fold cross-validation approach (Gutiérrez et al. 2013) was followed, with nonoverlapping test periods of five years each, covering the full period 1981–2000. To circumvent spurious trend effects, the five years forming each test period were randomly chosen.
a. Reanalysis differences in the predictor data
The top panel of Fig. 2 shows a comparison between ERA-Interim (taken as reference) and JRA-25 for the four predictor variables in Table 1 over the Coordinated Regional Climate Downscaling Experiment (CORDEX)-East Asian domain (http://wcrp-cordex.ipsl.jussieu.fr/images/pdf/cordex_regions.pdf) for the period 1981–2000. The left column shows the mean difference (bias) between both reanalyses, expressed as a percentage of ERA-Interim’s standard deviation. The middle column displays the ratio of variances (RV), defined as , where () is the variance of JRA-25 (ERA-Interim), respectively. In the right column, the Pearson correlation coefficient (r) between the two reanalyses is depicted.
As can be seen, there are appreciable differences (systematically lower for U850 and U300 than for Q850 and T850) between both reanalyses for the three validation measures considered, indicating that the perfect-prog assumption (reanalysis data reflecting “real” large-scale atmospheric conditions) does not hold for the area under study. Nevertheless, with respect to their application for SD, recall that the reanalysis time series are standardized to have zero mean and unit variance before “entering” the downscaling scheme (section 3). Consequently, differences in the mean and variance between the two reanalyses (left and middle columns) do not affect the SD results, whereas differences in the third- and fourth-order moments (i.e., skewness and kurtosis; see, e.g., Brands et al. 2011) and in day-to-day variations (right column) remain and are expected to affect them.
In the bottom panel of Fig. 2, the zonally averaged r between the predictor time series from JRA-25 and ERA-Interim is displayed for the specific case of the Philippines archipelago. The gridbox coordinates are mapped on the left-hand side and r as a function of latitude is displayed on the right-hand side. Noticeably, U850 exhibits values around 0.95 at all latitudes, which indicates that both reanalyses are in nearly perfect agreement for this variable. However, a north–south gradient is found for the remaining variables. In particular, correlations for T850 and Q850 drop from 0.95 to 0.70 and from 0.75 to 0.50, respectively, probably reflecting the increasing influence of subgrid processes—subject to reanalysis/model-dependent parameterization schemes—toward the equator.
b. Differences in cross-validation results
Figure 3 displays the Spearman correlation coefficient () between daily observed and downscaled precipitation time series over the period 1981–2010 for different predictor combinations—P1 (U850, U300), Q850, T850, and P4 (U850, U300, Q850, T850)—when considering predictor data from ERA-Interim and JRA-25 (solid and dashed lines, respectively). In each panel, the results for a specific CT are shown. Lines and error bars correspond to the mean and standard deviation of the cross-validation results (computed upon the four folds considered). Along the x axis, stations are sorted by decreasing latitude (from left to right). On the right-hand side, the CT-averaged are indicated. Points (asterisks) correspond to ERA-Interim (JRA-25) predictor data.
For both reanalyses, the combination of circulation, humidity, and temperature predictors (P4) yields highest correlation coefficients. The predictive potential is slightly lower for using circulation variables alone (P1) and further decreases if circulation is excluded from the predictor field, that is, for using Q850 and T850 separately or in combination (the latter not shown).
Moreover, for the sole use of circulation variables (P1), the downscaling results are generally not sensitive to reanalysis choice, except for the stations situated in the south (CT4). This is in agreement with the small differences found between ERA-Interim and JRA-25 for U300 and U850 as well as with the slight north–south uncertainty gradient detected for U300 (see Fig. 2). However, for Q850 and T850, appreciable reanalysis-induced differences are observed. In particular, Q850 from ERA-Interim yields better results than Q850 from JRA-25, whereas the opposite is the case for T850 (with the exception of CT1). This indicates that the “real” statistical relationship between Q850 (T850) and local-scale precipitation is more accurately captured by ERA-Interim (JRA-25). Moreover, when considering the “best” predictor combination (P4), results are systematically better for ERA-Interim than for JRA-25. Notably, the southward loss of predictive potential occurring in all CTs except CT2 is in agreement with the southward increase of reanalysis uncertainty (Figs. 2 and 3).
For a geographical overview of these results, Fig. 4 shows the mean pointwise cross-validation rs when considering ERA-Interim (left column) and JRA-25 (middle column) predictor data, with each row corresponding to a specific predictor combination. The corresponding differences—JRA-25 minus ERA-Interim—are displayed in the right column, so positive (negative) values indicate that JRA-25 (ERA-Interim) is more appropriate for SD. Because of the lower predictive potential described above, results for the single predictor variables (Q850 and T850) are not included in Fig. 4.
The spatial pattern of predictive potential is similar for the four predictor combinations. Highest rs values are obtained in the north and along the eastern coastline of the archipelago, whereas a gradual decrease is observed toward the south. For circulation predictors only, both reanalyses perform similarly (first row). However, if Q850 (T850) is added to circulation, better results are obtained for ERA-Interim (JRA-25) (second and third rows, respectively). Notably, for the case of including T850, the advantage of JRA-25 over ERA-Interim is most obvious along the eastern coastline. When considering the “full” predictor combination (P4), ERA-Interim systematically outperforms JRA-25 at all stations.
To further assess the increase in predictive potential from adding temperature and moisture information to circulation, Fig. 5 shows the difference in rs  obtained when adding Q850 and T850 separately (P2 and P3, respectively) and in combination (P3) to the “basic” circulation variables (P1). Results for calibrating with ERA-Interim and JRA-25 are given in the left and middle column, respectively. Additionally, the corresponding differences—JRA-25 minus ERA-Interim—are shown in the right column. Positive (negative) values indicate a larger increment for JRA-25 (ERA-Interim).
In congruence with Figs. 3 and 4, the performance improvement attained when adding Q850 (T850) is larger for ERA-Interim (JRA-25) than for JRA-25 (ERA-Interim). Moreover, when including Q850 + T850, the improvement is larger for ERA-Interim than for JRA-25. The previous results prove that, depending on the choice of reanalysis, up to 0.10 correlation points can be missed on the local scale for particular predictor combinations.
c. Differences in climate change projections
In this section it will be shown that local-scale climate projections obtained by SD are sensitive to the choice of reanalysis used for calibration. To this aim, the regression coefficients obtained from separately calibrating the GLMs with either ERA-Interim or JRA-25 are applied to predictor data from MPI-ECHAM5. This is done for the reference period 1981–2000 (using control run data) and for three different future periods (2011–40, 2041–70, and 2071–2100), using scenario run data (A1B, run 3). The underlying assumption of this procedure is that the predictor–predictand relationships obtained above remain stationary in time (Vrac et al. 2007).
The climate change projections are obtained by means of the delta method, that is, by subtracting the reference/control period’s mean from the mean of the corresponding target scenario period (Räisänen 2007). Deltas are shown as relative (%) deviations from the mean in the reference period (0% = no deviation).
Figure 6 shows, from left to right, three panels, one for each of the future periods considered. In each panel, the deltas projected by applying the coefficients learned from ERA-Interim (JRA-25) are shown in the left (middle) column, while the corresponding differences (JRA-25’s delta minus ERA-Interim’s delta) are provided in the right column; each row corresponds to a particular predictor combination. The numbers in each map indicate the spatial mean value for all stations (All) or those stations pertaining to a specific CT (CT1–CT4, respectively).
As can be seen, a negligible delta is found for any future period if precipitation is downscaled from circulation variables alone (first row). Note that this is in agreement with the time evolution of U850 and U300, which is virtually constant throughout the whole twenty-first century (first and second rows in Fig. 7), indicating that the large-scale circulation (as simulated by MPI-ECHAM5) over the target region is not sensitive to climate change.
However, if Q850 and/or T850 are added to circulation, the projected deltas increase as a function of lead time (i.e., are larger for the end of the century; see the second, third, and fourth rows). This holds valid for using Q850 and T850 as separate predictors (not shown). Remarkably, precipitation deltas for Q850 and T850 are larger than for P4, indicating that the inclusion of circulation damps the climate change signal (not shown).
The fact that Q850—either alone (not shown) or in combination with U850 and U300 (P2 in Fig. 6)—leads to the largest delta differences proves that the downscaling results are especially sensitive to reanalysis choice when this variable is included in the predictor field. For P2, reanalysis-induced delta differences reach 35% (45%) for the entire country (CT1) at the end of the century (2071–2100).
Finally, note that the reanalysis-induced differences in the downscaled time series are proportional to the climate change signal imposed by the GCM (cf. Figs. 6 and 7). Also, the magnitude of the projected deltas seems to be related to the cross-validation results of section 4b. In particular, larger deltas are obtained for the “better” performing reanalysis, that is, ERA-Interim (JRA-25) when Q850 (T850) is added to circulation.
5. Conclusions and discussion
In this study, a generalized linear model (GLM) approach is applied to downscale daily precipitation in the Philippines. To explore the effect of reanalysis uncertainty on statistical downscaling (SD), two distinct reanalysis datasets are used to obtain the link functions (regression coefficients) relating the local-scale predictands to the large-scale predictors. When comparing observed and downscaled daily precipitation time series over the period 1981–2000 using a cross-validation scheme, results are found to be sensitive to the reanalysis dataset selected for calibration, which is in agreement with the few previous studies addressing this issue (Koukidis and Berg 2009; Hofer et al. 2012; Park et al. 2013). However, with spatial average (local scale) correlation differences of 0.03 (0.10) at the utmost, this sensitivity is relatively small at this point.
The reanalysis-calibrated coefficients are subsequently applied to predictor data from a global climate model (GCM) in order to generate local-scale climate projections for the whole twenty-first century. In this case, the reanalysis-induced differences detected in present climate conditions are considerably amplified when signal-bearing variables—which are indispensable for capturing the correct climate change signal—are included in the predictor field. In particular, the projected deltas for the end of the century (2071–2100 minus 1981–2000) are found to differ by up to 35% (on average for the whole country) for the two reanalyses considered. Therefore, the choice of reanalysis used for calibration in perfect-prog SD is an important contributor to the uncertainty of local-scale climate change projections and, consequently, should be treated with as much care as other well-known uncertainty sources, such as the choice of GCM or downscaling method (Dibike and Coulibaly 2005; Chen et al. 2012).
Although these conclusions have been deduced for a specific region (the Philippines), they are likely to hold valid for the entire tropics since previous studies point out that reanalysis uncertainty is a general problem at low latitudes, especially for climate change signal–bearing predictor variables on daily time scale (Brands et al. 2012, 2013). Also, the presented results rely on a single GCM and therefore should be reconfirmed with alternative GCMs (with distinct model physics) in future studies. Besides, an exhaustive assessment on the predictive potential of alternative predictor variables in the tropics, such as velocity potential or streamfunction, and on the corresponding differences induced by reanalysis choice might be a useful future task.
Apart from being relevant for perfect-prog SD, reanalysis uncertainty is expected to be equally relevant for the model output statistics approach, in which GCMs are nudged to reanalysis data in order to force them to follow the “observed” large-scale variability (Eden et al. 2012). Here, it has been shown that the “real” large-scale atmospheric variability in the tropics is likely to be misrepresented by reanalyses and, consequently, also by the aforementioned nudged GCMs. Finally, since regional climate models can be nested into different reanalysis datasets, reanalysis uncertainty is also likely to affect the dynamical downscaling approach (Park et al. 2013).
The authors are grateful to the free distribution of the ECMWF ERA-Interim (http://www.ecmwf.int/en/research/climate-reanalysis/era-interim), JMA JRA-25 (http://jra.kishou.go.jp/JRA-25/index_en.html), and MPI-ECHAM5 data (http://cera-www.dkrz.de/WDCC/ui/Compact.jsp?acronym=ENSEMBLES_MPEH5_SRA1B_3_D) and acknowledge PAGASA for the observational data provided. This study was supported by the EU projects QWeCI and SPECS, funded by the European Commission through the Seventh Framework Programme for Research under Grant Agreements 243964 and 308378, respectively. RM also acknowledges the EU project EUPORIAS, funded by the European Commission through the Seventh Framework Programme for Research under Grant Agreement 308291. SB is grateful to the CSIC-JAE-Predoc Program for financial support.