1 Introduction
Computer simulation models of the physical world, such as numerical weather prediction (NWP) models, are imperfect and can only approximate the complex evolution of physical reality. Some of the errors are due to the uncertainty in the initial and boundary conditions, forcings, and model parameter values. Other errors, called structural model errors, are due to our incomplete knowledge about the true physical processes, and manifest themselves as missing dynamics in the model moosavi2017structual . Examples of structural errors include the misrepresentation of seaice in the spring and fall, errors affecting the stratosphere above polar regions in winter tremolet2007model , as well as errors due to the interactions among (approximatelyrepresented) physical processes.
Data assimilation improves model forecasts by fusing information from both model outputs and observations of the physical world in a coherent statistical estimation framework akella2009different ; le1986variational ; navon1992variational ; tremolet2007model . While traditional data assimilation reduces the uncertainty in the model state and model parameter values, no methodologies to reduce the structural model uncertainty are available to date.
In this study we consider the Weather Research and Forecasting (WRF) model wrf_website , a mesoscale atmospheric modeling system. The WRF model includes multiple physical processes and parametrization schemes, and choosing different model options can lead to significant variability in the model predictions fovell2006impact ; nasrollahi2012assessing .
Among different atmospheric phenomena, the prediction of precipitation is extremely challenging and is obtained by solving the atmospheric dynamic and thermodynamic equations nasrollahi2012assessing . Model forecasts of precipitation are very sensitive to physics options such as the microphysics, cumulus, long wave, and short wave radiation fovell2010influence ; nasrollahi2012assessing ; lowrey2008assessing . Other physics settings that can affect the WRF precipitation predictions include surface physics, planetary boundary layer (PBL), landsurface (LS) parameterizations, and lateral boundary condition. Selecting the right physical process representations and parameterizations is a challenge. In practice the values of physical parameters are empirically determined such as to minimize the difference between the measurements and model predictions wrf_website ; lowrey2008assessing .
Considerable effort has been dedicated to determining the best physical configurations of the weather forecast models such as to improve their predictions of precipitation. No single choice of physical parameters works perfectly for all times, geographical locations, or meteorological conditions gallus1999eta ; wang1997comparison . Lowrey and Yang lowrey2008assessing investigated the errors in precipitation predictions caused by different parameters including microphysics and cumulus physics, the buffer zone, the initialization interval, the domain size and the initial and boundary conditions. Jankov et al. jankov2007impact examined different combinations of cumulus convection schemes, microphysical options, and boundary conditions. They concluded that no configuration was the clear winner at all times, and the variability of precipitation predictions was more sensitive to the choice of the cumulus options rather than microphysical schemes. Another study conducted by Nasrollahi nasrollahi2012assessing showed that the best model ability to predict hurricanes was achieved using a particular cumulus parameterization scheme combined with a particular microphysics scheme. Therefore, the interactions of different physical parameterizations have a considerable impact on model errors, and can be considered as one of the main sources of uncertainty that affect the forecast accuracy.
This paper demonstrates the potential of machine learning techniques to help solve two important problems related to the structural/physical uncertainty in numerical weather prediction models. he first problem addressed herein is the estimation of systematic model errors in output quantities of interest at future times, and the use of this information to improve the model forecasts. The second problem considered is the identification of those specific physical processes that contribute most to the forecast uncertainty in the quantity of interest under specified meteorological conditions.
The application of machine learning techniques to problems in environmental science has grown considerably in recent years. In gilbert2010machine
a kernel based regression method is developed as a forecasting approach with performance close to Ensemble Kalman Filter (EnKF) and less computational resources. Krasnopol et al.
(krasnopol_sky2011development, ) employ an Artificial Neural Network technique for developing an ensemble stochastic convection parameterization for climate models. Attia et al. attia2016clusterHMC develop a new filtering algorithm called Cluster Hybrid Monte Carlo sampling filter (CLHMC) nonGaussian data assimilation which relaxes the Gaussian assumptions by employing a clustering step. Moosavi et al. moosavi2018localization use regression machine learning techniques for adaptive localization in ensemble based data assimilation.This study focuses on the uncertainty in forecasts of cumulative precipitation caused by imperfect representations of physics and their interaction in the WRF model. The total accumulated precipitation includes all phases of convective and nonconvective precipitation. Specifically, we seek to use the discrepancies between WRF forecasts and measured precipitation levels in the past in order to estimate in advance the WRF prediction uncertainty. The modelobservation differences contain valuable information about the error dynamics and the missing physics of the model. We use this information to construct two probabilistic functions. The first one maps the discrepancy data and the physical parameters onto the expected forecast errors. The second maps the forecast error levels onto the set of physical parameters that are consistent with them. Both maps are constructed using supervised machine learning techniques, specifically, using Artificial Neural Networks and Random Forests murphy2012machine . The two probabilistic maps are used to address the problems posed above, namely the estimation of model errors in output quantities of interest at future times, and the identification of physical processes that contribute most to the forecast uncertainty.
The remainder of this study is organized as follows. Section 2 covers the definition of the model errors. Section 3 describes the proposed approach of error modeling using machine learning. Section 4 reports numerical experiments with the WRF model that illustrate the capability of the new approach to answer two important questions regarding model errors. Conclusions are drawn in Section 5.
2 Model errors
Firstprinciples computer models capture our knowledge about the physical laws that govern the evolution of a real physical system. The model evolves an initial state at the initial time to states at future times. All models are imperfect, e.g., atmospheric model uncertainties are associated with subgrid modeling, boundary conditions, and forcings. All these modeling uncertainties are aggregated into a component that is generically called model error Glimm_2004 ; Orrell_2001 ; Palmer_2005 . In the past decade there has been a considerable scientific effort to incorporate model errors and estimate their impact on the best estimate in both variational and statistical approaches akella2009different ; cardinali2014representing ; hansen2002accounting ; rao2015posteriori ; tr2006accounting ; tremolet2007model ; zupanski2006model .
In what follows, we describe our mathematical formulation of the model error associated with NWP models. A similar formulation has been used in moosavi2017structual where the model structural uncertainty is studied based on the information provided by the discrepancy between the model solution and the true state of the physical system, as measured by the available observations.
Consider the following NWP computer model , that describes the timeevolution of the state of the atmosphere:
(1a)  
The state vector contains the dynamic variables of the atmosphere such as temperature, pressure, precipitation, tracer concentrations, at all spatial locations covered by the model, and at . All the physical parameters of the model are lumped into . 
Formally, the true state of the atmosphere can be described by a physical process with internal states , which are unknown. The atmosphere, as an abstract physical process, evolves in time as follows:
(1b) 
The model state seeks to approximates the physical state:
(1c) 
where the operator maps the physical space onto the model space, e.g., by sampling the continuous meteorological fields onto a finite dimensional computational grid moosavi2017structual .
Assume that the model state at has the ideal value obtained from the true state via (1c). The model prediction at will differ from the reality:
(2) 
where the discrepancy between the model prediction and reality is the structural model error. This vector lives in the model space.
Although the global physical state is unknown, we obtain information about it by measuring of a finite number of observables , as follows:
(3) 
Here is the observation operator that maps the true state of atmosphere to the observation space, and the observation error
is assumed to be normally distributed.
In order to relate the model state to observations we also consider the observation operator that maps the model state onto the observation space; the modelpredicted values of the observations (3) are:
(4) 
We note that the measurements and the predictions live in the same space and therefore can be directly compared. The difference between the observations (6b) of the real system and the model predicted values of these observables (4) represent the model error in observation space:
(5) 
For clarity, in what follows we make the following simplifying assumptions moosavi2017structual :

