Statistical Downscaling of ERA-Interim Forecast Precipitation Data in Complex Terrain Using LASSO Algorithm

1 College of Geographical Sciences, Fujian Normal University, Fuzhou 350007, China 2 Fujian Provincial Engineering Research Center for Monitoring and Assessing Terrestrial Disasters, Fuzhou 350007, China 3Department of Geography, Ludwig-Maximilians-Universität München, 80333 Munich, Germany 4 Institute of Water Management, Hydrology and Hydraulic Engineering, University of Natural Resources and Life Sciences, 1190 Vienna, Austria


Introduction
Precipitation is one of the most important variables for a large variety of environmental processes and its spatial and temporal variations directly influence the local as well as the global water, energy, and matter cycle (e.g., [1][2][3][4][5][6]).Precipitation is also an essential input parameter for land surface models in fields such as in hydrology, ecology, and climatology [7][8][9].Traditionally, precipitation is measured by rain gauges.However, most of them are located in homogeneous terrains and in valley regions.In complex terrains, few gauge stations exist due to difficulties in snow depth measuring and maintenance of stations which result in a lack of long-term and high resolution records [10].In order to overcome these limitations, in the past decades, GCMs have been widely applied to meet the specific needs of environmental impact models by providing time series of precipitation and plausible scenarios of change [8,11].
However, it is well known that the coarse spatial resolution (∼300 km) of GCMs limits the reliable use of these data in decision making and model based impact studies [12][13][14][15].Specifically, GCMs only provide precipitation characteristics that are based on the mean elevation of the pixel, thus not considering subgrid variability of topography and consequent atmospheric features.Local processes, such as orographic precipitation along mountain slopes, are not represented in the coarse grids [12].Even if RCMs nested in GCMs provide better descriptions of local scale characteristics, based on the finer spatial resolution of 10-50 km, they do not fulfill the requirement of hydrological and climatic impact models, which typically run on the scale of 0.1-1 km [11,[16][17][18][19][20][21].Thus, correction and downscaling are necessary for both GCM and RCM outputs before they are applied for environment impact models [22][23][24][25][26].
In a broad sense, downscaling includes two purposes: one is obtaining future emission scenarios (time series) for stations or regions or catchments from global climate models (e.g., HadAM3P) using the established statistical relationship between measurements and large-scale predictors.The other purpose is obtaining time series for nonmeasurement area.This also could be treated as spatial interpolation or disaggregation to some extent.Many previous studies have shown that a reproduction of station data and the generation of future times series at the respective locations is extremely useful (e.g., [19,22,24]).
In terms of downscaling methods, dynamical and statistical downscaling are the two main approaches.The latter approach establishes the statistical connections between large-scale circulation variables (predictors) and local observed variables (predictands) [27,28].Compared to dynamical downscaling, statistical downscaling methods have lower computational demands and allow for a fast application [27,29].Maraun et al. [11] comprehensively reviewed precipitation downscaling methods from an end user's point-of-view.For example, local intensity scaling (LOCI) and quantile-mapping (QM) are the standard methods used for correcting GCM or RCM outputs with respect to local observations.Other methods focus on investigating the relationship between local precipitation and large-scale atmosphere circulations, which vary from linear regression (e.g., [30][31][32][33][34][35]) to complex nonlinear models (e.g., [36][37][38][39][40]).Although numerous studies were carried out, a general standardized precipitation downscaling method still does not exist, especially for complex terrains.Furthermore, a separate predictor selection process (e.g., principal components analysis) is usually implemented in order to search for the most sensitive variables with regard to precipitation variations (e.g., [34,39,[41][42][43][44]).However, this procedure costs the additional computation time.Therefore, it is of particular interest to introduce new approaches, especially for the stations where benchmark methods do not work at all.
To this end, a new machine learning method, the "least absolute shrinkage and selection operator (LASSO)" algorithm, is introduced for downscaling ERA-Interim forecast precipitation data in complex terrain.Compared to standard linear downscaling approaches, LASSO is also well suited for possibly underdetermined linear regression problems, as well as for joint estimation and continuous variable selection.It is tested and validated against three different methods, local intensity scaling (LOCI), quantile-mapping (QM), and stepwise regression (Stepwise), using data from 50 meteorological stations located in the high mountainous region of the central Alps (Figure 1).This paper is structured as follows.Section 2 describes the ERA-Interim forecast precipitation data and meteorological observations in the study area.Section 3 describes the four downscaling methods as well as the evaluation criteria.The downscaling results and methods comparison are presented in Section 4, while finally a discussion and subsequent conclusions are given in Section 5.

Datasets
2.1.ERA-Interim.We make use of the ERA-Interim forecast precipitation data provided by the European Centre for Medium Range Weather Forecast (ECMWF) for the year 1979 onwards and continuing in real time [45,46].ERA-Interim shows some improvements when compared with ERA-40 in these aspects: a representation of the hydrological cycle, an improved description of the stratospheric circulation, and an enhanced handling of biases [46][47][48][49].Cycle 31r2 of ECMWF's Integrated Forecast System (IFS) was used in here.The model in this configuration comprises 60 vertical levels, with the top level at 0.1 hPa; it uses the T255 spectral harmonic representation for the basic dynamical fields and a reduced Gaussian grid (N128) with an approximately uniform spacing of 79 km [46,50].ERA-Interim assimilates four analyses per day at 00, 06, 12, and 18 UTC.Furthermore, two 10-day forecasts with a 3-hour resolution are initialized based on the 00:00 UTC and 12:00 UTC analyses.

Test Sites.
Daily total precipitation of the period 1983-2010 at 50 meteorological stations was made available through the interactive tools of IDAWEB, which is designed by MeteoSwiss (the Swiss Federal Office of Meteorology and Climatology) providing free, available, and extensive archive data of ground level monitoring networks.Table 2 lists the information about stations and Figure 1 shows the locations of test sites.The stations are located within a large range of altitude from 381 m to 3305 m.Among these stations, 16 are located below 500 m and 11 are situated between 500 m and 1000 m, 9 between 1000 m and 1500 m, 8 between 1500 m and 2000 m, and 6 above 2000 m.The observations and ERA-Interim data are processed for the same period.The available data is partitioned into two periods, 1983-1999 for calibration and 1999-2010 for validation.A 1 mm threshold was defined for defining a dry/wet day.It is necessary to note that the data of the stations GUE, PAY, and GVE are used within the ERA-interim data assimilation procedure, given their status as WMO SYNOP stations [46,51].According to the information of the ECMWF, it can be assumed that the majority of the stations (47 of 50 sites) are not used by ERA-Interim and therefore represent a fully independent dataset.

Local Intensity Scaling (LOCI).
LOCI is a robust method to directly correct GCM or RCM outputs for local observations.Although GCMs or RCMs are partly unrealistic due to their coarse resolution, they contain valuable information about the actual precipitation [11].
The assumption is realized by a so-called scaling factor, calculated from observation and climate model data of a reference period, which is then expanded to scenarios data.Here, LOCI is applied as the benchmark method for comparison with LASSO.LOCI was developed by Widmann and Bretherton [52].Widmann et al. [53] used it for scenario precipitation corrections.Not only GCMs, but also RCMs were corrected using the LOCI approach [54][55][56].Schmidli et al. [57] further modified LOCI for precipitation occurrence and amount correction, separately.In this study, LOCI is implemented based on a monthly scaling factor which is calculated in three steps as follows: where   is the target station precipitation,  Val ERA the undownscaled ERA-Interim data for validation,  thres ERA ERA-Interim precipitation threshold, and  thres obs the observation threshold and the brackets present the frequency condition judgment function.Here, 1 mm is used to define wet/dry days and SF is the scaling factor.In the first step, an adjusted threshold for ERA-Interim data is found that matches the occurrence of wet/dry days, based on the 1 mm threshold of observation.In a second step, the scaling factor is obtained and then, finally, the target station precipitation is calculated.

Quantile-Mapping (QM)
. QM introduced by Panofsky and Brier [58] is a popular statistical transformations approach to correct GCM and RCM outputs straightforwardly [25,[59][60][61][62][63][64][65].The distribution function (e.g., cumulative distribution function, CDF) of model precipitation is first adjusted to match the distribution of observations.Subsequently, this matched distribution is used for unbiased model (or future scenario) data.The mapping is usually implemented based on empirical quantiles or quantiles of gamma distributions [11,64].In this study, the corrected ERA-Interim can be obtained via where  −1 obs,cal and  −1 ERA,cal is the inverse CDF of observations and ERA-Interim for calibration, respectively, and  ERA,cal is the CDF of  Val ERA .

Stepwise Regression (Stepwise).
Stepwise Regression (Stepwise hereafter) is an automatic procedure where Stepwise combines an ordinary regression (3) with a predictor variable selection procedure.Three main approaches are used in Stepwise according to the relevant selection sequence: forward selection, backward elimination, or bidirectional elimination.The advantage of stepwise regression is easily explained and implemented.Several previous studies have used stepwise regression for different purposes.For example, Harpham and Wilby [39], Hessami et al. [42], and Huth [66] adopted Stepwise for predictor selection; Agnihotri and Mohapatra [30] applied it to occurrence estimation of daily summer monsoon precipitation.In this study, stepwise regression is adopted to test LASSO.Stepwise regression is implemented using backward elimination method with a significance level of 0.05, where  is the ×1 response vector,  is the × variable vector,  is the ×1 parameter vector, and  is the random errors.We used the same set of 20 variables (Table 1) for Stepwise and LASSO, for a better method comparison.For precipitation occurrence,  is defined as 1 for wet days (>1 mm) and 0 for dry days (<1 mm).All variables are standardized to make data fall between 0 and 1.

LASSO Algorithm.
Least absolute shrinkage and selection operator (LASSO) is an alternative regularized version of least squares, which is useful for feature selection and to avoid overfitting problems.LASSO shrinks the estimates of the regression coefficients towards zero to prevent overfitting problem and to reduce variables by using a penalty parameter [67,68].To simplify understanding, the history of LASSO is introduced briefly.The following equation presents the ordinary least squares regression (OLS) that tries to minimize the error RSS (Root of Sum of Squares), OLS is not always satisfactory for minimizing the RSS, especially when  contains a large number of variables.A penalty parameter  was added based on the normal OLS in LASSO (see the following equation): where V is the number of variables.LASSO imposes intentionally that some coefficients have to be zero, thus achieving a sparse model.Thus, the penalty parameter (regularizer)  controls the level of sparsity of the resulting model.In this study, we applied an efficient algorithm for solving LASSO [69].Also, we defined the value 1 for wet days (>1 mm) and 0 for dry days (<1 mm), and all variables are standardized the same as Stepwise.

Evaluation Criteria.
Precipitation downscaling procedure is implemented in two steps.Firstly, precipitation occurrence is modeled by the four methods, respectively.Secondly, precipitation amount conditional on the occurrence of wet days is modeled subsequently.Please note that the results of precipitation occurrence ranges from 0 to 1.We defined 0.5 as the threshold value to classify dry/wet days.The modeled precipitation amount also could be negative values.Therefore, we set these negative values to zero.The root mean square error (RMSE) and the mean absolute error (MAE) are used for the assessment of precipitation amount, and correspondence ratio (CR) is applied for the evaluation of dry/wet days classification accuracy ( 6)-( 7), where   and   are observed and modeled precipitation amount on wet days, respectively. is the number of wet days, where  dry and  wet are the numbers of dry or wet days that are correctly classified by downscaling methods.

Validation of the Original ERA-Interim Forecast
Precipitation Data.Table 3 shows the comparison of ERA-Interim daily precipitation forecasts with observations at 50 meteorological stations from 1983-2010.CR as well as the RMSE and MAE in mm is listed.CR varies from 0.68 to 0.86 and RMSE changes in the range of 4.40-11.13mm, while MAE ranges from 2.22 to 5.19 mm.The large errors show that there is a great need for the correction and downscaling of ERA-Interim data at local stations.The altitude differences change sharply in a large interval.We use the abbreviations for the characterization of the stations.Station COV is 1687 m higher than ERA-Interim grid height while station SCU is 515 m lower than grid height.ERA-Interim shows good agreement with the occurrence of observations (0.68 to 0.86), but large deviations with respect to the amount of precipitation on wet days.However, in contrast to temperature estimations, the difference in altitude is unable to explain the observed bias.For example, station DOL has the highest CR with the value of 0.86, but with a large elevation difference of 971 m.Station ZER has a CR of 0.74 but its elevation matches the grid height very well with an altitude gap of only 86 m.This again suggests that the general relationship between precipitation and elevation is not easy to obtain, due to the great variability in the interaction between atmospheric circulation and complex topographical characteristics.
Figure 2 illustrates the bias (in percent) between ERA-Interim data and station precipitation observations for the time period 1983-2010.Positive values indicate that annual precipitation is overestimated by ERA-Interim and vice versa.In general, ERA-Interim data overestimate observations the observed annual precipitation for the majority of test sites; only 7 of 50 sites are underestimated.Among them, station GSB has the largest negative bias with a value of −44.2%. 5 out of 50 sites show an overestimation of more than 100%.Station WFJ and DAV are located in the same ERA-Interim grid.However, ERA-Interim overestimated station WFJ with 7.5% and station DAV with a value of 45.6%.Besides, stations PIL and LUZ are also in the same ERA-Interim grid.ERA-Interim underestimates PIL with 25.2%.In the contrast, ERA-Interim overestimates LUZ with a value of 28.6%.The lower stations that are mainly located at the northern part of the Alps are underestimated by ERA-Interim.

Evaluation of Downscaling Methods in Precipitation
Occurrence.The performance of the four downscaling methods, as well as the ERA-Interim original data in precipitation occurrence in the validation period 1999-2010, is summarized in Table 4.The averaged CR showed that the four downscaling methods were much better than the original ERA-Interim data.QM and LASSO had the same range of CR from 0.82 to 0.88.LOCI ranged from 0.81 to 0.88 and Stepwise varied from 0.82 to 0.87.QM estimated a worse CR than the original ERA data at station GSB.LASSO predicted the worst occurrence at station DOL.The four methods had the same performance at 7 sites.For 11 sites, LASSO and Stepwise were the best methods in precipitation occurrence modeling.Figure 3 illustrates the averaged correspondence ratio for each month in the validation period 1999-2010.It reveals that LOCI, QM, Stepwise, and LASSO captured the general tendency of the monthly precipitation occurrence for the whole year.LASSO and Stepwise identified better dry and wet days than LOCI and QM in May, but they performed worse in October.LASSO classified dry/wet days similarly with Stepwise for all 12 months.The result shows that LASSO is quite suitable for dry/wet days classification.
Figure 4 shows the used variables in Stepwise and LASSO for dry/wet days classification.LASSO selected more variables than Stepwise.For LASSO, P ERA, T 700, and TCW were the most frequent (50 sites) variables.RH 500 was the least frequently (17 sites) applied by LASSO.P ERA and U 10 were the most frequent (more than 45 sites) variables in Stepwise.T 500 was the least frequent variable.

Evaluation of Downscaling Methods in Precipitation
Amount.The overall performance of four downscaling methods, as well as the ERA-Interim original data in precipitation amount conditional on the occurrence of a wet day for the validation period 1999-2010, is summarized in Table 5 Taking station WFJ as an example, Figures 5 and 6 illustrate the downscaled daily precipitation for the four downscaling methods, as well as the original ERA-Interim data in January and July, respectively, for the validation period 1999-2010.In general, four downscaling methods capture the tendency of the daily precipitation in January and July.However, LOCI and QM are generally inclined to misidentify the heavy precipitation events, in particular in January 2008 and July 2003.Stepwise and LASSO underestimated heavy precipitation events compared with LOCI and QM.
Figure 7 illustrates the comparison of observations with downscaled annual precipitation for the four methods, as well as the original ERA-Interim data in the validation period 1999-2010.LASSO and Stepwise predicted the annual precipitation best.In general, LOCI and QM modeled drier results, with mean value of −7.7% and −5.0%, respectively.31 out of 50 sites were underestimated by LASSO with the range from −19% (station KLO) to −0.1% (station RUE) while 19 sites were overestimated by 0.2% (station SCU) to 31.6% (station MLS).The majority of sites were underestimated by LOCI (48 out of 50 sites) and QM (40 out of 50 sites) while 30 of 50 sites were underestimated by Stepwise.Figure 8 compares the averaged sum of daily precipitation for each month between observations and the four downscaling methods, as well as the original ERA-Interim data in the validation period 1999-2010.It is easy to find that LASSO and Stepwise predicted more precipitation in the whole year than LOCI and QM with exception of August and November.LASSO reproduced best monthly precipitation for February, March, April, and November.Stepwise was the best method for May, June, July, and October.QM performed best in January, August, and September.LOCI outperformed the other methods in December.Figure 9 shows the used variables in Stepwise and LASSO for precipitation amount prediction on wet days.LASSO again selected more variables than Stepwise.For LASSO, P ERA, H 850, H 700, H 500, T 850, and MSLP were applied at all 50 sites.SH 500 was applied at 25 sites.P ERA was the sole variable used by all test sites by Stepwise.SH 850 and U 10 were the second most frequent (34 sites) variables.RH 700 and RH 500 were the least frequent variables (10 sites) by Stepwise.For the individual sites, station VAD used the least variables (11 variables) while 19 stations selected all 20 variables in LASSO.Station SBE applied the most variables (16 out of 20 vairables) while station WYN only applied 4 variables in Stepwise.

Discussion and Conclusion
The comparison between ERA-Interim and observations showed that ERA-Interim has a large error (5.79 mm of RMSE and 2.84 mm of MAE) in the central Alps (Table 3).Thus, there is a great need for the correction and downscaling of ERA-Interim data.This study compared four downscaling methods, LOCI, QM, Stepwise, and LASSO, for downscaling of ERA-Interim daily precipitation data in the central Alps.
As a frequent input variable for hydrological models, daily precipitation is always corrected or downscaled due to the limits of rain gauges.In the previous studies (e.g., [57,63]), LOCI and QM methods have been widely used for the bias correction with the advantages of maintaining the variation and distribution of historical data.Although LOCI and QM captured the best estimation of daily precipitation occurrence, the reduction of error is not significant by these two methods and even worse than the original ERA-Interim data for 13 out of 50 sites.It demonstrates that straightforward methods are not always suitable for downscaling local observations in complex terrain.
LASSO algorithm simulated the occurrence of daily precipitation as well as LOCI, QM, and Stepwise; it captured the occurrence for all test sites generally.Compared to the other three downscaling methods, it reduced the most amount error to 16.3% of MAE and 12.2% RMSE.LASSO also predicted the best annual and monthly precipitation.Although the LASSO algorithm has been developed for more than 15 years by statistician, the application in geosciences is still at the early-stage (e.g., [70]).A main practical challenge in applying LASSO for precipitation downscaling is that precipitation in heterogeneous terrain is such a complex process which tends to use more variables and to overestimate observations.It has shown that LASSO tends to use more variables than Stepwise.To avoid the overfitting problem, the penalty (or regularization) parameter  plays a key role.The cross-validation and hundreds of runs are necessary to find an appropriate penalty parameter.Therefore, LASSO is a little bit time-consuming than the other three methods.So far, the work presented herein has been limited to the central Alps with 50 meteorological stations providing calibration/validation data sets for testing LASSO algorithm.It would be necessary to extend LASSO method to other high mountainous areas around the world.Besides, this study focuses on the local stations at present; it is of special interest to extend LASSO method for nonstation areas.A limited number of variables (20 variables in this study) derived from ERA-Interim data were applied for LASSO; however, it should also be investigated whether other potential variables such as vapor pressure can be used in the presented approach.
This study focused on the daily total precipitation.Certainly, higher temporal resolution data, such as 3-hourly, would be of great interest in further investigations.Furthermore, applying other reanalysis data sets based on different land surface representations could also be valuable for validating LASSO algorithm.

Figure 2 :Figure 3 :
Figure 2: Annual precipitation bias (bars) between ERA-Interim and MeteoSwiss stations in percent for the period 1983-2010.Positive values indicate an overestimation of the annual precipitation by ERA-Interim and vice versa.The line with dots indicates altitude differences, Δℎ, defined as ERA-Interim grid height minus site elevation.

Figure 5 :
Figure 5: Comparison of observations with downscaled daily precipitation for the four downscaling methods, as well as the original ERA-Interim data in January at station WFJ during the validation period 1999-2010

Figure 6 :
Figure 6: Comparison of observations with downscaled daily precipitation for the four downscaling methods, as well as the original ERA-Interim data in July at station WFJ during the validation period 1999-2010.

Figure 7 :Figure 8 :
Figure 7: Percentage of annual precipitation bias between observations and the four downscaled results, as well as the original ERA-Interim data in the validation period 1999-2010.

Table 1 :
Predictors from ERA-Interim forecast dataset.All variables are aggregated from 3-hourly to daily averages.

Table 2 :
Test sites information.ERA-Interim grid height is also listed.

Table 3 :
Comparison of ERA-Interim forecast daily precipitation with observations at 50 meteorological stations from 1983 to 2010.CR, RMSE, and MAE in mm are listed.The elevations of sites are different with ERA-Interim grid heights, so the altitude differences (Δℎ, ERA-Interim grid height minus elevation) are also labeled.

Table 4 :
Comparison of downscaling methods, as well as the original ERA-Interim forecast data in precipitation occurrence in the validation period 1999-2010.
For the individual sites, station INT selected least variables (12 variables) while station COV, CHA, and CGI chose the most variables (20 variables) in LASSO.Stepwise selected 16 variables at stations GSB, CHA, FAH, and SMA, while it only applied 9 variables at stations SAM and PIO.LASSO applied more variables compared with Stepwise for the majority of stations (49 stations).Although, different variables are applied by the two methods, the performances of dry/wet days classification were similar.
. On the average, Stepwise and LASSO performed similarly in RMSE and MAE, which were much smaller than LOCI and QM.The reduction of error of Stepwise was 12.2% of RMSE and 15.7% of MAE, respectively.LASSO was slightly better than Stepwise in MAE (16.3%) and it had the same error reduction with Stepwise in RMSE (12.2%), whereas LOCI and QM were the worst methods.The reduction of RMSE was only 3.5% and 0.9% for LOCI and QM, respectively.In total, LOCI and QM did not reduce RMSE at 13 stations such as station WFJ.QM was the worst method at 24 stations while LOCI performed the worst in RMSE at stations ENG and CDF.LASSO outperformed the other methods at 9 stations.Stepwise reduced the most RMSE at 40 stations compared with other three methods.However, the four downscaling methods failed to reduce the errors at station CHA.In terms of MAE, LOCI outperformed the other methods at 26 stations.LASSO reduced the most MAE at 22 stations.QM method was not able to reduce MAE at 5 stations.It could be concluded that LOCI and QM method are not always suitable for local stations.It has to be noted that the four downscaling methods did not work for MAE reduction at station DOL, which has a large elevation difference of 971 m against ERA-Interim grid height.

Table 5 :
Comparison of observations with downscaled daily precipitation for the four downscaling methods, as well as the original ERA-Interim data in precipitation amount (mm) on wet days for the validation period 1999-2010.
Selected variables in Stepwise and LASSO for precipitation amount prediction on wet days.