Dynamic downscaling and bias correction of rainfall in the Pampanga River Basin, Philippines, for investigating flood risk changes due to global warming

: The Pampanga River Basin in Philippines suffers from floods due to typhoons or monsoonal rainfall every year. Assessment of changes in flood risk due to global warming is, therefore, an important issue for this flood-vulnerable basin. We studied possible changes in rainfall features under present (1979–2003) and future (2075–2099) climates using dynamic downscaling of MRI-AGCM experiments. The GCM projections were downscaled into a finer resolution of 5 km for hydrological simulation (catchment size of 10,500 km 2 ). The downscaled rainfall overestimated the number of weak rainfall events, which subsequently resulted in overestimation of monthly rainfall. Bias correction was carried out for downscaled rainfall with reference to raingauge rainfall to perform cumulative distribution mapping. The simulation results found that monthly rainfall would change slightly between present and future climate condi -tions, with extreme rainfall very likely to increase in the future. The annual maximum 48 h rainfall with a 50-year return period may increase from 320 mm under the present climate to 470 mm (MRI-AGCM 3.2S) or 530 mm (MRI-AGCM 3.2H) under the RCP8.5 future climate scenario. This increase in extreme rainfall in the future would have a significant impact on this vulnerable river basin.


INTRODUCTION
The Pampanga River Basin in central Luzon Island, the fourth largest basin in the Philippines, plays key roles in socio-economic activities of the country; it is the most important grain belt and provides Metropolitan Manila with 97% of water resource from the Angat dam in the basin. The Pampanga River Basin, however, suffers from floods due to typhoons or monsoonal rainfall at least once a year which cause disruption and damage; about 4,000 houses were totally destroyed in the worst year, and the cost of damage to infrastructure and agriculture was up to US$ 43 and 280 million in the 2000 and 2011 flood events, respectively . The Intergovernmental Panel on Climate Change 5th Assessment Report (IPCC, 2013) reports, "Future increase in precipitation extremes related to the monsoon is very likely in South America, Africa, East Asia, South Asia, Southeast Asia and Australia". Therefore, the impact of global warming on flood risk change is an important issue for this flood vulnerable basin, and assessment of flood risk change is essential to mitigate social damage and prepare measures to increase the resilience of local communities. Jose and Cruz (1999) studied global warming effects on the Angat water reservoir based on three General Circulation Models (GCMs) and a hydrological simulation model. They projected a decrease in water resources from rainfall due to global warming, and also stressed the importance of global warming on the socio-economy of the Philippines for the first time. More recent studies have investigated future changes in precipitation and temperature by downscaling GCM projections with different models: Chotamonsak et al. (2011) used the Weather Research and Forecasting (WRF) model, Wang et al. (2013) used the Consortium for Smallscale Modelling-Climate Limited-area Modelling (COSMO-CLM) in the framework of the Coordinated Regional Climate Downscaling Experiment (CORDEX), and Cruz et al. (2016) used the Non-Hydrostatic Regional Climate Model (NHRCM). However, all these studies were conducted with models having a downscale resolution of about 50 to 60 km and their discussion focused on monthly average precipitation. Because 50-km resolution models are only capable of resolving rain systems of more than 200 km, those models could not resolve extreme rainfall such as typhoon rainbands.
To examine flood risk and changes due to global warming, rainfall information with much finer horizontal and temporal resolutions is desired. In this study, dynamic downscaling into 5 km horizontal resolution was conducted for the first time in this area, since this finer downscaling has been found effective to improve the reproduction of heavy rainfall (Kanada et al., 2008). This level of resolution is also fine enough for streamflow and inundation simulation for this river basin (catchment size of 10,500 km 2 ). The boundary conditions were provided by the Meteorological Research Institute-Atmospheric General Circulation Model (MRI-AGCM) version 3.2. We further performed bias correction of the downscaled rainfall, because rainfall includes significant biases even after downscaling. The effect of global warming is discussed using extreme rainfall cases.