the physical system is finite dimensional ,

the model state lives in the same space as reality, i.e., and is the identity operator in (1c), and
These assumptions imply that the discretization errors are very small, and that the main source of error are the parameterized physical processes represented by and the interaction among these processes. Uncertainties from other sources, such as boundary conditions, are assumed to be negligible.
With these assumptions, the evolution equations for the physical system (1b) and the physical observations equation (3) become, respectively:
(6a)  
(6b) 
The model errors (2) are not fully known at any time , as having the exact errors is akin to having a perfect model. However, the discrepancies between the modeled and measured observable quantities (5) at past times have been computed and are available at the current time .
Our goal is to use the errors in observable quantities at past times, for , in order to estimate the model error at future times . This is achieved by unravelling the hidden information in the past values. Good estimates of the discrepancy , when available, could improve model predictions by applying the correction (6a) to model results:
(7) 
Our proposed error modeling approach constructs inputoutput mappings to estimate given aspects of model errors . The inputs to these mappings are the physical parameters of the model. The outputs to these mappings are different aspects of the error in a quantity of interest, such as the model errors over a specific geographical location, or the error norm of model error integrated over the entire domain.
Specifically, the aspect of interest (quantity of interest) in this study is the error in precipitation levels forecasted by the model. The parameters describe the set of physical processes that are essential to be included in the WRF model in order to produce accurate precipitation forecasts. The WRF model is modular and different combinations of the physical packages can be selected, each corresponding to a different value of .
We use the error mappings learned from past model runs to estimate the model errors . We also consider estimating what combination of physical processes leads to lower model errors, or reversely, what interactions of which physics cause larger errors in the prediction of the quantity of interest.
3 Approximating model errors using machine learning
We propose a multivariate inputoutput learning model to predict the model errors , defined in (2), stemming from the uncertainty in parameters . To this end, we define a probabilistic function that maps every set of input features to output target variables :
(8) 
and approximate the function using machine learning.
Different particular definitions of in (8) will be used to address two different problems related to model errors, as follows:

The first problem is to estimate the systematic model error in certain quantities of interest at future times, and to use this information in order to improve the WRF forecast. To achieve this one quantifies the model error aspects that correspond to running WRF with different physical configurations (different parameters ).

The second problem is to identify the specific physical processes that contribute most to the forecast uncertainty in the quantity of interest under specified meteorological conditions. To achieve this one finds the model configurations (physical parameters ) that lead to forecast errors smaller that a given threshold under specified meteorological conditions.
In what follows we explain in detail the function specification, the input features, and the target variables for each of these problems.
3.1 Problem one: estimating in advance aspects of interest of the model error
Forecasts produced by NWP models are contaminated by model errors. These model errors are highly correlated in time; hence historical information about the model errors can be used as an input to the learning model to gain insight about model errors that affect the forecast. We are interested in the uncertainty caused due to the interaction between the various components in the physics based model; these interactions are lumped into the parameter that is supplied as an input to the learning model. The learning model aims to predict the error of NWP model of next forecast window using the historical values of model error and the physical parameters used in the model. We define the following mapping:
(9) 
We use a machine learning algorithm to approximate the function . The learning model is trained using a dataset that consists of the following inputs:

