Intercomparison of Machine Learning Methods for Statistical Downscaling: The Case of Daily and Extreme Precipitation

TL;DR

This study compares traditional and machine learning statistical downscaling methods for daily precipitation, finding that simpler linear models like BCSD outperform more complex machine learning approaches in capturing climate features.

Contribution

The paper provides a comprehensive comparison of classical and machine learning methods for climate downscaling, highlighting the limited benefits of advanced ML techniques over traditional models.

Findings

Linear methods outperform non-linear approaches.

State-of-the-art ML methods do not improve over traditional models.

BCSD consistently provides reliable downscaling results.

Abstract

Statistical downscaling of global climate models (GCMs) allows researchers to study local climate change effects decades into the future. A wide range of statistical models have been applied to downscaling GCMs but recent advances in machine learning have not been explored. In this paper, we compare four fundamental statistical methods, Bias Correction Spatial Disaggregation (BCSD), Ordinary Least Squares, Elastic-Net, and Support Vector Machine, with three more advanced machine learning methods, Multi-task Sparse Structure Learning (MSSL), BCSD coupled with MSSL, and Convolutional Neural Networks to downscale daily precipitation in the Northeast United States. Metrics to evaluate of each method's ability to capture daily anomalies, large scale climate shifts, and extremes are analyzed. We find that linear methods, led by BCSD, consistently outperform non-linear approaches. The direct application of state-of-the-art machine learning methods to statistical downscaling does not provide improvements over simpler, longstanding approaches.

Tables

Table 1. Table 2: Large Scale Projection Results: After aggregating daily downscaled estimates to monthly and yearly time scales, RMSE and Skill are computed per location and averaged.

	RMSE (mm)		Skill Score		Correlation
Time-frame	Month	Year	Month	Year	Month	Year
BCSD	31.97	204.78	0.88	0.63	0.85	0.64
PCAOLS	50.01	362.73	0.75	0.27	0.63	0.41
PCASVR	92.17	414.40	0.83	0.69	0.29	0.15
ELNET	46.96	353.67	0.76	0.27	0.71	0.50
MSSL	62.63	597.80	0.56	0.05	0.67	0.40
BCSD-MSSL	31.24	155.04	0.88	0.87	0.82	0.60
CNN	112.21	1,204.27	0.01	0.00	0.59	0.54

Table 2. Table 3: Statistics for ClimDEX Indices: For each model’s downscaled estimate we compute four extreme indices, consecutive wet days (CWD), very heavy wet days (R20), maximum 5 day precipitation (RX5day), and daily intensity index (SDII), for each location. We then compare these indices to those extracted from observations to compute correlation and skill metrics.

	Correlation				Skill Score
Metric	CWD	R20	RX5day	SDII	CWD	R20	RX5day	SDII
model
BCSD	0.43	0.83	0.73	0.70	0.71	0.80	0.84	0.44
PCAOLS	0.25	0.65	0.44	0.67	0.69	0.60	0.65	0.44
PCASVR	0.24	0.81	0.19	0.25	0.78	0.89	0.80	0.65
ELNET	0.36	0.71	0.57	0.64	0.79	0.62	0.63	0.35
MSSL	0.33	0.84	0.56	0.52	0.90	0.63	0.57	0.16
BCSD-MSSL	0.25	0.83	0.70	0.69	0.41	0.75	0.84	0.08
CNN	0.07	—	0.33	-0.30	0.05	0.59	0.01	0.00

Equations