DATA
We used data from raingauge observation for validation of downscaled rainfall, and GCM projections for the boundary conditions of a regional climate model. We obtained hourly raingauge data collected in the Pampanga River Basin for 22 years over the periods of 1980-1982, 1993, and 1995-2012, including frequent periods with data missing. The number of raingauges was 9 at the beginning and increased up to 17 at the end.
The boundary conditions of dynamic downscaling were given by MRI-AGCM 3.2S (super high resolution in 20 km) and 3.2H (high resolution in 60 km) (Mizuta et al., ) for present (1979(Mizuta et al., -2003 and future (2075-2099) climate projections based on the RCP8.5 scenario. The model was designed to have quite a high resolution compared with other Coupled Model Intercomparison Project Phase 5 (CMIP5) GCMs, and brought the advantage of better reproduction of mean precipitation in the monsoon area (IPCC, 2013). Since MRI-AGCM is an atmospheric GCM, sea surface temperature (SST) was given by monthly mean observations (present) or observations plus future changes derived from the CMIP5. MRI-AGCM 3.2S was provided with four different SST distributions: multi-model ensemble (MME) and clusters 1, 2 and 3 (Kitoh and Endo, 2016). MRI-AGCM 3.2H was provided with only one future projection with MME SST. The surface variables were given in original resolutions of 20 and 60 km. However, three dimensional variables (U; zonal wind, V; meridional wind, T; temperature, Q; specific humidity, Z; geopotential height) were given in 1.25 degree resolution in 12 layers from 1000 to 100 hPa.

METHODS
Dynamic downscaling was performed using the Weather Research and Forecasting (WRF) model ver. 3.4.1, which is a next generation mesoscale numerical weather prediction model developed by the National Center for Atmospheric Research (NCAR) and several other research institutes in the United States (Skamarock et al., 2005). The model domains were set up with double nesting frames with 67 × 67 × 40 grids in both outer and inner frames at 15 or 5 km horizontal grid intervals and stretching vertical grids ( Figure 1). For parameterization schemes, the WRF single moment 3-class (WSM3) microphysics scheme, thermal diffusion surface physics, the Mellor-Yamada-Nakanishi-Niino level 2.5 planetary boundary layer scheme were adopted. The sea surface temperature used for the integration of GCMs was given as the lower boundary condition. To select a cumulus parameterization scheme, which is critical for production of precipitation, we tested three different types: the default Kain and Fritsch scheme (Kain and Fritsch, 1993), the Kain and Fritsch scheme tuned by the Japan Meteorological Agency (JMA) (Narita, 2008), and the Grell 3D ensemble scheme with a shallow convection option (Grell and Devenyi, 2002). After comparing rainfall simulated using the three schemes for a period of several years, we adopted the third scheme, as it demonstrated the best performance, producing rainfall having good agreement with extreme rainfall from raingauge observations. The downscaling calculations were carried out separately with respect to each year from January to December. We did not include any additional spinup periods. However, we assumed that the first week was negligible, since one-week spinup time was long enough and the beginning of each year was in the driest season. Figure 2 compares average rainfall from raingauge observations and WRF downscaling during four recent major flood events in the Pampanga River Basin. In the WRF simulation, the boundary condition was provided by ERAinterim reanalysis. The simulated rainfall variation in the September 2009 and September 2011 cases agreed well with that of the raingauge observations, although the rainfall peaks during the September 2009 event were overestimated. These two events were due to typhoon passages. In the June 2011 case, the simulated rainfall trend did not necessarily match the observed one, but it did reproduce rainfall peaks of similar magnitudes. In the August 2012 case, the simulation showed a fairly good reproduction of rainfall that was superimposed on diurnal variation. The latter two events were caused by the southwesterly monsoon. The comparison in Figure 2 shows that the WRF model is reliable in reproducing typhoon rainfalls, and it was also capable of representing the basic characteristics of monsoonal rainfalls, suggesting the validity of using the WRF model for our study on extreme rainfall.  Figure 3 shows seasonal variation in rainfall and frequency of appearance, before and after downscaling, and after bias correction. Figure 3a shows GCM rainfall where bias is corrected based on APHRODITE data (Yatagai et al.,Figure 3. Seasonal variation in basin average rainfall (a)-(c), and frequency of appearance of daily rainfall at raingauge sites (d)-(f). The first column shows rainfall from GCM, the second column shows rainfall after downscaling, and the third column shows rainfall after downscaling and bias corrected based on raingauge data. The GCM rainfall in (a) were bias corrected based on APHRODITE rainfall. In (a)-(c), the black line indicates observations (APHRODITE in (a), and raingauge in (b) and (c)), the blue line indicates present climate projections, and the red to yellow lines shows future climate projections in MRI-AGCM 3.2S (the solid line is MME SST, the long broken line is cluster 1, the short broken line is cluster 2, and the dot broken line is cluster 3). The vertical bars represent standard deviations. In (d)-(f), the black marks indicate raingauge, and the red and green marks indicates present climate projections of MRI-AGCM 3.2S and 3.2H, respectively  Figure 3b, on the other hand, the downscaled rainfall is overestimated compared to observed rainfall. Figure 3c shows rainfall with biases corrected based on raingauge observation, as explained below. The black line in Figure 3c is a basin average rainfall based on the raingauge, which is slightly smaller than the one in Figure 3a by APHRODITE, but shows quite similar variation. It was thus confirmed that the raingauge observation is equivalent with APHRODITE.

RESULTS
The dynamic downscaling dramatically improved the appearance of heavy rainfall from GCM, as illustrated in Figures 3d and 3e. Figure 3d shows that the GCM-based rainfall was quite a poor reproduction of frequency of appearance from raingauge data. However, Mizuta et al. (2012) have already pointed out the possibility of GCM's poor performance particularly in its application to tropical western Pacific islands (see Figure 16e in their study). The poor performance by MRI-AGCM in the case of the Philippine islands could be interpreted in this respect.
The downscaling overestimated the number of weak rainfall events with less than about 100 mm day -1 (Figure 3e), compared with raingauge observation. This overestimation of weak rainfall was a reason for overestimated monthly rainfall in Figure 3b. To correct this disagreement, we employed bias correction based on popular distribution mapping (De Troch et al., 2013). As shown in Rasmy et al. (2014), distribution mapping alone is appropriate for bias correction of downscaled rainfall, although a 3-step bias correction works better for GCM rainfall (Jaranilla- Sanchez et al., 2013;Nyunt et al., 2014). Because of the limited number of raingauge stations and incomplete data continuity, we mapped the downscaled basin-average rainfall onto basin average raingauge rainfall. The cumulative rainfall probability was used without any fitting on a Gamma function, since this non-parametric method performs the best in reproducing extreme rainfall. The mapping was performed for each month, but the top 0.5% rainfall events were mapped for each year since annual maximum extreme rainfall can appear in any month of the year (Inomata et al., 2009). The obtained correction rate was multiplied by the downscaled rainfall, and spatial distributions were recovered. The correction rates used for the present climate were also applied to future climate rainfall. After bias correction, seasonal variations ( Figure 3c) agreed well with observations, and the frequency appearance also showed consistency with raingauge observation (Figure 3f). Figure 3c compares the seasonal variation of rainfalls from observation and downscaling of present and future climate projections. Rainfall in the wet season (May to September), especially cluster 2 SST (orange short broken line), slightly increased due to global warming, but the increase was not significant. Cluster 2 SST is an El Niñotype pattern (Kitoh and Endo, 2016), in which SST cools slightly around the Philippines. Although it was not clear why cluster 2 SST projected the largest rainfall, we also found this tendency in other analyses. As Figure 3 shows, only a slight change in monthly rainfall is projected in the Pampanga River Basin due to global warming. Figure 4 indicates the results from frequency analysis of extreme rainfall. We chose the annual maximum 48 h rainfall for this analysis, since our flood inundation simulation area had the highest correlation with maximum 48 h rainfall. Figure 4a shows that the frequency distributions of raingauge rainfall, ERA-interim, and the present climate of MRI-AGCM 3.2S were consistent with each other. In Figure  4b, the present climate of MRI-AGCM 3.2H estimated a slightly larger rainfall for longer return periods without a significant degree of overestimation. This indicates that the bias corrected downscaled rainfall events represent the present climate rainfall features successfully. Figure 4. The frequency analysis of annual maximum 48 h rainfall from downscaling and bias corrected results of (a) MRI-AGCM 3.2S and (b) MRI-AGCM 3.2H. The lines are Gumbel distribution fitting. The points and lines in black are raingauge observations, the green ones are downscaled ERA-interim, the blue ones are present climate projections, the red ones are future climate projections with MME SST, the orange long broken line is future climate projection with cluster 1 SST, the orange short broken line is cluster 2 SST, and the orange dot-broken line is future climate projection with cluster 3 SST In Figure 4, the four lines of red to orange colors show the frequency distributions under the future climate. The frequency distributions with MME and cluster 2 SST in Figure  4a and MME SST in Figure 4b were significantly larger than those of the present climate. In the case of the 50-year return period, rainfall is estimated to be 320 mm in the present climate, while it is projected to increase by 45% to 470 mm (MRI-AGCM 3.2S, cluster 2) and by 65% to 530 mm (MRI-AGCM 3.2H) in the future climate. The results suggest a drastic increase in extreme rainfall in the Pampanga River Basin due to global warming. However, the two lines in MRI-AGCM 3.2S (clusters 1 and 3) in Figure 4a show little change due to global warming, suggesting uncertainty in the range of increase.
To confirm the credibility of the extreme rainfall changes under the future climate, we further examined those amounts by original MRI-AGCM 3.2S and 3.2H without downscaling and without bias correction in the same way as in Figure  4 (figures are not shown). The maximum 48 h rainfalls obtained in the 50-year return period were 329 mm and 317 mm under the present climate for MRI-AGCM 3.2S and 3.2H, respectively, and they increased to 551 mm and 461 mm, respectively, under the RCP8.5 scenario in MME SST. The rates of increase in extreme rainfall due to global warming by original MRI-AGCM 3.2S and 3.2H were 67% and 45%, respectively. Thus, the rates of increase in extreme rainfall, 45% and 65%, obtained from the downscaled rainfall are reasonable based on MRI-AGCM 3.2S and 3.2H.
Our analysis suggests a potential increase in extreme rainfall due to global warming. This is consistent with past studies (IPCC, 2013;Jaranilla-Sanchez et al., 2013;Kitoh and Endo, 2016). The advantage of this study is the use of a downscaling model with a higher resolution. Jaranilla-Sanchez et al. (2013) showed that, based on bias-corrected rainfall from six GCMs, the 10th percentile peak daily discharge may increase in the Pampanga River Basin from the years of 1981-2000 to 2046-2065 in the range of zero to six times. It is not appropriate to directly compare their results with ours. However, our analysis using a downscaling approach has contributed to reduction of the uncertainty concerning a potential increase in extreme rainfall due to global warming.
To discuss the mechanism of rainfall increase due to global warming, monsoon strength was compared between the present and future climate in Figure 5. The monsoon trough along 20°N is dominant in the present climate, driving the south-westerly monsoon in the summer time ( Figure  5a). In contrast, the monsoon trough may weaken in the future climate, resulting in a weaker monsoon flow ( Figure  5b) and a resultant decrease in west coast rainfall (figure not shown). The weakening of monsoon circulation is also mentioned in the IPCC AR5 report (IPCC, 2013). All this implies that, at least, the average monsoon flow may not strengthen the projected extreme rainfall.
In this river basin, the forcing of extreme rainfall is due to either typhoons or monsoon. Figure 6 compares the causes of annual maximum rainfall in each year in each scenario, which was determined subjectively based on temporal variation of rainfall distributions. Extreme rainfalls can be attributed to typhoons in 27 out of 50 years (54%) between 1979 and 2003 in the present climate, while in 30 out of 125 years (24%) between 2075 and 2099 in the future climate. This implies that the number of typhoon-induced rainfall events may significantly decrease due to global warming. This implication is consistent with many previous studies (e.g., Murakami et al., 2012) reporting a potential decrease in the number of typhoons in the Western Pacific due to global warming, although the number of very strong typhoons is projected to increase. Nevertheless, very intense extreme rainfall events tend to be attributed to typhoons in most cases, as Figure 6 indicates, which suggests that typhoon rainfall would still remain important in forcing very extreme rainfall in the future climate.

SUMMARY AND CONCLUSIONS
We investigated the impact of global warming on changes in extreme rainfall in a flood prone river basin of the Philippines under the present  and future (2075-2099) climate, using dynamic downscaling of MRI-AGCM 3.2S and 3.2H and the RCP 8.5 scenario. Since downscaled rainfall overestimates the number of weak rainfall events, bias correction was applied based on distribution mapping with data from raingauge observations. The corrected present climate rainfall showed reasonable agreement with raingauge observation with regard to monthly rainfall, frequency appearance of daily rainfall, and frequency analysis. The monthly rainfall in the future climate did not show any significant changes from that in the present climate. However, the frequency analysis found that extreme rainfall may increase in the future climate, depending on SST conditions. With four SST conditions in MRI-AGCM 3.2S and one SST condition in MRI-AGCM 3.2H, future climate projections suggested a possible increase with large uncertainty, having a wide range from 0 to 65%. The increase in extreme rainfall projected in this study is consistent with past studies. In addition, the study contributed to specifying an increase in range. This rainfall data will be used in future hydrological simulations to examine changes in flood inundation extent due to global warming.
Extreme rainfall events over the Philippines are caused by either typhoons or monsoon flows. The strength of the East Asian monsoon in summer may decrease in the future. In contrast, 54% of annual extreme rainfall events were attributable to typhoons in the present climate, compared to only 24% in the future climate. This means that the monsoon may become weaker, and the number of typhoons may decrease, in the future climate. Further study is necessary to elucidate the mechanism of a potential increase in extreme rainfall. Figure 6. Annual maximum 48 h rainfall. The yellow bars indicate typhoon attributions and the green bars indicate monsoon attributions. The two panels in the blue rectangle are from the present climate, and the five panels in the red rectangle are from the future climate