WRF physical packages that affect the physical quantity of interest (),

historical WRF forecasts ( for ),

historical model discrepancies ( for ),

WRF forecast at the current time (),

the available model discrepancy at the current time () since we have access to the observations from reality at the current time step.
In supervised learning process, the learning model identifies the effect of physical packages, the historical WRF forecast, the historical model discrepancy, and the WRF forecast at the current time on the available model discrepancy at the current time. After the model get trained on the historical data, it yields an approximation to the mapping
. We denote this approximate mapping by .During the test phase the approximate mapping is used to estimate the model discrepancy in advance. We emphasize that the model prediction (WRF forecast) at the time of interest () is available, where as the model discrepancy is an unknown quantity. In fact the run time of WRF is much smaller than the time interval between and , or in other way, the time interval is large enough to run the WRF model and obtain the forecast for next time window, estimate the model errors for next time window and finally improve the model forecast by combining the model forecast and model errors.
At the test time we predict the future model error as follows:
As explained in moosavi2017structual , the predicted error in the observation space can be used to estimate the error in the model space. In order to achieve this one needs to use additional information about the structure of the model and the observation operator. For example, if the error represents the projection of the full model error onto the observation space, we have:
(10a)  
where we use the linearized observation operator at the current time, . A more complex approach is to use a Kalman update formula:  
(10b) 
where is the covariance of observation errors. The Kalman update approach requires estimates of the covariance matrices between model variables; such covariances are already available in an ensemble based data assimilation system. Once we estimate the future model error , we can improve the NWP output using equation (7).
3.2 Problem two: identifying the physical packages that contribute most to the forecast uncertainty
Typical NWP models incorporate an array of different physical packages to represent multiple physical phenomena that act simultaneously. Each physical package contains several alternative configurations (e.g., parameterizations or numerical solvers) that affect the accuracy of the forecasts produced by the NWP model. A particular scheme in a certain physical package best captures the reality under some specific conditions (e.g., time of the year, representation of seaice, etc.). The primary focus of this study is the accuracy of precipitation forecasts, therefore we seek to learn the impacts of all the physical packages that affect precipitation. To this end, we define the following mapping:
(11) 
that estimates the configuration of the physical packages such that the WRF run generates a forecast with an error consistent with the prescribed level (where defined in equation (5) is the forecast error in observation space at time .)
We train the model to learn the effect of the physical schemes on the mismatch between WRF forecasts and reality. The input data required for the training process is obtained by running the model with various physical package configurations , and comparing the model forecast against the observations at all past times to obtain the corresponding errors for and . The output data is the corresponding physical combinations that leads to the input error threshold.
In order to estimate the combinations of physical process configuration that contribute most to the uncertainty in predicting precipitation we take the following approach. The dataset consisting of the observable discrepancies during the current time window is split into a training part and a testing part. In the test phase we use the approximated function to estimate the physical process settings that are consistent with the observable errors . Here we select for each . Note that in this case, since we know what physics has been used for the current results, one can take to be the real parameter values used to generate the test data. However, in general, one selects in an applicationspecific way and the corresponding parameters need to be estimated.
Next, we reduce the desired forecast error level to , and use the approximated function to estimate the physical process setting that corresponds to this more accurate forecast. To identify the package setting that has the largest impact on the observable error we monitor the variability in the predicted parameters . Specifically, the number of times the setting of a physical process in is different from its setting in is an indicator of the variability in model prediction when that package is changed. A higher variability in predicted physical packages implies a larger contribution towards the model errors  as estimated by the ML model.
3.3 Machine learning algorithms
In order to approximate the functions (9) and (11) discussed earlier we use regression machine learning methods. Choosing a right learning algorithm to use is challenging as it largely depends on the problem and the data available attia2016clusterHMC ; asgari2017Latent ; moosavi2017multivariate ; moosavi2018localization . Here, we use Random Forests (RF) and Artificial Neural Networks (ANN) as our learning algorithms murphy2012machine . Both RF and ANN algorithms tan handle nonlinearity in regression and classification. Given that the physical phenomena governing precipitation are highly nonlinear, and and atmospheric dynamics is chaotic, we believe that RF and ANN approaches are well suited to capture the associated features. We briefly review these techniques next.
3.3.1 Random forests
A random forest breiman2001random
is an ensemble based method that constructs multiple decision trees. The principle idea behind ensemble methods is that a group of weak learners can come together to form a strong learner
breiman1996bagging ; breiman2001random . The decision tree is built topdown from observations of target variables. The observation dataset is partitioned, smaller subsets are represented in branches, and decisions about the target variables are represented in the leaves.There are many specific decisiontree algorithms available, including ID3 (Iterative Dichotomiser 3) quinlan1986induction , C4.5 (successor of ID3) quinlan2014c4 , CART (Classification And Regression Tree), CHAID (CHisquared Automatic Interaction Detector), and conditional inference trees strobl2008conditional . If the dataset has multiple attributes, one can decide which attribute to place at the root or at different levels of the tree by considering different criteria such as information gain or the gini index ceriani2012origins .
Trees can be nonrobust, with small changes in the tree leading to large changes in regression results. Moreover, trees tend to overfit the data segal2004machine
. The random forest algorithm uses the bagging technique for building an ensemble of decision trees which are accurate and powerful at handling large, high dimensional datasets. Moreover, the bagging technique greatly reduces the variance
dietterich2000ensemble . For each tree in the forest, a bootstrap sample breiman1996bagging ; dietterich2000ensemble is selected from the dataset and instead of examining all possible featuresplits, some subset of the features is selected liaw2002classification . The node then splits on the best feature in the subset. By using a random sample of features the correlation between trees in the ensemble decreases, and the learning for each tree is much faster by restricting the features considered for each node.3.3.2 Artificial neural networks
ANN is a computational model inspired by human brain’s biological structure. ANN consist of neurons and connections between the neurons (weights) which are organized in layers. At least three layers of neurons (an input layer, a hidden layer, and an output layer) are required for construction of a neural network, where the input layer distributes the input signals to the first hidden layer. The feedforward operation in a network passes information to neurons in a subsequent hidden layer. The neurons combine this information, and the output of each layer is obtained by passing the combined information through a differentiable transfer function that can be logsigmoid, hyperbolic tangent sigmoid, or linear transfer function.
In supervised learning the network is provided with samples from which it discovers the relations of inputs and outputs. The learning problem consists of finding the optimal parameters of network such that the error between the desired output and the output signal of the network is minimized. The network first is initialized with randomly chosen weights and then the error is backpropagated through the network using a gradient descent method. The gradient of the error function is computed and used to modify weights and biases such that the error between the desired output and the output signal of the network is minimized funahashi1989approximate ; rumelhart1985learning . This process is repeated iteratively until the network output is close to the desired output haykin2009neural .
4 Numerical experiments
We apply the proposed learning models to the Weather Research and Forecasting model wrf_website in order to:

