Regional climate models (RCMs) should be evaluated with respect to their ability to downscale large-scale climate information to the local scales, which are sometimes strongly modulated by surface conditions. This is the case for La Réunion (southwest Indian Ocean) because of its island context and its complex topography. Large-scale atmospheric configurations such as tropical cyclones (TCs) may have an amplifying effect on local rainfall patterns that only a very high-resolution RCM, forced by the large scales and resolving finescale processes, may simulate properly.
This paper documents the capability of the Weather Research and Forecasting Model (WRF) RCM to regionalize rainfall variability at very high resolution (680 m) over La Réunion island for daily to seasonal time scales and year-to-year differences. Two contrasted wet seasons (November–April) are selected: 2000–01 (abnormally dry) and 2004–05 (abnormally wet). WRF rainfall is compared to a dense network of rain gauge records interpolated onto the WRF grid through the regression-kriging (RK) technique. RK avoids the point-to-grid comparison issue, but produces imperfect estimates due to sampling, so its quality also needs to be tested.
Seasonal rainfall amounts and contrasts produced by WRF are fairly realistic. At intraseasonal and daily time scales, differences to RK are more sizable. These differences are not easy to interpret in sectors where the rain gauge network is less dense and the quality of RK more uncertain, as over the eastern slopes of Piton de la Fournaise volcano where WRF seems to simulate more realistic rainfall than RK. Finally, the heavy rainfall associated with TC Ando on 6 January 2001, is documented. WRF shows weak disagreements with RK, indicating its capability to regionalize rainfall during extreme events.
La Réunion island (21°S, 55°E: 800 km off the east coast of Madagascar) is a French overseas department located in the southwest Indian Ocean (SWIO). It is a mountainous island (2512 km2) with very complex terrain (Fig. 1a). La Réunion is also a hotspot volcano with three major central depressed areas (the Cirques of Mafate in the northwest, Salazie in the northeast, and Cilaos in the southwest: see locations in Fig. 1a), which are enormous collapses of La Réunion’s oldest volcano, with edges that are practically vertical walls, more than 1000 m high (Oehler et al. 2008). La Réunion has also two coalescent volcanic peaks: the Piton des Neiges (3070 m) and the Piton de La Fournaise (2632 m, Fig. 1a).
La Réunion island has a tropical maritime climate marked by two seasons depending on the behavior of the southern Hadley cell and the Walker circulation (Baldy et al. 1996). During austral winter from May to October/November, the Hadley–Walker cell circulation strengthens and produces steady low-level easterly winds and upper-level westerly winds (Barcelo and Coudray 1996; Taupin et al. 1999). This vertical shear coincides with a temperature inversion that behaves as a barrier for the development of clouds. During austral summer in contrast, the trade winds weaken and their inversion almost disappears. La Réunion experiences wet conditions, especially from January through March. April and December are transition months, likely to be equally dry or wet.
The rainfall field over the island shows a huge variability in both time and space. Spatially, average annual amounts present a very marked west–east gradient, reaching values larger than 10–12 m in the elevated sectors facing the dominant moisture fluxes associated with the trade winds. Temporally, the climate of La Réunion is one of the most abrupt and violent in the world, the island being placed under the influence of tropical cyclones (TCs) during the summer season. La Réunion has world record rainfall amounts for all time scales ranging from 12 h to 15 days, due to the heavy rainfall associated with TCs, their effects being strongly modulated locally by the topography (Barcelo et al. 1997; Mayes 2007). Such extreme events recurrently cause deaths, floods, and significant damages to infrastructure and crops.
Given these characteristics of the local climate, with frequent TCs having potentially devastating effects, together with the impact of topographic forcing on rainfall amounts at very fine scales, a detailed examination of the space–time variability of rainfall at high resolution over the island is then needed. This necessitates the development of high-resolution gridded rainfall products with continuous spatial and temporal coverage. Information, as given by these products, is of great interest for applications in meteorology, climate, and hydrology. One application, for example, would be to assess the effects of large-scale atmospheric configurations (TCs, synoptic situations, or more generally circulation systems over the SWIO) on the local rainfall patterns.
Météo-France, the French weather forecasting service, is maintaining a dense network of observations covering most parts of the island. In recent years, it was, for example, composed of roughly a hundred of rain gauge stations, making La Réunion an almost unique location to analyze the anisotropy of spatial rainfall distributions in the tropics. Yet, as for all observational networks, it cannot document the variability of the rainfall field continuously in time and even more clearly in space. Indeed, the cumulative area of rain gauge (4 × 10−2 m2 per gauge) leads to a spatial sampling of roughly 1.6 × 10−8 (obtained as the ratio between the surface of the rain gauges and that of the island).
One way for obtaining time–space-evolving information for rainfall at local scales is to apply dynamical downscaling with a high-resolution regional climate model (RCM) forced by global climate model boundary conditions and resolving the finescale surface processes (e.g., Wang et al. 2004). Over the neighboring regions of southern Africa and the SWIO, many studies documented the skill of current RCMs to simulate the regional climate mean state and variability. Investigated time scales spun from the case studies of rain-bearing systems [e.g., mesoscale convective systems: Blamey and Reason (2009), or tropical-temperate troughs: van den Heever et al. (1997)], seasonal scale [e.g., Crétat et al. (2011), (2012a); Crétat and Pohl (2012); Ratnam et al. (2013)], and more recently climate scales [Williams et al. (2010); Vigaud et al. (2012); Boulard et al. (2013); Ratna et al. (2014); Pohl et al. (2014)]. However, none of these studies reached very high spatial resolutions (a few hundreds of meters) needed to take properly into account the anisotropic and heterogeneous surface conditions over an island as small as La Réunion (roughly 50 km × 50 km). For example, Ratna et al. (2014) used 9-km grids over South Africa; Pohl et al. (2014) used 18-km grids over South Africa and 3.6-km grids over its northeastern part. Such resolutions are still too coarse to simulate the local climate of La Réunion island properly. To that regard, Yu et al. (2014) used a high-resolution numerical model at 4-km, 2-km, 1-km, and 500-m horizontal resolutions to simulate a 4-day heavy rainfall event over La Réunion. They showed that the grid spacing below 1 km was necessary to represent well the spatial distributions of the 48- and 24-h cumulative rainfall fields and the extreme rainfall amounts at particular locations. However, this study, the only work on the modeling of finescale rainfall that was done on La Réunion itself, solely focused on the short-term rainfall extreme.
In this study, we propose to assess the capability, usefulness, and limitations of a current state-of-the-art nonhydrostatic RCM to regionalize rainfall over La Réunion island at very high resolution (a few hundreds of meters). Time scales considered in this work range from seasonal differences between two contrasted rainy seasons on the one hand, to daily rainfall including one extreme event caused by a TC on the other hand. Seasonal means and intraseasonal variability are also documented. An extensive evaluation of the RCM skill is made possible by the dense network of rain gauge stations provided by Météo-France. Because a direct comparison on a “grid by point” basis may provide intrinsic errors due to the potentially high local variability of daily rainfall on oceanic tropical islands (e.g., Stow and Dirks 1998), a “grid to grid” comparison, which accounts for this local variability, is preferred. To that end, observations are projected onto the same grid as that used by the RCM, using geostatistics. The approach retained here is the so-called regression-kriging (RK) technique (e.g., Hengl 2009; Joly et al. 2011), which allows establishing statistical relationships between large-scale climate conditions and local surface properties in order to estimate the spatial distribution of a given field (e.g., rainfall). The gridded dataset thus created is to become a new alternative “reality” to use for comparison with the model instead of the rain gauge measurements. As a result of the quality of RK fields being intrinsically dependent on the density and representativeness of the observational network, the realism of gridded fields needs to be tested beforehand.
This study should be seen as a preliminary work, based on hindcast simulations that aim to document how large-scale circulation patterns over the SWIO modulate the amounts and spatial distribution of rainfall, at local-terrain scale, over La Réunion. The eventual goal of our work will be to predict the effects of the large-scale circulation patterns on rainfall at the local scale (i.e., fill the gap between the terrain reality and the coarse resolution provided by the general circulation models used for climate studies).
This paper is organized as follows. Section 2 presents the datasets and methodologies used in this work, including the dynamical downscaling and geostatistical upscaling methods. Section 3 is dedicated to the presentation of the results, which are then discussed in section 4 and briefly summarized in section 5.
2. Data and methods
a. Observational data
Two contrasted austral summer seasons are considered: November through April 2000–01 (abnormally dry) and 2004–05 (abnormally wet). Daily rainfall amounts are provided by the Météo-France weather station network. In total 86 and 96 stations are retained for 2000–01 and 2004–05, respectively (Fig. 1a). Data contain a larger amount of missing values in 2004–05 (around 3.1%) than in 2000–01 (less than 0.5%).
b. Dynamical downscaling
Regional experiments are performed using the nonhydrostatic Advanced Research Weather Research and Forecasting Model (ARW-WRF), version 3.3.1 (hereafter WRF; Skamarock et al. 2008). The physical package includes the Kain–Fritsch scheme for atmospheric convection (Kain 2004) using the trigger function developed by Ma and Tan (2009), the WRF single-moment 6-class microphysics scheme (WSM6) for cloud microphysics (Hong and Lim 2006), and the Yonsei University parameterization of the planetary boundary layer (PBL; Hong et al. 2006). Radiative transfer is parameterized with the Rapid Radiative Transfer Model scheme (Mlawer et al. 1997) for longwaves and Dudhia (1989) scheme for shortwaves. Over the continents WRF is coupled with the four-layer Noah land surface model (Chen and Dudhia 2001a,b). Except for elevation data, which are from the Shuttle Radar Topography Mission (SRTM), all other surface data are taken from U.S. Geological Survey (USGS) database, which describes a 24-category land-use index based on climatological averages, and a 17-category Food and Agriculture Organization of the United Nations soil data, both available at 10 arc-min.
Simulations are carried out over two contrasted austral summer seasons: November through April 2000–01 and 2004–05. The regional domains consist of four two-way nested domains (Fig. 1b), at 43.45-km, 10.86-km, 2.72-km, and 680-m horizontal resolutions, all centered on La Réunion island (21.0°S, 55.0°E), and with 28 sigma levels on the vertical. Deep atmospheric convection is explicitly resolved for domains 3–4.
Lateral forcings are provided every 6 h by the Interim European Centre for Medium-Range Weather Forecasts (ECMWF) Re-Analysis (ERA-Interim; Simmons et al. 2007; Dee et al. 2011) at a 1.5° horizontal resolution and 19 pressure levels. SST fields are prescribed every 24 h after a linear interpolation of monthly ERA-Interim SST. Integrations are initialized on 15 October for each year, allowing for a 15-day-long spinup period. Adaptative time stepping (Hutchinson 2007) is used to increase model efficiency. Data are archived every 24 h from 15 October to 30 April. Because our study is meant to be a downscaling exercise, a relaxation (nudging) approach is considered in order to reduce the model’s biases and internal variability (Pohl and Crétat 2014). Spectral nudging toward ERA-Interim horizontal winds (every 12 h), and temperature and specific humidity (every 24 h), is applied throughout the larger domain (domain 1: 39.5°–2.5°S, 8.5°–101.5°E; covering southern Africa, Madagascar, and the southern Indian Ocean) and above the PBL. It retains wavelengths larger than ~1000 km, for which reanalyses are thought to provide realistic solutions. By construction, WRF-simulated large scales are thus close to ERA-Interim, especially above the PBL (Diaconescu and Laprise 2013). In the following, ERA-Interim solutions are thus solely used to characterize the large-scale atmospheric configurations associated with rainfall variability at local scales over La Réunion.
c. Geostatistical upscaling
In this study, simulated rainfall is compared against observations. However, to avoid the grid-to-point comparison issue, observed rainfall data are interpolated onto WRF’s finest grid (domain 4; Fig. 1a) using the RK method (Odeh et al. 1995). RK, mathematically equivalent to kriging with external drift (KED; Hengl et al. 2007), is suitable as a nonstationary model for mapping rainfall (Lloyd 2010). This method generally performs well in areas of complex topography (e.g., Prudhomme and Reed 1999). We briefly introduce here the main features of RK; for further details, see Hengl (2009).
A starting point for the method is the idea of Matheron (1969) of a universal model of spatial variation, where a value of a target variable (e.g., rainfall) at some location can be modeled as the sum of a deterministic component (trend) and a stochastic component (residuals). In the case of RK, estimations are made separately for the trend and the residuals and then added back together.
Let a set of observations of a target variable be denoted as z(s1), z(s2), …, z(sn), where si = (xi, yi) is a location and xi and yi are the coordinates in geographical space, and n is the number of observations. Thus, the value of the target variable at a new location s0 is estimated as follows:
where the trend, (s0), is fitted using multiple linear regression analysis, and the interpolated residual, (s0), is estimated using ordinary kriging (Cressie 1993).
If βk are the regression model coefficients, qk are the predictors, p is the number of predictors, λi are the kriging weights, and e(si) = z(si) − (si) is the residual at location si, the prediction (s0) is made by
Traditionally, regression coefficients βk are fitted using the ordinary least squares (OLS) procedure, which works well as long as the predictor variables are fairly few and uncorrelated. When some of the predictor variables are highly correlated, OLS tends to become inefficient. A stepwise selection process (Hocking 1976) can be used to reduce the colinearity of predictors by limiting their number depending on their statistical significance in the regression (e.g., Oettli and Camberlin 2005; Joly et al. 2011).
In this study, linear regression analysis is carried out between rainfall and 10 candidate predictor variables. These are six commonly used, temporally invariant predictors: longitude, latitude, elevation and its derivatives (slope, aspect), along with distance to the coastline. To incorporate the conjoined effects of lower-atmosphere characteristics, we also use four time-dependent low-atmosphere state variables: horizontal wind components at the 10-m level, specific humidity integrated from 1000 to 700 hPa, and vertical wind component due to orographic lifting. These variables are derived from ERA-Interim reanalysis daily products based on an original approach similar to that used by Kyriakidis et al. (2001). In this work, we could have used local surface observed winds from about 40 weather stations. However, because these in situ measurements are representative of local flows forced by topography, they are not the important factor controlling rainfall amounts. ERA-Interim reanalysis variables are available at a resolution of 0.75° × 0.75°. We retain reported values of specific humidity and horizontal wind components at 225 reanalysis nodes within a 10° × 10° domain centered on La Réunion roughly corresponding to WRF domain 2 (Fig. 1b). These values are then interpolated onto the same digital elevation model as used in WRF’s domain 4 grid (shown in Fig. 1a). The vertical wind component is finally computed as the inner product of directional elevation gradients with the interpolated 700-hPa horizontal winds [see Kyriakidis et al. (2001) for details]. Although domain 4 only occupies a little bit less than four ERA-Interim grid points, the predictors thus derived provide a smooth picture of the large-scale state of the lower atmosphere over the study region, which is expected to be related to observed rainfall at the local scale. This hypothesis is formulated because rainfall amounts over La Réunion are mainly linked to the large-scale low-atmosphere moisture fluxes associated with the dominant trade winds (see the introduction). Though the regression is purely spatial, because ERA-Interim reanalysis variables are daily, the temporal variability can be addressed via averaging over the time periods of interest (a season, a subset of days or cluster, or a single day), allowing for the predictors to change over time. To the authors’ knowledge, the inclusion of dynamic predictors in this way has not been reported previously.
Kriging weights λi are derived from the semivariogram model defining the spatial autocorrelation structure (Matheron 1971; Cressie 1993; Isaaks and Srivastava 1989; Goovaerts 2000). The semivariogram model involves semivariance γ (h), which is estimated from the experimental data as half the average squared difference between the paired data values:
where zi = z(si) is the value of the target field at some number of locations si. The sum is over all pairs of points a distance h apart, and n(h) is the number of such pairs.
As mentioned previously, rainfall data from 86 and 96 stations are used for 2000–01 and 2004–05. The RK procedure is applied only to stations with complete time series. Finally, errors associated with the RK-derived rainfall are computed from leave-one-out cross validation (LOOCV).
d. Choice of the predictors for seasonal rainfall
RK maps of seasonal rainfall are generated based on stepwise selection of predictors in the spatial variation of the regression model of each period. This approach leads to combining 7 to 8 among the 10 available explanatory variables, among which some are linearly related (e.g., altitude and distance to the sea are correlated at roughly +0.78, giving thus the simplified picture of an imperfect, slightly elongated cone for the island). Four out of the six geographic–topographic variables (latitude, longitude, distance to the sea, and slope) are selected in the two seasonal regression models. Aspect is never selected. Seasonal rainfall is more sensitive to longitude; depending on the period, the Pearson correlation coefficient r varies from 0.77 (2000–01) to 0.64 (2004–05). The effect of latitude is secondary, with inconstant signs. The slope predictor amplifies seasonal rainfall. Distance to the sea and altitude form a binomial: the presence (absence) of altitude in the 2004–05 model (2000–01) impacts the sign associated with distance to the sea. Atmospheric variables also influence rainfall, especially specific humidity, selected as an explanatory variable for both periods. There are also strong associations between seasonal rainfall and the zonal and meridional components of the 10-m wind (u10 and v10, respectively). The presence of the pair u10–v10 is not systematic and the associated signs can change (this is the case for v10). This suggests that the wind speed is more important (in term of impacts on the seasonal rainfall) than its direction. Finally, the resulting equations are
The adjusted R2, quantifying the fraction of the variance explained by the predictors in the regression models, reaches 82.4% (71.4%) for 2000–01 (2004–05). The next step consists in kriging the residuals. Given that a Gaussian distribution of data is preferred for kriging (Schuurmans et al. 2007), we first verify that the residuals are normally distributed and then use a Gaussian function to adjust the experimental semivariograms (not shown). The root-mean-square error (RMSE) in LOOCV is lower for RK than pure regression; the decrease in RMSE is higher in 2004–05 (from 0.658 to 0.475 m: −27.8%) than in 2000–01 (from 0.418 to 0.391 m: −6.5%). This is not surprising given the higher R2 value from regression for 2000–01, indicating that not much gain was to be expected by combining regression and kriging for this period.
e. Choice of the predictors for daily rainfall
An attempt is made to use the same procedure in order to produce maps of high-resolution daily rainfall. This is a more difficult task for several reasons: 1) daily rainfall patterns are often weakly coherent spatially, especially in regions where convective rainfall is prevalent (Moron et al. 2007); 2) internal variability in an RCM is larger at local temporal and spatial scales (Crétat et al. 2011); and 3) probability density functions of daily rainfall are often biased in climate models (Crétat et al. 2012a).
For each day of the period, the same stepwise selection of predictors is performed. Figure 2 shows box-and-whisker plots of regression coefficients series for 2004–05, along with the number of days (out of 181 days) in which each of them was retained in the stepwise selection process. Roughly similar results are found for 2000–01 (not shown). The total number of predictors selected each day varies from 0 (i.e., pure kriging is applied) to at most 9, with mean value around 5. Predictors the most frequently selected include atmospheric and geographic predictor variables, in contrast with topographic variables (e.g., elevation is only selected in 40% of the cases). In terms of number of occurrences, 10-m level υ wind (68% of the cases), latitude (67%), longitude (64%), specific humidity (60%), and 10-m level u wind (58%) are the most important. Although no work has attempted to evaluate the relative skill of different atmospheric predictors for rainfall amounts over La Réunion to date, results here confirm our a priori choice strategy, partly based on dynamic predictors (section 2c).
The influence of the selected factors may be very different from one day to another. For a first group of predictors comprising the altitude and the topographic variables, the sign of the underlying regression coefficients is the same most of the days; only their amplitude strongly varies from one day to another. For a second group comprising the latitude, the longitude, and the atmospheric variables, the sign of the regression coefficients is strongly variable (positive in about half the cases and negative otherwise). The first group indicates the relief effects independently of the synoptic-scale airflow direction, and the second group, the directional nature of the rainfall pattern in relation to the overall configuration of the atmospheric flows. Depending on the number of selected variables, the explained variance varies significantly, from 0% to roughly 90% at the most, with mean value of 37%.
f. Daily rainfall classification procedure
To document intraseasonal variability without producing two series of 181 daily maps, the most resembling days are objectively regrouped together using a hierarchical agglomerative technique (Crétat et al. 2012b), hereby forming statistically homogeneous classes. The ascending technique means that, at the start of the algorithm, each day constitutes a separate class. Then, iteratively, it merges two by two all the classes (i.e., days or groups of days), by identifying the two most similar ones. The algorithm can be stopped at any step, before the classes are finally merged to form a single class. This nonparametric clustering technique, for which there is no criterion for the choice of the number of clusters, also has the advantage over other techniques based on k-means that it is applicable to non-Gaussian series like daily rainfall. Cheng and Wallace (1993), Martineu et al. (1999), Casola and Wallace (2007), or more recently Crétat et al. (2012b) have demonstrated the usefulness of this technique for the identification of weather regimes.
The hierarchical classification is based on the distance matrix Ψ, calculated as the Euclidean distance between each observation (day) and all others. Each day is defined by N variables (i.e., the rain gauge records, with N = 86 for 2000–01 and 96 for 2004–05). Ward’s clustering method (Ward 1963) is then used to regroup the days by minimizing the intraclass inertia (i.e., the heterogeneity between patterns of a given class). The most similar individual events, noted A and B, are first identified using the minimum distance in Ψ, and grouped together to make a first cluster R. The distance d(R, Q) between R and any other cluster Q (individual event or group of events), constituted by mR and mQ events, respectively, is a linear function of the distances of Q with respect to the original clusters A and B merged into R. It is defined as
where mA and mB are the number of events of A and B. Distances d(A, Q), d(B, Q), and d(A, B) are obtained directly from the distance matrix Ψ (or, at further stages of the classification procedure, from distances computed at earlier steps).
a. Seasonal rainfall
Figure 3a shows the daily area-averaged rainfall time series of the observations and WRF for 2000–01 and 2004–05. The observed seasonal rainfall is larger in 2004–05 (9.8 mm day−1) than 2000–01 (6.2 mm day−1), in line with a higher (lower) number of rainy (dry) days. Four isolated maximum values exceeding 80 mm day−1 in 2004–05 can be identified in Fig. 3a, versus only one in 2000–01. However, this rainfall peak reaches a much higher value (191.8 mm day−1 on 6 January 2001) than its counterpart for 2004–05 (113.7 mm day−1 on 23 March 2005). Its timing is well reproduced by WRF, although associated amounts are underestimated by ~30%. This event will attract a case study, further documented in section 3c below.
The spatial distribution of observed seasonal rainfall amounts is similar for the two years (Figs. 3b,c), with generally more rainfall on the east than on the west. Except for very few stations, 2004–05 is wetter than 2000–01, especially over the northern to southwestern parts of the island. Cumulative rainfall averaged over all stations is 1766 (1116) mm for 2004–05 (2000–01).
Figure 4 compares the cumulative summertime (November–April) rainfall over the two periods produced with 680-m inner WRF nest simulations (Figs. 4a,c) and RK interpolation (Figs. 4b,d). From Figs. 5a,b showing the differences between nearest RK grid points and observations, one can appreciate the good performance of the interpolation method used for the two periods. Because RK is constrained to fall within the range of observed values, only small local differences between the two are found, resulting in weak negative biases with low values of RMSE. Figures 5c,d display the differences between WRF and RK for the two periods, thus avoiding the point-to-grid comparison issue.
Overall, WRF reproduces reasonably well the observed west–east rainfall gradient (Fig. 4). Rainfall is sensibly lower on the western part of the island and is highest at midlevel altitude in the eastern areas. In the western areas, rainfall is low (<1 m) at the coast and increases with altitude (up to 2 m at the highest levels in 2004–05). In the northeast, rainfall increases with altitude up to about 3 (4) m between 1000 and 1500 m in 2000–01 (2004–05). In RK, high rainfall is also found on the eastern side of the Piton de La Fournaise (location in Fig. 1a), reaching roughly 4 (6) m at 1000 m in 2000–01 (2004–05). WRF tends to produce even larger rainfall amounts farther southeast, with values greater than 5 (10) m over La Fournaise massif in the altitude range 1000–1500 m for 2000–01 (2004–05); the maximum is around 9 (12) m in 2000–01 (2004–05).
On the whole, WRF produces slightly wetter conditions than RK (Fig. 5), as indicated by a positive mean bias error (MBE) of 0.073 and 0.269 m in 2000–01 and 2004–05, respectively (Figs. 5c,d). Spatially, the differences are more pronounced. WRF tends to produce lower rainfall amounts in the northeastern areas from low- to midlevel altitudes (roughly below 1000 m). The amplitude of the differences between WRF and RK is quite large there, reaching up to −2 m at midlevel altitude in the north for 2004–05. In sharp contrast, WRF produces much larger rainfall amounts than RK on the eastern slopes of the Piton de La Fournaise volcano, reaching a maximum positive difference of 6.9 m for 2000–01 and 8.1 m for 2004–05 between 1000 and 1500 m. WRF also produces larger amounts than RK in the western highs. In particular, large positive differences (1–2 m) are found between 1000–1500 m from the northwest to the south for 2000–01, and in many places above 1500 m at the west and in the center of the island for 2004–05. The fact that WRF is drier in the flatter areas of the island and too wet elsewhere may be related to the orographic forcing features possibly too strong in the model. For 2004–05, fine contrasts with localized small patches of higher rainfall are also observed along the edges of the depressed areas (Cirques) inside the island with the RK method (Fig. 4d), producing there noisy patterns of WRF minus RK differences (Fig. 5d). This result, hypothesized to be an artefact of the RK method, is also consistent with the stronger importance of the “slope” predictor for 2004–05 than for 2000–01 [sections 2d and 2e; Fig. 2; Eqs. (4)–(5)].
Although the rain gauge network over La Réunion is dense enough for RK to produce reasonable rainfall estimates, it is, as all in situ measurements, far from perfect because it tends to undersample steep, elevated, and inaccessible sectors in the hinterland parts of the island. This is especially the case over the eastern slopes of the volcano and along the edges of the Cirques, where large differences between WRF and RK have been reported (Figs. 5c,d). Because of the lack of observations at these locations, it is difficult then to assess which method achieves better performance, although the presence of localized rainfall maximums along the steepest sections of the slopes in RK for 2004–05 (Fig. 4d), for example, seems to be questionable physically.
Figure 6a–c shows the seasonal rainfall amount differences between the two years as given by the observations, RK and WRF. Overall, 2004–05 is wetter than 2000–01 everywhere except for five stations in the southeast: two at the coast and three at high-level altitude (Plaine des Cafres, location in Fig. 1a). The largest differences (exceeding 1 m) are concentrated in the northeastern half of the island. These seasonal contrasts can be explained by the mean moisture fluxes over the SWIO region in the lower layers of the troposphere (Figs. 6d–f), which present a marked southerly (northerly) component in 2000–01 (2004–05) inhibiting (favoring) moisture advections toward La Réunion. Differences between the two years additionally confirm that moisture convergence is larger in 2004–05 in the neighboring of the island (Fig. 6f).
By construction, rainfall differences in RK are close to observed ones (Fig. 6b), although the area of negative differences gets artificially larger, and even reaches the northwestern tip of the island, without any counterpart in the rain gauge records. WRF fails at reproducing the observed pattern, especially in the western areas, where unrealistic negative values are found. In addition, marked differences also prevail in the southeast, although this part of the island is poorly documented. This result is rather surprising since large-scale thermodynamics and dynamics is relaxed toward ERA-Interim solutions. They could result from recurrent errors at the daily or intraseasonal time scales impacting the overall seasonal amounts.
b. Intraseasonal rainfall variability
In this section we refine analyses at the intraseasonal time scale. Instead of considering full time series constituted of 181 days for each season of interest, we adopt a more compact approach derived from an objective classification technique (section 2e) that regroups together the most similar daily rainfall fields into a limited number of statistically coherent classes. The methodology was applied separately for both years, but we present here the results for 2004–05 only (results for 2000–01 are very similar, except that intraseasonal variability is less marked and the number of dry days is larger). The classification tree (Fig. 7a) shows the successive regroupings performed by the algorithm for 2004–05. It is here truncated at six classes because this appears as a good compromise between compactness and intracluster heterogeneity: since no statistical test exists to determine the optimal number of classes in this methodology, this choice mostly depends on the level of details expected from the clustering (Pohl and Camberlin 2014). Figure 7b shows that, at the scale of the whole island, the six clusters strongly separate dry and wet days (ascribed to clusters 1–2 and 5–6, respectively). The “dry” clusters clearly concentrate the largest number of days: 126 days out of 181 are regrouped into clusters 1–2.
The spatial patterns associated with the six clusters are shown in Fig. 8. Corresponding moisture flux composites are shown in Fig. 9: they are represented over WRF’s largest domain in order to document how large-scale circulations are likely to modulate the spatial repartition of rainfall amounts over La Réunion.
Cluster 1 is uniformly dry (1.2 mm day−1 on average: Fig. 8), confirming Fig. 7b. Clusters 2–3 are wet (dry) in the northeastern (southwestern) part of the island, while the reverse prevails for cluster 4. Cluster 6 shows a northward rainfall gradient. All rain gauge stations received large rainfall amounts during the occurrences of cluster 5 (ranking between 15 and 260 mm day−1, with an average of nearly 74 mm day−1 and a median of 50 mm day−1).
Moisture fluxes are physically consistent with these rainfall patterns (Fig. 9). The relative dry conditions found for clusters 1–2 occur concomitantly with southeasterly or easterly moisture fluxes associated with no clear convergence over the neighboring of La Réunion. Cluster 3 is associated with stronger easterly fluxes and larger moisture convergence, in line with the larger rainfall amounts over the eastern part of the island. During occurrences of cluster 4 (6), a strong cyclonic vortex is located east (west) of La Réunion, favoring southerly (northerly) fluxes likely to explain once again the spatial distribution of local rainfall. Cluster 5 is associated with strong northeasterly fluxes that convey moisture and converge over the island, explaining the heavy rainfall recorded there.
Figure 10 shows the composite rainfall fields obtained by the RK methodology (Fig. 10a) and simulated by WRF (Fig. 10b). For WRF, they are computed as the average of daily fields ascribed to each cluster. For RK, they are produced based on the interpolation of mean observed rainfall in each specific cluster. Left-hand panels in Fig. 11 present the differences between RK estimates for the nearest grid points and rain gauges. As for the seasonal time scale, differences between nearest RK grid points and observations are still small (relative MBE < 2%–3% for all clusters) and noisy spatially, confirming again the overall good performance of the interpolation method used. Right-hand panels in Fig. 11 show the differences and between WRF and RK methods.
For RK (Fig. 10a), all clusters except cluster 1 show to some extent the localized rainfall maxima along the slopes of the inner depressions (Cirques), already discernible at the seasonal time scale (Fig. 4d) and hypothesized to be an artifact of the method. Such rainfall distribution is not simulated by WRF (Fig. 10b), resulting in marked differences between the two approaches locally (Fig. 11, right-hand panels). Figure 10a illustrates that these patterns are also prevalent at the daily and intraseasonal time scales and exhibit some recurrence. These results differ from the drier 2000–01 year, for which such patterns are not obtained at the seasonal time scale (Fig. 4b) and only concern one cluster formed by 6 days (not shown). Similarly, WRF produces nonnegligible rainfall amounts along the eastern slopes of La Fournaise volcano during all clusters, including the generally dry cluster 1. The positive differences between WRF and RK there, already found at the seasonal scale (Figs. 5c,d), are also valid at the daily–intraseasonal time scale (Fig. 11, right-hand panels). Of course it is not possible to state which approach produces the best results because of the lack of meteorological stations in these sectors of the island.
Temporally, WRF may give the impression to smooth (i.e., to underestimate, the magnitude of day-to-day variability), because it is clearly associated with generalized wet (dry) biases during the occurrences of dry (wet) clusters (Fig. 11: see dry clusters 1–2 vs wet clusters 5–6). This statement is also verified as far as the spatial distribution of rainfall is concerned: for almost all clusters, WRF under- (over-) estimates rainfall amounts where observed amounts are largest (lowest). Several reasons may be invoked to explain this moderate performance of WRF. (i) At finer scales its internal variability increases (Crétat et al. 2011), due to the chaotic behavior of the regional atmosphere inside domains 2–4 (domain 1 being nudged toward the reanalyses). This could result in weakly reproducible rainfall peaks due to enhanced uncertainties in the location and propagation of rain-bearing systems and/or convective cells, the inner WRF domains (2–4) being not relaxed toward reanalyses. (ii) The methodology used to project the observed clusters (Fig. 8) onto WRF-simulated rainfall is based on co-occurrences, which do not take into account the possible timing errors produced by the model (Fig. 3a). Another methodology, based on the resemblance between observed and simulated rainfall fields, could have led to better results, albeit an imperfect timing.
c. Extreme events
La Réunion is known to have world records in terms of rainfall amount for time scales comprised between 12 h and 15 days (Barcelo et al. 1997). They are mostly due to the influence of TCs over the SWIO region that recurrently affect the local climate and may contribute, within a few hours or days, to sizable fractions of the overall seasonal rainfall. During the two years of interest here, one can clearly identify an extreme event that occurred on 6 January 2001 (Fig. 3a), which was associated with an average of nearly 190 mm within 24 h over the island. This rainfall peak was due to TC Ando, which was located roughly 215 km northwest of the island, between La Réunion and the eastern coast of Madagascar (Fig. 12a). Ando was the first TC in the SWIO sector for 2000–01. It was associated with wind speeds of ~220 km h−1 (62 m s−1) at the surface, with gusts up to 75 m s−1 (corresponding to a category-5 TC according to the Saffir–Simpson scale), and a central sea level pressure (SLP) of 922 hPa. WRF simulation (Fig. 12b) shows wind speeds of ~140 km h−1 (38 m s−1) at the most, and a central SLP (~953 hPa) at the least, on 7 January. This was obtained from the model outputs over domain 1 only. Given that data are archived every day, the maximum values reached effectively each day may be lacking. In addition, the intermediate resolution over domain 1 (~43 km) may give smoother results than expected. On 6 January, the vortex associated to Ando caused strong northerly fluxes and moisture convergence over La Réunion (Fig. 12c), explaining the heavy rainfall recorded over the island. On the same day, one can also note the presence of Tropical Cyclone Bindu, located further northeast over the tropical SWIO, which was sensibly weaker than Ando and migrated nearly 500 km east of the Mascarene Islands without causing significant rainfall there.
Rain-gauge records over La Réunion confirm that the rainfall peak associated with Ando was mostly concentrated on 6 January (Fig. 13a), although weak-to-moderate daily amounts were also recorded on 5 and 7 January inland. This corresponds indeed to the maximum strength of the cyclone, and also its location closest to the island. Daily rainfall records reach up to 1200 mm in elevated sectors, maxima corresponding to the slopes of both volcanoes. Of course, corresponding RK fields give a very similar picture of these 3 days. As already suggested by Fig. 3a, WRF does a very good job, both in terms of timing and amounts. Figure 13b shows that the spatial pattern of daily rainfall amounts is also rather realistic. The main differences with RK (Fig. 13b) consist in weaker rainfall amounts over the whole island, except on the eastern side of La Fournaise volcano and the two Cirques of Mafate and Salazie (see Fig. 1) near the Piton des Neiges, in the center of the island. Roughly similar patterns, albeit weaker in magnitude, are also found for 5 and 7 January.
This good performance of WRF can be partly attributed to the relaxation of some of its prognostic variables performed over its largest domain (section 2b), which ensures that transient variability in WRF is phased with the observations. This causes the TC to be at the correct location and at the correct time (Figs. 12b,c). Yet, localized effects of Ando over La Réunion are calculated in domain 4, which is not relaxed toward the reanalyses, hereby suggesting that dynamical downscaling methodologies can be successfully used to regionalize the rainfall field at high resolutions during extreme events, such as those associated with TCs.
4. Discussion: RK-upscaled versus WRF-downscaled altitude–rainfall gradients
At the fine resolutions such as those used for this work, RCM evaluations against observations require an observational network as representative as possible of the region of interest. This is even more important in the case of a terrain with a complex topography, such as La Réunion. Of course, the “perfect” network does not exist: most of the time in situ measurements over- (under-) sample low, accessible and rather flat (elevated, inaccessible, and steep) sectors. In the results presented above, two singular areas emerge, where the network is less dense and where RK is hypothesized to produce potentially spurious rainfall estimates: the sharp, almost vertical edges of the three inner depressions (Cirques, Fig. 1a), and the volcano site (La Fournaise). Cross-validation techniques allow us quantifying RK biases in the location of the stations deliberately removed from the statistical model, but the question of the rainfall amounts where no station exists remains nonetheless unresolved.
Along the eastern slopes of the volcano, WRF produces much more rainfall than RK over all time scales. According to the model, this is the wettest part of the island, a result intuitively consistent with its location facing the dominant moisture fluxes and its elevation causing orographic uplift. Lava often flows there during volcanic eruptions, about once a year, explaining the lack of a permanent weather station. This makes it difficult to assess whether WRF produces realistic rainfall amounts. Similar issues, albeit of lesser importance, concern the Cirques.
To further document the effects of the weak rain gauge density in the Fournaise sector, the scatterplots in Fig. 14b illustrate the relationship between elevation and seasonal rainfall amounts as derived from observations, RK and WRF, over the eastern volcano site (see Fig. 1a). Obvious differences between WRF and RK are found in the elevation–rainfall relationship. According to WRF, a positive association exists between altitude and seasonal rainfall up to 1500 m, above which rainfall declines. This rainfall maximum (10 and 12 m in 2000–01 and 2004–05, respectively) lies well above the highest rain gauge there (820 m, Hauts de Sainte-Rose station). In contrast, and by construction, the heaviest rainfall totals estimated by RK are close to the highest rainfall value observed at any rain gauge station over the island: in 2004–05, the observed maximum value (~6 m) is indeed obtained over the volcano at the Hauts de Sainte-Rose station; in 2000–01, this maximum (~4 m) is rather obtained in the east at a lower altitude (660 m, Takamaka station). The only possible evaluation of the model quality these is through indirect comparisons with Barcelo and Coudray (1996), who completed the Météo-France rain gauge network over La Fournaise volcano with 12 additional but temporary rain gauges spanning altitudes from 220 to 2490 m during one single year (1993). They found a maximum rainfall zone over 12 m yr−1 on the windward eastern slopes of the massif between 1300 and 1800 m, and locally up to 2000 m. Therefore, their results suggest that RK strongly underestimates rainfall amounts, and produces a strongly biased elevation–rainfall relationship, at the volcano site. The statistical relationships between rainfall and the 10 predictors identified at the scale of the whole island are not appropriate to the microscale on the eastern slopes of the volcano. In this respect WRF seems to produce more coherent and probably more realistic results locally.
Physically, the presence of this maximum is directly related to the presence, altitude, and strength of the thermal inversion layer. When it is strong and low, rainfall is restricted to low elevations. During the warm (wet) season, when the thermal layer inversion is at a higher altitude or even absent (mainly when tropical storms pass close to the island), heavy rains occur near the summit. This is for instance the case for 6 January 2001, when a different elevation–rainfall relation is obtained (Fig. 14a) at the time TC Ando was in close proximity to La Réunion (section 3c). As noted previously, WRF produces similar patterns to RK with large rainfall amounts over the summits of the island and no maximum at midlevels. Over the whole island and under the influence of a strong TC, the elevation–rainfall relationship is also rather similar between RK and WRF (Fig. 14a), although RK tends to produce heavier rainfall between 0 and 2000 m.
This work aims at exploring to what extent a nonhydrostatic regional climate model (RCM) can regionalize rainfall over a complex terrain at very high resolutions (680 m), for various time scales including year-to-year differences between contrasted rainy seasons, seasonal means, intraseasonal variability, and daily extreme events. The region of interest is La Réunion island, located in the southwest Indian Ocean. The period analyzed includes the austral summer wet season (November–April) for 2000–01 (abnormally dry) and 2004–05 (abnormally wet).
To better sample topography (rain gauge stations being often located in weakly elevated, relatively flat and accessible areas) and avoid the grid-to-point comparison issue, RCM outputs are compared to regression-kriging (RK) products. RK consists in establishing statistical relationships between a combination of predictors and observed rainfall amounts using a multiple linear regression, and then kriging the residuals. The originality of this work is the introduction of atmospheric predictors, in combination with the more traditional morphometric ones.
At seasonal time scales WRF produces rainfall patterns rather similar to those observed. At intraseasonal and daily time scales, however, WRF biases are generally larger, due to some timing errors (rainy events shifted by a few hours, leading to wrong daily amounts concerning two contiguous days) and spatial errors (denoting an increased difficulty to simulate the correct amount of rainfall at the correct location at fine temporal and spatial scales).
The capability of WRF to regionalize an extreme event is assessed though the case study of TC Ando, which caused heavy rainfall over La Réunion (up to 1200 mm day−1 in the upper parts of the island) on 6 January 2001. Results show that WRF can be successfully used to regionalize the rainfall field during extreme events. This suggests that both dynamical downscaling and geostatistical interpolation are useful to fill the gap between the large-scale atmospheric configurations associated with the synoptic situation in the region, and the local scales that more directly concern the local populations. This promising result encourages us to pursue this work and consider a larger number of TCs and evaluate how their effects on La Réunion are modulated, on the one hand by their location and intensity, and on the other hand, by local surface conditions.
This work is a contribution to the LEFE/IDAO VOASSI program funded by CNRS. WRF was provided by the University Corporation for Atmospheric Research website (http://www.mmm.ucar.edu/wrf/users/download/get_source.html). ERA-Interim data were provided by the ECMWF Meteorological Archival and Retrieval System (MARS). Calculations were performed using HPC resources from DSI-CCUB, Université de Bourgogne.