$$\text{E}[Y]=R\ast E[Y|R]\text{ where }R=\{\begin{array}{ll}0,& \text{if}\text{P}(Rainy)<0.5\\ 1,& \text{otherwise}\end{array}$$

$$\hat{\beta }=\arg \underset{\beta }{\min }\parallel Y-X\beta \parallel$$

$$\displaystyle\hat{\beta}=\arg\!\min_{\beta}\big{(}\parallel Y-X\beta\parallel_{2}^{2}+\lambda_{1}\parallel\beta\parallel_{1}+\lambda_{2}\parallel\beta\parallel_{2}^{2}\big{)}$$

$$f(x)=\sum _{i=1}^d{w}_i\times K(x_i,x)+b$$

$$\displaystyle\min_{\boldsymbol{W},\boldsymbol{\Omega}\succ 0}\bigg{\{}\dfrac{1}{2}\sum_{k=1}^{K}\parallel X\boldsymbol{W}_{k}-Y_{k}\parallel_{2}^{2}-\dfrac{K}{2}\text{log}|\boldsymbol{\Omega}|+Tr(\boldsymbol{W}\boldsymbol{\Omega}\boldsymbol{W}^{T})+\lambda\parallel\boldsymbol{\Omega}\parallel_{1}+\gamma\parallel\boldsymbol{W}\parallel_{1}\bigg{\}}$$

$$\boldsymbol{W}^{t+1}|\boldsymbol{\Omega}^{t}=\min_{\boldsymbol{W}}\bigg{\{}\dfrac{1}{2}\sum_{k=1}^{K}\parallel X_{k}\boldsymbol{W}_{k}-Y_{k}\parallel_{2}^{2}+Tr(\boldsymbol{W}\boldsymbol{\Omega}\boldsymbol{W}^{T})+\gamma\parallel\boldsymbol{W}\parallel_{1}\bigg{\}}\\$$

$$\boldsymbol{\Omega}^{t+1}|\boldsymbol{W}^{t+1}=\min_{\boldsymbol{\Omega}}\bigg{\{}Tr(\boldsymbol{W}\boldsymbol{\Omega}\boldsymbol{W}^{T})-\dfrac{K}{2}log|\boldsymbol{\Omega}|+\lambda\parallel\boldsymbol{\Omega}\parallel_{1}\bigg{\}}\\$$

Full text

1 Introduction

The sustainability of infrastructure, ecosystems, and public health depends on a predictable and stable climate. Key infrastructure allowing society to function, including power plants and transportation systems, are built to sustain specific levels of climate extremes and perform optimally in it’s expected climate. Studies have shown that the changing climate has had, and will continue to have, significant impacts on critical infrastructure [13, 31]. Furthermore, climate change is having dramatic negative effects to ecosystems, from aquatic species to forests ecosystems, caused by increases in greenhouse gases and temperatures [43, 32, 20]. Increases in frequency and duration of heat waves, droughts, and flooding is damaging public health [18, 12].

Global Circulation Models (GCMs) are used to understand the effects of the changing climate by simulating known physical processes up to two hundred years into the future. The computational resources required to simulate the global climate on a large scale is enormous, limiting models to coarse spatial and temporal scale projections. While the coarse scale projections are useful in understanding climate change at a global and continental level, regional and local understanding is limited. Most often, the critical systems society depends on exist at the regional and local scale, where projections are most limited. Downscaling techniques are applied to provide climate projections at finer spatial scales, exploiting GCMs to build higher resolution outputs. Statistical and dynamical are the two classes of techniques used for downscaling. The statistical downscaling (SD) approach aims to learn a statistical relationship between coarse scale climate variables (ie. GCMs) and high resolution observations. The other approach, dynamical downscaling, joins the coarse grid GCM projections with known local and regional processes to build Regional Climate Models (RCMs). RCMs are unable to generalize from one region to another as the parameters and physical processes are tuned to specific regions. Though RCMs are useful for hypothesis testing, their lack of generality across regions and extensive computational resources required are strong disadvantages.

1.1 Statistical Downscaling

Statistical downscaling methods are further categorized into three approaches, weather generators, weather typing, and transfer functions. Weather generators are typically used for temporal downscaling, rather than spatial downscaling. Weather typing, also known as the method of analogues, searches for a similar historical coarse resolution climate state that closely represents the current state. Though this method has shown reasonable results [11], in most cases, it is unable to satisfy the non-stationarity assumption in SD. Lastly, transfer functions, or regression methods, are commonly used for SD by learning functional relationships between historical precipitation and climate variables to project high resolution precipitation.

A wide variety of regression methods have been applied to SD, ranging from Bias Correction Techniques to Artificial Neural Networks. Traditional methods include Bias Correction Spatial Disaggregation (BCSD) [46] and Automated Regression Based Downscaling (ASD) [21] and are the most widely used. BSCD assumes that the climate variable being downscaled is well simulated by GCMs, which often is not the case with variables such as precipitation [35]. Rather than relying on projections of the climate variable being downscaled, regression methods can be used to estimate the target variable. For instance, precipitation can be projected using a regression model with variables such as temperature, humidity, and sea level pressure over large spatial grids. High dimensionality of covariates leads to multicollinearity and overfitting in statistical models stemming from a range of climate variables over three dimensional space. ASD improves upon multiple linear regression by selecting covariates implicitly, using covariate selection techniques such as backward stepwise regression and partial correlations. The Least Absolute Shrinkage and Selection Operator (Lasso), a widely used method for high dimensional regression problems through the utilization of a $l_{1}$ penalty term, is analogous to ASD and has shown superior results in SD [42, 19]. Principle component analysis (PCA) is another popular approach to dimensionality reduction in SD [39, 14, 24], decomposing the features into a lower dimensional space to minimize multicollinearity between covariates. PCA is disadvantaged by the inability to infer which covariates are most relevant to the problem, steering many away from the method. Other methods for SD include Support Vector Machines (SVM) [14], Artificial Neural Networks (ANNs) [40, 8], and Bayesian Model Averaging [49].

Many studies have aimed to compare and quantify a subset of the SD models presented above by downscaling averages and/or extremes at a range of temporal scales. For instance, Burger et al. presented an intercomparison on five state-of-the-art methods for downscaling temperature and precipitation at a daily temporal resolution to quantify extreme events [5]. Another recent study by Gutmann et al. presented an intercomparison of methods on daily and monthly aggregated precipitation [17]. These studies present a basis for comparing SD models by downscaling at a daily temporal resolution to estimate higher level climate statistics, such as extreme precipitation and long-term droughts. In this paper we follow this approach to test the applicability of more advanced machine learning models to downscaling.

1.2 Multi-task Learning for Statistical Downscaling

Traditionally, SD has focused on downscaling a locations independently without accounting for clear spatial dependencies in the system. Fortunately, numerous machine learning advances may aid SD in exploiting such dependencies. Many of these advancements focus on an approach known as multi-task learning, aiming to learn multiple tasks simultaneously rather than in isolation. A wide variety of studies have shown that exploiting related tasks through multi-task learning (MTL) greatly outperforms single-task models, from computer vision [48] to biology [27]. Consider the work presented by [10] in which increasing the number of tasks leads to more significant feature selection and lower test error through the inclusion of task relatedness and regularization terms in the objective function. MTL has also displayed the ability to uncover and exploit structure between task relationships [50, 6, 2].

Recently Goncalves et al. presented a novel method, Multi-task Sparse Structure Learning (MSSL),  [16] and applied it to GCM ensembles in South America. MSSL aims to exploit sparsity in both the set of covariates as well as the structure between tasks, such as set of similar predictands, through alternating optimization of weight and precision (inverse covariance) matrices. The results showed significant improvements in test error over Linear Regression and Multi-model Regression with Spatial Smoothing (a special case of MSSL with a pre-defined precision matrix). Along with a lower error, MSSL captured spatial structure including long range teleconnections between some coastal cities. The ability to harness this spatial structure and task relatedness within a GCM ensembles drives our attention toward MTL in other climate applications.

Consider, in SD, each location in a region as a task with an identical set of possible covariates. These tasks are related through strong unknown spatial dependencies which can be harnessed for SD projections. In the common high dimensional cases of SD, sparse features learned will provide greater significance as presented by [10]. Furthermore, the structure between locations will be learned and may aid projections. MSSL, presented by [16], accounts for sparse feature selection and structure between tasks.

In this study we aim to compare traditional statistical downscaling approaches, BCSD, Multiple Linear Regression, Lasso, and Support Vector Machines, against new approaches in machine learning, Multi-task Sparse Structure Learning and Convolutional Neural Networks (CNNs). During experimentation we apply common training architectures as part of the automated statistical downscaling framework. Results are then analyzed with a variety of metrics including, root mean Square error (RMSE), bias, skill of estimating underlying distributions, correlation, and extreme indices.

2 Statistical Downscaling Methods

2.1 Bias Corrected Spatial Disaggregation

BCSD [46] is widely used in the downscaling community due to its simplicity [1, 5, 45, 30]. Most commonly, GCM data is bias corrected followed by spatial disaggregation on monthly data and then temporally disaggregated to daily projections. Temporal disaggregation is performed by selecting a month at random and adjusting the daily values to reproduce it’s statistical distribution, ignoring daily GCM projections. Thrasher et al. presented a process applying BCSD directly to daily projections [41], removing the step of temporal disaggregation. We the following steps with overlapping a reanalysis dataset and gridded observation data.

1. Bias correction of daily projections using observed precipitation. Observed precipitation is first remapped to match the reanalysis grid. For each day of the year values are pooled, $\pm$ 15 days, from the reanalysis and observed datasets to build a quantile mapping. With the quantile mapping computed, the reanalysis data points are mapped, bias corrected, to the same distribution as the observed data. When applying this method to daily precipitation detrending the data is not necessary because of the lack of trend and is therefore not applied.

2. Spatial disaggregation of the bias-corrected reanalysis data. Coarse resolution reanalysis is then bilinearly interpolated to the same grid as the observation dataset. To preserve spatial details of the fine-grained observations, the average precipitation of each day of the year is computed from the observation and set as scaling factors. These scaling factors are then multiplied to the daily interpolated GCM projections to provide downscaled GCM projections.

2.2 Automated Statistical Downscaling

ASD is a general framework for statistical downscaling incorporating covariate selection and prediction [21]. Downscaling of precipitation using ASD requires two essential steps: 1. Classify rainy/non-rainy days ($\geq$ 1mm), 2. Predict precipitation totals for rainy days. The predicted precipitation can then be written as:

[TABLE]

Formulating $R$ as a binary variable preserves rainy and non-rainy days. We test this framework using five pairs of classification and regression techniques.

2.2.1 Multiple Linear Regression

The simplicity of Multiple Linear Regression (MLR) motivated its use in SD, particularly as part of SDSM [44] and ASD [21]. To provide a baseline relative to the following methods, we apply a variation of MLP using PCA. As discussed previously, PCA is implemented to reduce the dimensionality of a high dimensional feature space by selecting the components that account for a percentage (98% in our implementation) of variance in the data. These principle components, $X$, are used as inputs to classify and predict precipitation totals. We apply a logistic regression model to classify rainy versus non-rainy days. MLP is then applied to rainy days to predict precipitation amounts, $Y$:

[TABLE]

This particular formulation will aid in comparison to  [14] where PCA is coupled with an SVM.

2.2.2 Elastic-Net

Covariate selection can be done in a variety of methods, such as backward stepwise regression and partial correlations. Automatic covariate selection through the use of regularization terms, such as the $L_{1}/L_{2}$ norms in the statistical methods Lasso [42], Ridge [23], and Elastic-Net [51]. Elastic-Net uses a linear combination of $L_{1}/L_{2}$ norms which we will apply in this intercomparison. Given a set of covariates $X$ and observations $Y$, Elastic-net is defined as:

[TABLE]

The $L_{1}$ norm forces uninformative covariate coefficients to zero while the $L_{2}$ norm enforces smoothness while allowing correlated covariates to persist. Cross-validation is applied with a grid-search to find the optimal parameter values for $\lambda_{1}$ and $\lambda_{2}$. High-dimensional Elastic-Net is much less computational than stepwise regression techniques and most often leads to more generalizable models. A similar approach is applied to the classification step by using a logistic regression with an $L_{1}$ normalization term. Previous studies have considered the use of Lasso for SD [19] but to our knowledge, none have considered Elastic-Net.

2.2.3 Support Vector Machine Regression

Ghosh et al. introduced a coupled approach of PCA and Support Vector Machine Regression (SVR) for statistical downscaling [15, 14]. The use of SVR for downscaling aims to capture non-linear effects in the data. As discussed previously, PCA is implemented to reduce the dimensionality of a high dimensional feature space by selecting the components that account for a percentage (98% in our implementation) of variance in the data. Following dimensionality reduction, SVR is used to define the transfer function between the principle components and observed precipitation. Given a set of covariates (the chosen principle components) $X\in\mathbb{R}^{n\times m}$ and $Y\in\mathbb{R}^{n}$ with $d$ covariates and $n$ samples, the support vector regression is defined as [37]:

[TABLE]

where $K(x_{i},x)$ and $w_{i}$ are the kernel functions and their corresponding weights with a bias term $b$. The support vectors are selected during training by optimizing the number of points from the training data to define the relationship between then predictand ($Y$) and predictors ($X$). Parameters $C$ and $\epsilon$ are set during training, which we set to $1.0$ and $0.1$ respectively, corresponding to regularization and loss sensitivity. A linear kernel function is applied to limit overfitting to the training set. Furthermore, support vector classifier was used for classification of rainy versus non-rainy days.

2.2.4 Multi-task Sparse Structure Learning

Recent work in Multi-task Learning aims to exploit structure in the set of predictands while keeping a sparse feature set. Multi-task Sparse Structure Learning (MSSL) in particular learns the structure between predictands while enforcing sparse feature selection ([16]). Goncalves et al. presented MSSL’s exceptional ability to predict temperature through ensembles of GCMs while learning interesting teleconnections between locations ([16]). Moreover, the generalized framework of MSSL allows for implementation of classification and regression models. Applying MSSL to downscaling with least squares regression (logistic regression for classification), we denote $K$ as the number of tasks (observed locations), $n$ as the number of samples, and $d$ as the number of covariates with predictor $X\in\mathbb{R}^{n\times d}$, and predictand $Y\in\mathbb{R}^{n\times K}$. As proposed in [16], optimization over the precision matrix, $\boldsymbol{\Omega}$, is defined as

[TABLE]

where $\boldsymbol{W}\in\mathbb{R}^{d\times K}$ is the weight matrix and $\boldsymbol{\Omega}\in\mathbb{R}^{K\times k}$ is an inverse precision matrix. The $L_{1}$ regularization parameters $\lambda$ and $\gamma$ enforce sparsity over $\boldsymbol{\Omega}$ and $\boldsymbol{W}$. $\boldsymbol{\Omega}$ represents the structure contained between the high resolution observations. Alternating minimization is applied to (5)

1. Initialize $\boldsymbol{\Omega}^{0}=I_{k},\boldsymbol{W}^{0}=\textbf{0}_{d\boldsymbol{X}k}$

2. for t=1,2,3,.. do

[TABLE]

[TABLE]

6 and 7 are independently approximated through Alternating Direction Method of Multipliers (ADMM). Furthermore, by assuming the predictors of each task is identical (as it is for SD), 6 is updated using Distributed-ADMM across the feature space [4].

MSSL enforces similarity between rows of $\boldsymbol{W}$ by learning the structure $\boldsymbol{\Omega}$. For example, two locations which are nearby in space may tend to exhibit similar properties. MSSL will the exploit these properties and impose similarity in their corresponding linear weights. By enforcing similarity in linear weights, we are encouraging smoothness of SD projections between highly correlated locations. $L_{1}$ regularization over $\boldsymbol{W}$ and $\boldsymbol{\Omega}$ jointly encourages sparseness and does not force structure. The parameters encouraging sparseness, $\gamma$ and $\lambda$, are chosen from a validation set using the grid-search technique. These steps are applied for both regression and classification.

2.2.5 Convolutional Neural Networks

Artificial Neural Networks (ANN) have been widely applied to SD with mixed results [40, 36, 5], to name a few. In the past, ANNs had difficulty converging to a local minimum. Recent progress in deep learning has renewed interested in ANNs and are beginning to have impressive results in many applications, including image classification and speech recognition [29, 22, 3]. In particular, Convolutional Neural Networks (CNNs) have greatly impacted computer vision applications by extracting, representing, and condensing spatial properties of the image [29]. SD may benefit from CNN advances by learning spatial representations of GCMs. Though CNNs rely on a high number of samples to reduce overfitting, dropout has been shown to be an effective method of reducing overfitting with limited samples [38]. We note that the number of observations available to daily statistical downscaling may cause overfitting.

CNNs rely on two types of layers, a convolution layer and a pooling layer. In the convolution layer, a patch of size $3\times 3$ is chosen and slid with a stride of 1 around the image. A non-linear transformation is applied to each patch resulting in 8 filters. Patches of size $2\times 2$ are then pooled by selecting the maximum unit with a stride of $2$. A second convolution layer with a $3\times 3$ patch to $2$ filters is followed by a max pooling layer of size $3\times 3$ with stride $3$. The increase of pooling size decreases the dimensionality further. The last pooling layer is then vectorized and densely connected to each high resolution location. This architecture is presented in Figure 1.

Multiple variables and pressure levels from our reanalysis dataset are represented as channels in the CNN input. Our CNN is trained using the traditional back propagation optimization with a decreasing learning rate. During training, dropout with probability 0.5 is applied the densely connected layer. This method aims to exploit the spatial structure contained in the GCM. A sigmoid function is applied to the output layer for classification. To our knowledge, this is the first application of CNNs to statistical downscaling.

2.3 Bias Corrected Spatial Disaggregation with MSSL

To further understand the use of BCSD in Statistical Downscaling, we propose a technique to estimate the errors introduced in BCSD. As presented above, BCSD utilizes a relatively simple quantile mapping approach to statistical downscaling following by interpolation and spatial scaling. Following the BCSD estimates of the observed climate, we compute the presented errors, which may be consistent and have a predictive signal. Modeling such errors using the transfer function approaches above, such as MSSL, may uncover this signal and improve BCSD projections. To apply this technique, the following steps are taken:

Apply BCSD to the coarse scale climate variable and compute the errors. 2. 2.

Excluding a hold out dataset, use MSSL where they predictand is the computed errors and the predictands are from a different set of climate variables, such as Temperature, Wind, Sea Level Pressure, etc. 3. 3.

Subtract the expected errors modeled by step 2 from BCSD projections in step 1.

The transfer function learned in step 2 is then applicable to future observations.

3 Data

The Northeastern United States endures highly variable season and annual weather patterns. Variable climate and weather patterns combined with diverse topology provides difficulty in regional climate projection. Precipitation in particular varies heavily in frequency and intensity seasonally and annually [26]. We choose this region to provide an in-depth comparison of statistical downscaling techniques for daily precipitation and extremes.

3.1 United States Unified Gauge-Based Analysis of Precipitation

High resolution gridded precipitation datasets often provide high uncertainties due to a lack of gauge based observations, poor quality control, and interpolation procedures. Fortunately, precipitation gauge data in the continental United States is dense with high temporal resolution (hourly and daily). The NOAA Climate Prediction Center CPC Unified Gauge-Based Analysis of Precipitation exploits the dense network of rain gauges to provide a quality controlled high resolution (0.25*∘* by 0.25*∘*) gridded daily precipitation dataset from 1948 to the current date. State of the art quality control [7] and interpolation [47] techniques are applied giving us high confidence in the data. We select all locations within the northeastern United States watershed.

3.2 NASA Modern-Era Retrospective Analysis for Research and Applications 2 (MERRA-2)

Reanalysis datasets are often used as proxies to GCMs for statistical downscaling when comparing methods due to their low resolution gridded nature with a range of pressure levels and climate variables. Uncertainties and biases occur in each dataset, but state-of-the-art reanalysis datasets attempt to mitigate these issues. NASA’s MERRA-2 reanalysis dataset [34] was chosen after consideration of NCEP Reanalysis I/II [25] and ERA-Interm [9] datasets. [28] showed the reduced bias of MERRA and ERA-Interm over NCEP Reanalysis II, which is most often used in SD studies. MERRA-2 provides a significant temporal resolution from 1980 to present with relatively high spatial resolution (0.50*∘* by 0.625*∘*). Satellite data provided by NASA’s GEOS-5 project in conjunction with NASA’s data assimilation system when producing MERRA-2 [34].

Only variables available from the CCSM4 GCM model are selected as covariates for our SD models. Temperature, vertical wind, horizontal wind, and specific humidity are chosen from pressure levels 500hpa, 700hpa, and 850hpa. At the surface level, temperature, sea level pressure, and specific humidity are chosen as covariates. To most closely resemble CCSM4, each variable is spatially upscaled to 1.00*∘* to 1.25*∘* at a daily resolution. A large box centralized around the Northeastern Region ranging from 35*∘* to 50*∘* latitude and 270*∘* to 310*∘* longitude is used for each variable. When applying the BCSD model, we use a spatially upscaled Land Precipitation MERRA-2 Reanalysis dataset at a daily temporal resolution. Bilinear interpolation is applied over the coast to allow for quantile mapping of coastal locations as needed.

4 Experiments and Evaluation

In-depth evaluation of downscaling techniques is crucial in testing and understanding their credibility. The implicit assumptions in SD must be clearly understood and tested when applicable. Firstly, SD models assume that the predictors chosen credibly represent the variability in the predictands. This assumption is partially validated through the choice of predictors presented above, which physically represents variability of precipitation. The remainder of the assumption must be tested through experimentation and statistical tests between downscaled projections and observations. The second assumption then requires the statistical attributes of predictands and predictors to be valid outside of the data using for statistical modeling. A hold out set will be used to test the feasibility of this assumption at daily, monthly, and annually temporal resolutions. Third, the climate change signal must be incorporated in the predictands through GCMs. Predictands chosen for this experiment are available through CMIP5 CCSM4 simulations. It is understood that precipitation is not well simulated by GCMs and therefore not used in ASD models [35].

To test these assumptions, we provide in-depth experiments, analysis, and statistical metrics for each method presented above. The years 1980-2004 are used from training and years 2005-2014 are used for testing, taken from the overlapping time period of MERRA-2 and CPC Precipitation. For each method (excluding the special case of BCSD), we chose all covariates from each variable, pressure level, and grid point presented above, totaling 12,781 covariates. Each method applies either dimensionality reduction or regularization techniques to reduce complexity of this high dimensional dataset. Separate models are trained for each season (DJF, MAM, JJA, SON) and used to project the corresponding observations.

Analysis and evaluation of downscaled projections aim to cover three themes:

Ability to capture daily anomalies. 2. 2.

Ability to respond to large scale climate trends on monthly and yearly temporal scales. 3. 3.

Ability to capture extreme precipitation events.

Similar evaluation techniques were applied in recent intercomparison studies of SD [5, 17]. Evaluation of daily anomalies are tested through comparison of bias (Projected - Observed), Root Mean Square Error (RMSE), correlations, and a skill score [33]. The skill score presented by [33] measures how similar two probability density functions are from a range of 0 to 1 where 1 corresponds to identical distributions. Statistics are presented for winter (DJF), summer (JJA), and annually to understand season credibility. Statistics for spring and fall are computed but not presented in order to minimize overlapping climate states and simply results. Each of the measures are computed independently in space then averaged to a single metric. Large scale climate trends are tested by aggregating daily precipitation to monthly and annual temporal scales. The aggregated projections are then compared using Root Mean Square Error (RMSE), correlations, and a skill score as presented in [33]. Due to the limited number of data points in the monthly and yearly projections, we estimate each measure using the entire set of projections and observations.

Climate indices are used for evaluation of SD models’ ability to estimate extreme events. Four metrics from ClimDEX (http://www.clim-dex.org), chosen to encompass a range of extremes, will be utilized for evaluation, as presented by Bürger [5].

CWD - Consecutive wet days $\geq$ 1mm 2. 2.

R20 - Very heavy wet days $\geq$ 20mm 3. 3.

RX5day - Monthly consecutive maximum 5 day precip 4. 4.

SDII - Daily intensity index = Annual total / precip days $\geq$ 1m

Metrics will be computed on observations and downscaled estimates followed by annual (or monthly) comparisons. For example, correlating the maximum number of consecutive wet days per year between observations and downscaled estimates measures each SD models’ ability to capture yearly anomalies. A skill score will also be utilized to understand abilities of reproducing statistical distributions.

5 Results

Results presented below are evaluated using a hold-out set, years 2005-2014. Each model’s ability to capture daily anomalies, long scale climate trends, and extreme events are presented. Our goal is to understand a SD model’s overall ability to provide credible projections rather than one versus one comparisons, therefore statistical significance was not computed when comparing statistics.

5.1 Daily Anomalies

Evaluation of daily anomalies depends on a model’s ability to estimate daily precipitation given the state of the system. This is equivalent to analyzing the error between projections and observations. Four statistical measures are used to evaluate these errors: bias, Pearson Correlation, skill score, and root mean square error (RMSE), as presented in Figure 2, Figure 3, and Table 1). All daily precipitation measures are computed independently in space and averaged to provide a single value. This approach is taken to summarize the measures as simply as possible. Figure 2 shows the spatial representation of annual bias in Table 1.

Overall, methods tend underestimate precipitation annually and seasonally with only PCASVR overestimating. BCSD-MSSL shows the lowest annual and summer bias and second lowest winter bias. BCSD is consistently under projects daily precipitation, but by modeling the possible error with MSSL, bias is reduced. PCAOLS and ELNET are less biased compared to MSSL. CNN has a strong tendency underestimate precipitation. Figure 2 shows consistent negative bias through space for BCSD, ELNET, PCAOLS, MSSL, and CNN while PCASVR shows no discernible pattern.

Correlation measures in Table 1 presents a high linear relationship between projections and observations for the models ELNET (0.64 annually) and MSSL (0.62 annually). We find that BCSD has a lower correlation even in the presence of error correction in BCSD-MSSL. PCASVR provides low correlations, averaging 0.33 annually, but PCAOLS performs substantially better at 0.55.

The skill score is used to measure a model’s ability to reproduce the underlying distribution of observed precipitation where a higher value is better between 0 and 1. BCSD, MSSL, and PCASVR have the largest skill scores, 0.93, 0.92, and 0.91 annually. We find that modeling the errors of BCSD decreases the ability to replicate the underlying distribution. The more basic linear models, PCAOLS and ELNET, present lower skill scores. The much more complex CNN model has difficulty replicating the distribution.

RMSE, presented in Figure 3 and Table 1, measures the overall ability of prediction by squaring the absolute errors. The boxplot in Figure 3, where the box present the quartiles and whiskers the remaining distributions with outliers as points, shows the distribution of RMSE annually over space. The regularized models of ELNET and MSSL have similar error distributions and outperform others. CNN, similar to its under performance in bias, shows a poor ability to minimize error. The estimation of error produced by BCSD-MSSL aids in lowering the RMSE of plain BCSD. PCAOLS reasonably minimizes RMSE while PCASVR severely under-performs compared to all other models. Regression models applied minimize error during optimization while BCSD does not. Seasonally, winter is easier to project with summer being a bit more challenging.

5.2 Large Climate Trends

Analysis of a SD model’s ability to capture large scale climate trends can be done by aggregating daily precipitation to monthly and annual temporal scales. To increase the confidence in our measures, presented in Table II and Figures 4 and 5, we compare all observations and projections in a single computation, rather than separating by location and averaging.

Table 2 and Figure 4 show a wide range of RMSE. A clear difficulty in projecting precipitation in the fall, October in particular, is presented by each time-series in Figure 4. The difference in overall predictability relative to RMSE between the models is evident. BCSD and BCSD-MSSL have significantly lower monthly RMSEs compared to the others. Annually, BCSD-MSSL reduced RMSE by 25% compared to plain BCSD. The linear models, ELNET, MSSL, and PCAOLS, have similar predictability while the non-linear models suffer, CNN being considerably worse.

The skill scores in Table 2 show more difficulty in estimating the annual distribution versus monthly distribution. On a monthly scale BCSD and BCSD-MSSL skill scores outperform all other models but BCSD suffers slightly on an annual basis. However, BCSD-MSSL does not lose any ability to estimate the annual distribution. PCAOLS annual skill score is remarkably higher than the monthly skill score. Furthermore, the three linear models outperform BCSD on an annual basis. PCASVR’s skill score suffers on an annual scale and CNN has no ability to estimate the underlying distribution.

Correlation measures between the models and temporal scales show much of the same. BCSD has the highest correlations in both monthly ( 0.85) and yearly ( 0.64) scales while BCSD-MSSL are slightly lower. CNN correlations fall just behind BCSD and BCSD-MSSL. PCASVR fails with correlation values of 0.22 and 0.18. ELNET has slightly higher correlations in relation to MSSL and PCAOLS.

5.3 Extreme Events

A SD model’s ability to downscale extremes from reanalysis depends on both the response to observed anomalies and ability to reproduce the underlying distribution. Resulting correlation measures present the response to observed anomalies, shown in Figure 6 and Table 3. We find that BCSD has higher correlations for three metrics, namely consecutive wet days, very heavy wet days, and daily intensity index along with a similar results from 5-day maximum precipitation. Furthermore, modeling BCSD’s expected errors with BCSD-MSSL decreases the ability to estimate the chosen extreme indices. Non-linear methods, PCASVR and CNN, suffer greatly in comparison to more basic bias correction and linear approaches. The linear methods, PCAOLS, ELNET, and MSSL, provide similar correlative performance.

A skill score is used to quantify each method’s ability to estimate an indices statistical distribution, presented in Table 3. Contrary to correlative results, PCASVR outperforms the other methods on two metrics, very heavy wet days and daily intensity index, with better than average scores on the other two metrics. BCSD also performs reasonably well in terms of skill scores while BCSD-MSSL suffers from the added complexity. MSSL estimates the number of consecutive wet days well but is less skilled on other metrics. The very complex CNN model has little ability to recover such distributions.

Figure 6 displays a combination of correlative power and magnitude estimate of the daily intensity index. The SDII metric is computed from total annual precipitation and number of wet days. A low SDII metric corresponds to either a relatively large number of estimated wet days or low annual precipitation. We find that the on average methods underestimate this intensity. Based on Figure 5 we see that CNN severely underestimates annual precipitation, causing a low SDII. In contrast, PCASVR overestimates annual precipitation and intensity.

Inconsistent results of PCASVR and CNN indicates that capturing non-linear relationships is outweighed by overfitting. However, BCSD and linear methods are more consistent throughout each metric.

6 Discussion and Conclusion

The ability of statistical downscaling methods to produce credible results is necessary for a multitude of applications. Despite numerous studies experimenting with a wide range of models for statistical downscaling, none have clearly outperformed others. In our study, we experiment with the off-the-shelf applicability of machine learning advances to statistical downscaling in comparison to traditional approaches.

Multi-task Sparse Structure learning, an approach that exploits similarity between tasks, was expected to increase accuracy beyond automated statistical downscaling approaches. We find that MSSL does not provide improvements beyond ELNET, an ASD approach. Furthermore, the parameter set, estimated through cross-validation, attributed no structure aiding prediction.

The recent popularity in deep learning along with it’s ability to capture spatial information, namely Convolutional Neural Networks, motived us to experiment with basic architectures for statistical downscaling. CNNs benefit greatly by implicitly learning abstract non-linear spatial features based on the target variable. This approach proved to poorly estimate downscaled estimates relative to simpler methods. We hypothesize that implicitly learning abstract features rather than preserving the granular feature spaced caused poor performance. More experimentation with CNNs in a different architecture may still provide valuable results.

BCSD, a popular approach to statistical downscaling, outperformed the more complex models in estimating underlying statistical distributions and climate extremes. In many cases, correcting BCSD’s error with MSSL increased daily correlative performance but decreased skill of estimating the distribution. From this result, we can conclude that a signal aiding in prediction was lost during quantile mapping, interpolation, or spatial scaling. Future work may study and improve each step independently to increase overall performance.

Of the seven statistical downscaling approaches studied, the traditional BCSD and ASD methods outperformed non-linear methods, namely Convolutional Neural Network and Support Vector regression, while downscaling daily precipitations. We find that BCSD is skilled at estimating the statistical distribution of daily precipitation, generating better estimates of extreme events. The expectation of CNN and MSSL, two recent machine learning advances which we found most applicable to statistical downscaling, to outperform basic modeled proved false. Improvements and customization of machine learning methods is needed to provide more credible projections.

Acknowledgments

This work was funded by NSF CISE Expeditions in Computing award 1029711, NSF CyberSEES award 1442728, and NSF BIGDATA award 1447587.

MERRA-2 climate reanalysis datasets used were provided by the Global Modeling and Assimilation Office at NASA’s Goddard Space Flight Center. The CPC Unified Gauge-Based Analysis was provided by NOAA Climate Prediction Center.