predict the bias in precipitation forecast caused by structural model errors,

predict the statistics associated with the precipitation errors, and

identify the specific physics packages that contribute most to precipitation forecast errors for given meteorological conditions.
4.1 The WRF model
In this study we use the nonhydrostatic WRF model version 3.3. The simulation domain, shown in Fig. 1, covers the continental United States and has dimensions of horizontal grid points in the westeast and southnorth directions respectively, with a horizontal grid spacing of wang2014downscaling . The grid has 60 vertical levels to cover the troposphere and lower part of the stratosphere between the surface to approximately . In all simulations, the 6hourly analysis from the National Centers for Environmental Prediction (NCEP) are used as the initial and boundary conditions of the model NCEP_data . The stage IV estimates are available at an hourly temporal resolution over continental United States. For experimental purposes, we use the stage IV NCEP analysis as a proxy for the true state of the atmosphere. The simulation window begins at UTC (Universal Time Coordinated) on May 1st 2017, and the simulation time is a six hour window time the same day. The “true” states of the atmosphere are available using the NCEP analysis data hourly. All the numerical experiments use the NCEP analysis data to run WRF model on May 1st 2017.
The model configuration parameters represent various combinations of microphysics schemes, cumulus parameterizations, short wave, and long wave radiation schemes. Specifically, each process is represented by the schema values of each physical parameter it uses, as detailed in WRF model physics options and references wrf_physics . The microphysics option provides atmospheric heat and moisture tendencies in atmosphere which also accounts for the vertical flux of precipitation and the sedimentation process. The cumulus parameterization is used to vertically redistribute heat and moisture independent of latent heating due to precipitation. The long wave radiation considers clearsky and cloud upward and downward radiation fluxes and the short wave radiation considers clearsky and cloudy solar fluxes.
A total number of 252 combinations of the four physical modules are used in the simulations. The microphysics schemes include: Kessler kessler1995continuity , Lin lin1983bulk , WSM3 Hong hong2004revised , WSM5 Hong hong2004revised , Eta (Ferrier), WSM6 hong2006wrf , Goddard tao1989ice , Thompson thompson2008explicit , Morrison morrison2009impact . The cumulus physics schemes applied are: KainFritsch kain2004kain , BettsMillerJanjic janjic1994step , Grell Freitasgrell2013scale . The long wave radiation physics include: RRTM mlawer1997radiative , Cam conley2012description . Short wave radiation physics include: Dudhia dudhia1989numerical , Goddard chou1999solar , Cam conley2012description .
For each of the 252 different physics combinations, the effect of each physics combination on precipitation is investigated. The NCEP analysis grid points are , while the WRF computational model have
grid points. For obtaining the discrepancy between the WRF forecast and NCEP analysis we linearly interpolate the analysis to transfer the physical variables onto the model grid. Figure
1(a) and 1(b) shows the NCEP analysis at and on 5/1/2017 which are used as initial condition and “true” (verification) state, respectively. The WRF forecast corresponding to the physics microphysics: Kessler, cuphysics: KainFritsch, ralwphysics: Cam , raswphysics: Dudhia is illustrated in Figure 1(c). Figure 2 shows contours of discrepancies at discussed in equation (5) for two different physical combinations, which illustrates the effect that changing the physical schemes has on the forecast.4.2 Experiments for problem one: predicting pointwise precipitation forecast errors over a small geographic region
We demonstrate our learning algorithms to forecast precipitation in the state of Virginia on May 1st 2017 at . Our goal is to use the learning algorithms to correct the bias created due to model errors and hence improve the forecast for precipitation. As described in section 3.1, we learn the function of equation (9) using the training data from the previous forecast window ( to ):
We use two learning algorithms to approximate the function , namely, the RF and ANN using Scikitlearn, machine learning library in Python scikitlearn
. The RF with ten trees and CART learning tree algorithm is used. The ANN with six hidden layers and hyperbolic tangent sigmoid activation function in each layer and linear activation function at last layer is employed. The number of layers and number of neurons in each layer are tuned empirically. For training purposes, we use the NCEP analysis of the May 1st 2017 at
as initial conditions for the WRF model. The forecast window is 6 hours and the WRF model forecast final simulation time is . The input features are:
The physics combinations ().

The hourly WRF forecasts projected onto observation space , . The WRF state () includes all model variables such as temperature, pressure, precipitation, etc. The observation operator extracts the precipitation portion of the WRF state vector, . Accordingly, is the discrepancy between WRF precipitation forecast and the observed precipitation .

The observed discrepancies at past times (, ).
The output variable is the discrepancy between the NCEP analysis and the WRF forecast at , i.e., the observable discrepancies for the current forecast window (). In fact, for each of the different physical configurations, the WRF model forecast as well as the difference between the WRF forecast and the analysis are provided as inputoutput combinations for learning the function . The number of grid points over the state of Virginia is . Therefore for each physical combination we have grid points, and the total number of samples in the training data set is with features.
Both ANN and RF are trained with the above inputoutput combinations described above and during the training phase, the learning model learns the effect of interaction between different physical configurations on the WRF forecast and model error and obtains the approximation to the function which we denote by . The goal is to have more accurate forecast in the future time windows. We don’t have the analysis data of future time windows but we can run WRF for future time windows and also predict the future model error using the approximated function . Once we obtain the predicted model error we can use that information in order to raise the accuracy of WRF forecast. In the testing phase we use the function to predict the future forecast error given the combination of physical parameters as well as the WRF forecast at time as input features.
To quantify the accuracy of the predicted error we calculate the Root Mean Squared Error (RMSE) between the true and predicted discrepancies at :
(12) 
where is the number of grid points over Virginia,. is the predicted discrepancy in the grid point, and is the actual discrepancy in the grid point. The actual discrepancy is obtained as the difference between the NCEP analysis and WRF forecast at time . This error metric is computed for each of the different configurations of the physics. The minimum, maximum and average RMSE over the runs is reported in Table 1.
ANN  

RF 
The predicted discrepancy in the observation space can be used to approximate the discrepancy in the model space using equation (10). Here all the grid points are observed and therefore the error in the model space equal to the error in the observation space. Next, the estimate forecast error can be used to correct the forecast bias caused by model errors using (7), and hence to improve the forecast at : . Figure 3(a) shows the WRF forecast for for the state of Virginia using the following physics packages (the physics options are given in parentheses):

Microphysics (Kessler),

Cumulusphysics (Kain),

Shortwave radiation physics (Dudhia),

Longwave radiation physics (Janjic).
Figure 3(b) shows the NCEP analysis at time , which is our proxy for the true state of the atmosphere. The discrepancy between the NCEP analysis and the raw WRF forecast is shown in the Figure 4(a). Using the model error prediction we can improve the WRF result by adding the predicted bias to the WRF forecast. The discrepancy between the corrected WRF forecast and the NCEP analysis is shown in the Figure 4(b). The results show a considerable reduction of model errors as compared to the uncorrected forecast of Figure 4(a). Table 2 shows the minimum and average of original model error vs the improved model errors.
Original forecast  

Improved forecast 
4.3 Experiments for problem one: predicting the norm of precipitation forecast error over the entire domain
We now seek to estimate the twonorm of precipitation model error over the entire continental U.S., which gives a global metric for the accuracy of the WRF forecast, and helps provide insight about the physics configurations that result in more accurate forecasts. To this end the following mapping is constructed:
To build the training dataset, we run WRF with each of the different physical configurations. The forecast window is 6 hours and the WRF model forecast final simulation time is at . The hourly WRF forecast and discrepancy between the analysis and WRF forecast is used as training features.
The input features are:

different physics schemes (),

the norms of the WRF model predictions at previous time windows, as well as at the current time (, ), and

the norms of past observed discrepancies (, ).
The output variable is the norm of the discrepancy between WRF precipitation prediction and the NCEP precipitation analysis for the current time window ().
We use two different learning algorithms, namely, RF with ten trees in the forest and ANN with four hidden layers, the hyperbolic tangent sigmoid activation function in each layer and linear activation function at last layer. The number of layers and neurons at each layer is tuned empirically. The total number of samples in the training set is with of features. During the training phase the model learns the effect of interaction of different physical configurations on model error and obtains the approximated function .
In the test phase we feed the approximated function the model information from to the endpoint of the next forecast window to predict the norm of the model error .
Validation of the learned error mapping
Table 3 shows the RMSE between the actual and predicted norms of discrepancies for ANN and RF. The RMSE is taken over the 252 runs with different physics combinations. Both learning models are doing well, with the ANN giving slightly better results than the RF.
ANN  

RF 
Analysis of the best combination of physical packages
Based on our prediction of the norm of model error, the best physics combination that leads to lowest norm of precipitation error over the entire continental U.S. for the given meteorological conditions is:

the BMJ cumulus parameterization, combined with

the WSM5 microphysics,

Cam long wave, and

Dudhia short wave radiation physics.
According to the true model errors, the best physics combination leading to the lowest norm of model error is achieved using the BMJ cumulus parameterization, combined with the WSM5 microphysics, Cam long wave, and Cam short wave radiation physics.
4.4 Experiments for problem two: identify the physical processes that contribute most to the forecast uncertainty
The interaction of different physical processes greatly affects precipitation forecast, and we are interested in identifying the major sources of model errors in WRF. To this end we construct the physics mapping (11) using the norm and the statistical characteristics of the modeldata discrepancy (over the entire U.S.) as input features:
Statistical characteristics include the mean, minimum, maximum, and variance of the filed across all grid points over the continental U.S. Note that this is slightly different than (11) where the inputs are the raw values of these discrepancies for each grid point. The output variable is the combination of physical processes that leads to model errors consistent with the input pattern and .
To build the dataset, the WRF model is simulated for each of the different physical configurations, and the mismatches between the WRF forecasts and the NCEP analysis at the end of the current forecast window are obtained. Similar to the previous experiment, the initial conditions used in the WRF model is the NCEP analysis for the May 1st 2017 at . The forecast window is 6 hours and the WRF model forecast is obtained for time . The discrepancy between the NCEP analysis at and WRF forecast at forms the observable discrepancy for the current forecast window . For each of the 252 different physical configurations, this process is repeated and statistical characteristics of the WRF forecast model error , and the norm of model error are used as feature values of the function .
Validation of the learned physics mapping
From all the collected data points, (202 samples) are used for training the learning model, and the remaining (50 samples) are used for testing purposes.
The RF has default ten trees in the forest and ANN has four hidden layers and hyperbolic tangent sigmoid activation function in each layer with linear activation function at last layer. The number of layers and neurons at each layer is tuned empirically. The learning model uses the training dataset to learn the approximate mapping . This function is applied to the each of the test samples to obtain the predicted physical combinations . In order to evaluate these predictions, we run the WRF model again with the physical setting and obtain the new forecast , and the corresponding observable discrepancy . The RMSE between the norm of actual observable discrepancies and the norm of predicted discrepancies are shown in Table 4. The small values of the difference demonstrates the performance of the learning algorithm.
ANN  

RF 
Analysis of variability in physical settings
We repeat the test phase for each of the test samples with the scaled values of observable discrepancies as inputs, and obtain the predicted physical combinations . Large variability in the predicted physical settings indicate that the respective physical packages variability have a strong influence on the WRF forecast error. We count the number of times the predicted physics is different from when the input data spans the entire test data set.
The results shown in Figure 5 indicate that microphysics and cumulus physics are not too sensitive to the change of input data, whereas shortwave and longwave radiation physics are quite sensitive to changes in the input data. Therefore our learning model indicates that having an accurate shortwave and longwave radiation physics package will aid in greatly reducing the uncertainty in precipitation forecasts due to missing/incorrect physics.
5 Conclusions
This study proposes a novel use of machine learning techniques to understand, predict, and reduce the uncertainty in the WRF model precipitation forecasts due to the interaction of several physical processes included in the model.
We construct probabilistic approaches to learn the relationships between the configuration of the physical processes used in the simulation and the observed model forecast errors. These relationships are then used to solve two important problems related to model errors, as follows: estimating the systematic model error in a quantity of interest at future times, and identifying the physical processes that contribute most to the forecast uncertainty in a given quantity of interest under specified conditions.
Numerical experiments are carried out with the WRF model using the NCEP analyses as a proxy for the real state of the atmosphere. Ensembles of model runs with different parameter configurations are used to generate the training data. Random forests and Artificial neural network models are used to learn the relationships between physical processes and forecast errors. The experiments validate the new approach, and illustrates how it is able to estimate model errors, indicate best model configurations, and pinpoint to those physical packages that influence most the WRF prediction accuracy.
While the numerical experiments are done with WRF, and are focused on forecasting precipitation, the methodology developed herein is general and can be applied to the study of errors in other models, for other quantities of interest, and for learning additional relationships between model physics and model errors.
Acknowledgments
This work was supported in part by the projects AFOSR DDDAS 15RT1037 and AFOSR Computational Mathematics FA95501710205 and by the Computational Science Laboratory at Virginia Tech. The authors would like to thank Dr. Răzvan Ştefănescu for his valuable assistance and suggestions regarding WRF runs and the NCEP dataset.
References
References
 [1] Santha Akella and Ionel M Navon. Different approaches to model error formulation in 4DVar: a study with highresolution advection schemes. Tellus A, 61(1):112–128, 2009.
 [2] Elham Asgari and Kaveh Bastani. The utility of Hierarchical Dirichlet Processfor relationship detection of latent constructs. In Academy of Management Proceedings, 2017.
 [3] Ahmed Attia, Azam Moosavi, and Adrian Sandu. Cluster sampling filters for nonGaussian data assimilation. arXiv preprint arXiv:1607.03592, 2016.
 [4] Leo Breiman. Bagging predictors. Machine learning, 24(2):123–140, 1996.
 [5] Leo Breiman. Random forests. Machine learning, 45(1):5–32, 2001.
 [6] Carla Cardinali, Nedjeljka Žagar, Gabor Radnoti, and Roberto Buizza. Representing model error in ensemble data assimilation. Nonlinear Processes in Geophysics, 21(5):971–985, 2014.
 [7] Lidia Ceriani and Paolo Verme. The origins of the Gini index: extracts from variabilità e mutabilità (1912) by Corrado Gini. The Journal of Economic Inequality, 10(3):421–443, 2012.
 [8] MingDah Chou and Max J Suarez. A solar radiation parameterization (cliradsw) for atmospheric studies. NASA Tech. Memo, 10460:48, 1999.
 [9] Andrew J Conley, Rolando Garcia, Doug Kinnison, JeanFrancois Lamarque, Dan Marsh, Mike Mills, Anne K Smith, Simone Tilmes, Francis Vitt, Hugh Morrison, et al. Description of the NCAR community atmosphere model (CAM 5.0). NCAR technical note, 2012.

[10]
Thomas G Dietterich et al.
Ensemble methods in machine learning.
Multiple classifier systems
, 1857:1–15, 2000.  [11] Jimy Dudhia. Numerical study of convection observed during the winter monsoon experiment using a mesoscale twodimensional model. Journal of the Atmospheric Sciences, 46(20):3077–3107, 1989.
 [12] Robert G Fovell. Impact of microphysics on hurricane track and intensity forecasts. In Preprints, 7th WRF Users’ Workshop, NCAR, 2006.
 [13] Robert G Fovell. Influence of cloudradiative feedback on tropical cyclone motion. In 29th Conference on Hurricanes and Tropical Meteorology, 2010.
 [14] KenIchi Funahashi. On the approximate realization of continuous mappings by neural networks. Neural networks, 2(3):183–192, 1989.
 [15] WA Gallus Jr. Eta simulations of three extreme rainfall events: Impact of resolution and choice of convective scheme. Wea. Forecasting, 14:405–426, 1999.
 [16] Robin C Gilbert, Michael B Richman, Theodore B Trafalis, and Lance M Leslie. Machine learning methods for data assimilation. Computational Intelligence in Architecturing Complex Engineering Systems, pages 105–112, 2010.
 [17] J. Glimm, S. Hou, Y.H. Lee, D.H. Sharp, and K. Ye. Sources of uncertainty and error in the simulation of flow in porous media. Computational & Applied Mathematics, 23:109–120, 2004.
 [18] Georg A Grell and Saulo R Freitas. A scale and aerosol aware stochastic convective parameterization for weather and air quality modeling. Atmospheric Chemistry & Physics Discussions, 13(9), 2013.
 [19] James A Hansen. Accounting for model error in ensemblebased state estimation and forecasting. Monthly Weather Review, 130(10):2373–2391, 2002.
 [20] S.S. Haykin. Neural Networks and Learning Machines. Number v. 10 in Neural networks and learning machines. Prentice Hall, 2009.
 [21] SongYou Hong, Jimy Dudhia, and ShuHua Chen. A revised approach to ice microphysical processes for the bulk parameterization of clouds and precipitation. Monthly Weather Review, 132(1):103–120, 2004.

[22]
SongYou Hong and JeongOck Jade Lim.
The WRF singlemoment 6class microphysics scheme (wsm6).
J. Korean Meteor. Soc, 42(2):129–151, 2006.  [23] Zaviša I Janjić. The stepmountain eta coordinate model: Further developments of the convection, viscous sublayer, and turbulence closure schemes. Monthly Weather Review, 122(5):927–945, 1994.
 [24] Isidora Jankov, Paul J Schultz, Christopher J Anderson, and Steven E Koch. The impact of different physical parameterizations and their interactions on cold season QPF in the American River basin. Journal of Hydrometeorology, 8(5):1141–1151, 2007.
 [25] John S Kain. The Kain–Fritsch convective parameterization: an update. Journal of Applied Meteorology, 43(1):170–181, 2004.
 [26] E Kessler. On the continuity and distribution of water substance in atmospheric circulations. Atmospheric research, 38(14):109–145, 1995.
 [27] VM Krasnopol sky, Michael FoxRabinovitz, Alexei Belochitski, Philip J Rasch, Peter Blossey, and Yefim Kogan. Development of neural network convection parameterizations for climate and NWP models using Cloud Resolving Model simulations. US Department of Commerce, National Oceanic and Atmospheric Administration, National Weather Service, National Centers for Environmental Prediction, 2011.
 [28] FrançoisXavier Le Dimet and Olivier Talagrand. Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects. Tellus A: Dynamic Meteorology and Oceanography, 38(2):97–110, 1986.
 [29] Andy Liaw, Matthew Wiener, et al. Classification and regression by random forest. R news, 2(3):18–22, 2002.
 [30] YuhLang Lin, Richard D Farley, and Harold D Orville. Bulk parameterization of the snow field in a cloud model. Journal of Climate and Applied Meteorology, 22(6):1065–1092, 1983.
 [31] Marla R Knebl Lowrey and ZongLiang Yang. Assessing the capability of a regionalscale weather model to simulate extreme precipitation patterns and flooding in central Texas. Weather and Forecasting, 23(6):1102–1126, 2008.
 [32] Eli J Mlawer, Steven J Taubman, Patrick D Brown, Michael J Iacono, and Shepard A Clough. Radiative transfer for inhomogeneous atmospheres: RRTM, a validated correlatedk model for the longwave. Journal of Geophysical Research: Atmospheres, 102(D14):16663–16682, 1997.
 [33] Azam Moosavi, Ahmed Attia, and Adrian Sandu. A machine learning approach to adaptive covariance localization. arXiv preprint arXiv:1801.00548, 2018.
 [34] Azam Moosavi and Adrian Sandu. A statespace approach to analyze structural uncertainty in physical models. Metrologia, 2017.
 [35] Azam Moosavi, Razvan Stefanescu, and Adrian Sandu. Multivariate predictions of local reducedordermodel errors and dimensions. arXiv preprint arXiv:1701.03720, 2017.
 [36] Hugh Morrison, Gregory Thompson, and V Tatarskii. Impact of cloud microphysics on the development of trailing stratiform precipitation in a simulated squall line: Comparison of oneand twomoment schemes. Monthly Weather Review, 137(3):991–1007, 2009.
 [37] Kevin P Murphy. Machine learning: A probabilistic perspective. MIT press, 2012.
 [38] Nasrin Nasrollahi, Amir AghaKouchak, Jialun Li, Xiaogang Gao, Kuolin Hsu, and Soroosh Sorooshian. Assessing the impacts of different WRF precipitation physics in hurricane simulations. Weather and Forecasting, 27(4):1003–1016, 2012.
 [39] I Michael Navon, Xiaolei Zou, J Derber, and J Sela. Variational data assimilation with an adiabatic version of the NMC spectral model. Monthly weather review, 120(7):1433–1446, 1992.
 [40] National Oceanic and Atmospheric Administration (NOAA). https://www.ncdc.noaa.gov/dataaccess/modeldata/modeldatasets/globalforcastsystemgfs.
 [41] D. Orrell, L. Smith, J. Barkmeijer, and T.N. Palmer. Model error in weather forecasting. Nonlinear Processes in Geophysics, 8:357–371, 2001.
 [42] T.N. Palmer, G.J. Shutts, R. Hagedorn, F.J. DoblasReyes, T. Jung, and M. Leutbecher. Representing model uncertainty in weather and climate prediction. Annu. Rev. Earth Planet. Sci, 33:163–93, 2005.
 [43] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay. Scikitlearn: Machine learning in Python. Journal of Machine Learning Research, 12:2825–2830, 2011.
 [44] J. Ross Quinlan. Induction of decision trees. Machine learning, 1(1):81–106, 1986.
 [45] J Ross Quinlan. C4. 5: programs for machine learning. Elsevier, 2014.
 [46] Vishwas Rao and Adrian Sandu. A posteriori error estimates for the solution of variational inverse problems. SIAM/ASA Journal on Uncertainty Quantification, 3(1):737–761, 2015.
 [47] David E Rumelhart, Geoffrey E Hinton, and Ronald J Williams. Learning internal representations by error propagation. Technical report, DTIC Document, 1985.
 [48] Mark R Segal. Machine learning benchmarks and random forest regression. Center for Bioinformatics & Molecular Biostatistics, 2004.
 [49] Carolin Strobl, AnneLaure Boulesteix, Thomas Kneib, Thomas Augustin, and Achim Zeileis. Conditional variable importance for random forests. BMC bioinformatics, 9(1):307, 2008.
 [50] WeiKuo Tao, Joanne Simpson, and Michael McCumber. An icewater saturation adjustment. Monthly Weather Review, 117(1):231–235, 1989.
 [51] Gregory Thompson, Paul R Field, Roy M Rasmussen, and William D Hall. Explicit forecasts of winter precipitation using an improved bulk microphysics scheme. part ii: Implementation of a new snow parameterization. Monthly Weather Review, 136(12):5095–5115, 2008.
 [52] Yannick Tr’emolet. Accounting for an imperfect model in 4DVar. Quarterly Journal of the Royal Meteorological Society, 132(621):2483–2504, 2006.
 [53] Yannick Trémolet. Modelerror estimation in 4DVar. Quarterly Journal of the Royal Meteorological Society, 133(626):1267–1280, 2007.
 [54] Jiali Wang and Veerabhadra R Kotamarthi. Downscaling with a nested regional climate model in nearsurface fields over the contiguous united states. Journal of Geophysical Research: Atmospheres, 119(14):8778–8797, 2014.
 [55] Wei Wang and Nelson L Seaman. A comparison study of convective parameterization schemes in a mesoscale model. Monthly Weather Review, 125(2):252–278, 1997.
 [56] Weather Research Forecast Model. https://www.mmm.ucar.edu/weatherresearchandforecastingmodel.
 [57] WRF Model Physics Options and References. http://www2.mmm.ucar.edu/wrf/users/phys_references.html.
 [58] Dusanka Zupanski and Milija Zupanski. Model error estimation employing an ensemble data assimilation approach. Monthly Weather Review, 134(5):1337–1354, 2006.
Comments
There are no comments yet.