DEnKF – Variational Hybrid Snow Cover Fraction Data Assimilation for Improving Snow Simulations with the Common Land Model

This study assesses the analysis performance of a hybrid DEnKF-variational data assimilation (DA) method (DEnVar) for assimilating the MODIS snow cover fraction (SCF) into the Common Land Model (CoLM). Coupling a deterministic ensemble Kalman filter (DEnKF) with a one-dimensional variational DA method (1DVar), DEnVar without observation perturbations is a two-step DA method. That is, the analysis ensemble mean and analysis error covariance of DEnKF are introduced into the 1DVar hybrid cost function, and the analysis mean of DEnKF is replaced by the 1DVar analysis. The analysis performance of DEnVar was experimentally compared with DEnKF, 1DVar, and EnVar (hybrid ensemble-variational DA) at five sites in the Altay region of China from November 2008 to March 2009. From our results, it is shown that the four DA experiments can improve snow simulations at most sites when the available MODIS SCF is assimilated. The DEnVar experiment using the hybrid error covariance shows the best analysis performance among the four DA experiments at most sites. Furthermore, sensitivity tests show that DEnVar is slightly sensitive to the weighting coefficient, which controls the respective weights of ensembleand (National Meteorological Center) NMC-based error covariances, but is highly sensitive to the observation error. DEnVar obtains better analysis performance when using the ensemble-based analysis error covariance rather than the hybrid error covariance coupling ensemble-based analysis and static NMC-based OPEN ACCESS Remote Sens. 2014, 6 10613 error covariances. The inaccurate distribution of observation error may invalidate the DEnVar method.


Introduction
Seasonal snow cover plays a crucial role in the Earth's hydrological processes, land surface energy balance, and climate [1][2][3] due to its high albedo and low thermal conductivity [4,5].The snow cover fraction (SCF) is a key parameter in the studies of both hydrology and climatology [6].Normally, SCF is derived from snowpack state such as snow depth (SD) or snow water equivalent (SWE).SD and SWE are important internal variables of hydrological and climatic models.SWE is the product of SD and snow density.They have significant impacts on climatic and hydrological simulations and snowmelt runoff predictions [7].Substantial snow accumulation frequently causes severe snow disasters, such as frostbite, death for both animals and people, and traffic disruption [8].Massive and rapid snowmelt may cause flooding.Therefore, the accurate estimates of SD and SWE are important for reduction, mitigation, and prevention of snow-related disasters [9].
In recent decades, a hybrid DA method of utilizing advanced features of both ensemble-based and variational DA methods has been proposed (e.g., [28][29][30][31][32]) for numerical weather prediction (e.g., [33][34][35][36]).For the hybrid method, the ensemble covariance generated from the ensemble-based DA system is easily incorporated into the hybrid cost function of variational DA.The hybrid DA system can conform well to the existing variational framework and be more flexible than the conventional ensemble-based methods for exploring sampling errors and model errors [36].Recent studies have also demonstrated the potential of the hybrid method as compared with the stand-alone EnKF and variational methods (e.g., [32,37,38]).
Given that these hybrid DA methods have been successfully used for improving numerical weather prediction, we hypothesize that a hybrid DA method is certainly beneficial to snow simulations.The stochastic EnKFs refer to perturbing observations for exploring observation uncertainties [39], which may introduce sampling errors and lead to non-optimal DA, especially for small ensembles [40].In the SCF assimilation, perturbing the full SCF observations may decrease the SCF observations, and therefore systematically and substantially underestimate the SWE [15].A feasible alternative to ensemble square root filter (EnSRF) [41] that has been developed is a relatively robust deterministic ensemble Kalman filter (DEnKF) without observation perturbations.As compared with the EnSRF [42], the DEnKF is easier to implement.The DEnKF method has been used to assimilate MODIS SCF into the Common Land Model (CoLM) for SD simulation improvements.However, the assimilated results have shown that the assimilation of MODIS SCF only using the DEnKF method (a sample ensemble size of 30) cannot effectively improve SD simulations [20].To avoid underestimation of analysis ensemble mean and background error covariance of DEnKF with smaller ensembles, a hybrid method (DEnVar) of DEnKF-variational DA without observation perturbations that couples DEnKF and 1DVar methods is proposed to assimilate MODIS SCF observations into the CoLM for improving SD and SWE simulations.As a comparison with the other three DA methods-DEnKF, 1DVar, and EnVar (hybrid ensemble-variational DA)-the DEnVar method was used at five sites in the Altay region of China from November 2008 to March 2009.The EnVar method is similar to the method of Zhang et al. [36].We performed sensitivity tests on the hybrid error covariance weight and the observation errors.
The CoLM and related data assimilation methods are described in Section 2. The experimental design and data sets are introduced in section 3. Section 4 presents a comparison of DEnVar and other DA methods and sensitivity tests on DEnVar.Concluding remarks are given in Section 5.

Common Land Model
A column (1-D) model, the Common Land Model (CoLM), is the latest land surface model for climate studies [43].The CoLM combines the advanced features of three outstanding and well-documented land surface models, including the Land Surface Model (LSM) [44], the Biosphere-Atmosphere Transfer Scheme (BATS) [45], and the 1994 version of the Chinese Academy of Science Institute of Atmospheric Physics LSM (IAP94) [46].In the CoLM, SD and SWE are forecast variables.As a diagnostic variable, snow cover calculated by total SD directly affects energy fluxes and grid-averaged albedo in the CoLM.The CoLM has the five-layer snow scheme which accounts for layer-based liquid water retention, thawing-freezing, snow-melt, and heat energy.Snow compaction accounts for the destructive, pressure, and melt metamorphisms.At each time step, snow layers can be combined or divided in terms of the snow accumulation or melting conditions, like layer thicknesses.More details about snow processes in the CoLM can be found in the works of Dai et al. [47] and Dai et al. [43].

Variational Method
In traditional variational assimilation, the analysis is performed by minimizing a cost function J [28]: (1) where x = [snowdp, swe] stands for the analysis state vector (e.g.snow depth and snow water equivalent in the CoLM), and x b for the background state vector.B stands for the static background error covariance matrix, y o for the MODIS SCF observation vector, R for the error variance of MODIS SCF observations, and H for the nonlinear observation operator mapping model variables to observation variables.A snow cover depletion curve (SDC) that parameterizes the relationship between SD and SCF in the CoLM is considered as the observation operator of snow DA [45]: ( where stands for the predicted SCF observation, snowdp for the CoLM simulated snow depth, and z0g for the roughness length for bare soil with a default value of 0.01 m. By defining the analysis increment = − and linearizing H(x) about x b , Equation ( 1) is rewritten as follows [48]: (3) where = − stands for the innovation vector, and = ⁄ at x b for the Jacobian matrix of the nonlinear observation operator H.The gradient of the cost function from Equation ( 3) with respect to is given as follows: (4) Next, we convert a minimum function problem into a problem of finding the solution of ∇ .The conjugate-gradient method, which offers a good tradeoff between convergence rates and computer memory requirements, is used to minimize the quadratic cost function in our work.

Deterministic Ensemble Kalman Filter
Without observation perturbations, a deterministic ensemble Kalman filter (DEnKF) has been proposed by Sakov and Oke [41].The updated equations of the DEnKF are given as follows: (5) (6) (7) (8) where ̅ = , and ̅ stand for the analysis and background ensemble means of snow depth (snowdp) and snow water equivalent (swe) states, = , for the CoLM simulated snow state, i for the ensemble index, and N for the ensemble size.K stands for the Kalman gain matrix and y° for the MODIS SCF observation.Here the MODIS SCF observation is not perturbed.H stands for an observation operator as Equation (2).P b stands for the background error covariance matrix and R for the error variance of MODIS SCF observations.= − ̅ , − ̅ , ⋯ , − ̅ stands for the background ensemble perturbation.The DEnKF maintains the ensemble spread by halving the Kalman gain in the analysis perturbation update stage.That is, ( ) where A a stands for the analysis ensemble perturbation.This formulation can be thought of as an inherent covariance inflation similar to the relaxation-to-prior algorithm of Zhang et al. [49].The covariance inflation can avoid the collapse of ensemble spread and further reduce the risk of filter divergence [41,50].
Finally, the analysis ensemble is generated by adding perturbations to the analysis ensemble mean.It is denoted by . (10)

Coupled Method
For snow data assimilation (DA), a two-way coupled DEnKF-1DVar hybrid method (DEnVar) patterned after Zhang et al. [36] is proposed.Two-way coupling re-centers the DEnKF analysis ensemble mean with the 1DVar analysis.Such re-centering may help to prevent the divergence of the DEnKF analyses so that the ensemble perturbations empirically represent the distribution of the forecast errors.Unlike the EnVar method (hybrid ensemble-variational DA) of Zhang et al. [36], the two-step DEnVar method without observation perturbations in Figure 1 includes three steps: (1) the DEnKF analysis is implemented to generate the analysis ensemble mean ̅ and the analysis error covariance , (2) the analysis ensemble mean ̅ is used as the first guess for the 1DVar and the analysis error covariance is introduced into the 1DVar hybrid cost function, and (3) the analysis ensemble mean ̅ of DEnKF is replaced with the 1DVar analysis for the next ensemble forecast.The hybrid error covariance is included in the hybrid cost function by separating the Jb in Equation (1) into two parts [28,51]: (11) where Jb1 stands for the traditional 1DVar background term as in Equation ( 1) and Jb2 for the ensemble-based term.BNMC stands for the static background error covariance calculated by the NMC (National Meteorological Center) method [52], and for the analysis error covariance valid at the time of analysis.β stands for the weighting coefficient that is used to control the respective weights of static NMC-based background and ensemble-based analysis error covariances.If β approaches 0, then the hybrid formulation is functionally equivalent to the standard 1DVar method.If β is 1, the hybrid error covariance only contains the ensemble-based analysis error covariance.If β is nonzero, then the hybrid method takes the mixed value of both static NMC and ensemble-based error covariances.

Error Evaluation Methods
To evaluate the analysis performance of these assimilation methods, we selected three indices of bias, root mean squared error (RMSE), and normalized error reduction percent (NERP) [53] as follows: (12) (13) (14) where ys,i and yo,i stand for the i-th simulated (simulation without assimilation or assimilation) and observed snow states, respectively, and n for the number of time points.The NERP index is used to measure the improved magnitudes of these assimilation methods.RMSEs and RMSEa are the root-mean square errors of simulated and analyzed snow states, respectively.NERP ranges from negative infinity to 100%, where 100% means that the assimilated results are identical to the observations.A negative NERP represents a percent of deterioration for the analyzed states as compared with the simulations without DA.

Experimental Sites
The snow DA methods were used at five sites in the Altay region of Northern Xinjiang, China (Figure 2).These five sites are located in a homogeneous land cover area.The geographic locations and elevations of five sites are given in Table 1.Influenced by the Siberian High, these five sites experience a long, cold and snowy winter.This leads to an approximate 120-day duration of snow cover from November to March of the next year [54].Only the daily in situ SD and air temperature measurements were collected at all sites from January 2004 to March 2009.Huang et al. [54] and Xu and Shu [20] have pointed out that it is feasible to compare in situ SD observations with MODIS snow cover or the CoLM-simulated SD over a variety of environmental conditions.For example, it is possible to make a comparison between in situ SD observations and the CoLM-simulated SD with a 0.01° spatial resolution.In situ SD observations from November 2008 to March 2009 were considered as validation data sets for assessing the performance of snow DA methods.The mean SDs at five sites during this period are shown in Table 1. ) )

MODIS Snow Cover Fraction
For snow DA experiments, the Terra MODIS daily snow cover fraction product at an approximately 10:00 am local overpass time (MOD10A1; Daily L3 Global 500 m Grid, Collection 5) from 2004 to 2009 was downloaded from the National Snow and Ice Data Center [55].A previous study has shown that the MODIS SCF product has a mean absolute error of less than 10% in the fractional snow cover of 0% to 100% [56].However, this SCF error of 10% may be underestimated in some areas of dense forests and slope-induced shadow [56,57], where incorrect MODIS SCF observation errors [25] exist.To be consistent with the CoLM's spatial resolution, the daily 500 m SCF product is resampled to the daily 0.01° SCF product by the nearest-neighbor resampling of the NASA's MODIS Reprojection Tool (MRT).

AMSR-E Snow Water Equivalent
The AMSR-E, launched on the NASA Aqua platform in May 2002, is a multi-frequency, dual-polarized passive microwave radiometer with frequencies of 6.9 GHz to 89 GHz.The AMSR-E is aimed to provide spatially and temporally continuous SWE products for global change analyses.The daily AMSR-E SWE product with a spatial resolution of 25 km during the period of 2004 to 2009 was downloaded from NSIDC [58].Due to homogeneous grassland land cover at five sites in the plain area of the Altay region, the AMSR-E-retrieved SWE is relatively accurate from November to February.The AMSR-E-retrieved SWE during the melting period of March may be less accurate for the reason of liquid water content [59], but overall it is relatively reliable.The AMSR-E-retrieved SWEs are downscaled to the spatial resolution of 0.01° in consistence with the resolution of the CoLM, using a nearest-neighbor resampling method in the snow DA system.The SWEs from November 2008 to March 2009 are considered as the "reference" values to assess the analysis performance of snow DA methods, and the SWEs from other times are used to calculate the MODIS SCF observation error.

Experimental Design
The China Meteorological Forcing Dataset (CMFD), with a 0.1° spatial resolution and 3-h temporal resolution provided by the Environmental and Ecological Science Data Center for West China (WestDC), includes downward short-wave solar and long-wave radiation, precipitation rate, air temperature, wind speed, specific humidity, and surface pressure (Pa).The atmospheric forcing dataset is used to drive the CoLM, which is a physically based land surface model.The CoLM can be used to forecast the energy and water balance.The transfer between energy and water of the soil-vegetation-snow-atmosphere is formalized to be a function of atmospheric forcing data, terrain elevation, and vegetation and soil characteristics.Simulated and assimilated experiments are implemented at a single grid point with a 0.01° (~1 km) spatial resolution at five sites.
A simulation in the CoLM from January 2004 to March 2009 was implemented in a continuous mode.The simulation of the first four and half years served as a spin-up of the model so that SD and SWE could reach a steady state.The previous year's simulation, without assimilation, from November 2008 to March 2009 was considered as a benchmark in testing the analysis results of different DA experiments.
We performed four snow DA experiments denoted 1DVar, DEnKF, EnVar-Beta 0.5, and DEnVar-Beta 0.5 in Table 2.In the four snow DA experiments, if the daily SCF and simulated snow states (SD and SWE) take positive values, then the Terra MODIS SCF observations are assimilated into the CoLM at the local satellite overpass time of about 10:00 am using different DA methods.Using the method of Arsenault et al. [25], the MODIS SCF observation errors (σobs) are calculated as the difference between MODIS SCF observations and SCF validation "truth".The SD data is calculated by AMSR-E SWE with a bulk snow density of 250 kg/m 3 (SWE = SD × snow density), then the SCF validation "truth" is calculated with the SD data with the observation operator, SCF = snowdp / (snowdp + 0.1).The MODIS SCF observation errors are characterized with the observation error standard deviation in the DA system.Nonzero snow observation values of AMSR-E SWE and MODIS SCF from January 2004 to October 2008 were used to calculate the MODIS SCF observation errors, and the number (Num) of nonzero snow observations is given in Table 1.The period of November 2008 to March 2009 was selected as the experimental time for assimilation and validation.The final SCF observation errors at five sites are given in Table 1.In Table 1, the observation standard errors (σobs) at five sites are approximately 23%-30%, which may be jointly caused by the uncertainty of the AMSR-E-retrieved SWE, observation operator, and snow density (e.g. a bulk snow density of 250 kg/m 3 ).Particularly in the snow melting stage, the AMSR-E-retrieved SWE may have errors introduced by liquid water content.
In the 1DVar experiment, the traditional 1DVar method of Equation ( 3) and ( 4) is used.Using the NMC method [52], the differences between 48 h and 24 h CoLM snow forecasts starting from the same time at 0900 UTC from November 2007 to March 2008 were used to estimate the static background error covariance B through the following Equation ( 15): (15) where x t stands for the true state of CoLM SD and SWE and ε b for the background error.The overbar denotes an average over time.The valid pairs (Pairs) in calculations of B during this period at five sites are given in Table 1.The cross covariance between SD and SWE is also considered in B.
In the DEnKF experiment, an ensemble size of 25 is selected as in previous studies [15,17,25].During the assimilation, the perturbations are introduced into the atmospheric forcing dataset and the snow state variables (SD and SWE) in a continuous cycle mode to represent distribution of these model input errors empirically.Following the method of De Lannoy et al. [15] and Reichle et al. [60], multiplicative perturbations with the mean of 1 and standard deviations of 0.5 and 0.1 are imposed on precipitation (P) and downward shortwave radiation (SW), respectively.Additive perturbations with the mean of 0 and standard deviations of 0.5 K and 15 W m −2 are imposed on air temperature (T) and downward longwave radiation (LW), respectively.The SD and SWE state variables are subject to multiplicative perturbations with the mean 1 and standard deviation 0.01.Moreover, the cross correlations between perturbations to state variable and forcing dataset are also calculated by using the methods of De Lannoy et al. [15] and Reichle et al. [60].These perturbations do not have spatial correlation because all DA experiments are implemented at a single grid point with a 0.01° (~1 km) spatial resolution.
To remove sampling errors that resulted from MODIS SCF observation perturbations, the EnVar-Beta 0.5 experiment, adapted from the method of Zhang et al. [36], was implemented by coupling the 1DVar and DEnKF methods.The DEnKF provides the forecast ensemble mean ( ̅ ) and background error covariance ( ) for 1DVar.The analysis ensemble mean in the DEnKF is replaced by the 1DVar analysis.The analysis perturbation in the DEnKF is updated by Equation ( 9) for the next ensemble forecast.In the EnVar-Beta 0.5 experiment, the static NMC-based background error covariance (BNMC, Equation ( 15)) and the ensemble-based background error covariance ( , Equation ( 8)) are equally weighted, that is, Unlike the EnVar-Beta 0.5 experiment, the DEnVar-Beta 0.5 experiment with the weighting coefficient of 0.5 needs a two-step process.That is, the analysis ensemble mean ( ̅ ) and analysis error covariance ( ) calculated in DEnKF are introduced into 1DVar, and the 1DVar analysis is used to update the analysis ensemble mean in DEnKF.
Furthermore, additional sensitivity tests of DEnVar to weighting coefficients and MODIS SCF observation errors (Table 2) were performed.We wrote the computer programs of these DA algorithms in FORTRAN and integrated related modules into the CoLM.
Table 2. Summary of snow DA experiments.Note that the SCF observation standard errors of these DA experiments except for the DEnVar-Beta 1.0-2 experiment are given in Table 1.

Experiments Description Simulation
No snow DA 1DVar Snow DA using 1DVar with static covariance from NMC method DEnKF Snow DA using DEnKF with ensemble size of 25 EnVar-Beta 0.5 Snow DA using a hybrid method like Zhang, Zhang and Poterjoy [36] with ensemble size of 25, weighting coefficient of 0.5 DEnVar-Beta 0.5 Snow DA using the proposed hybrid method with ensemble size of 25, weighting coefficient of 0.5 in Equation ( 11

Comparison with In situ SD Observations
In Figure 3, the CoLM running at most sites produces late snow accumulations and early snowpack melt-offs as compared with in situ SD observations, except at the Fuyun site, partly due to underestimation of albedo and SCF in the CoLM [20].At the Fuyun site, the CoLM simulation shows time series trends consistent with in situ SD observations during the snow accumulation period.However, there is still a large discrepancy between the CoLM-simulated SD and in situ observations during the melting period.Figure 3 shows that the four DA experiments have similar time series trends.
In the experiments of DEnVar-Beta 0.5, DEnKF, and EnVar-Beta 0.5, the analyzed SD values during the snow accumulation and melting periods are elevated to a value that is very close to in situ observations at the Aletai, Buerjin, and Qinghe sites.In contrast to 1DVar experiment (see Figure 3A,B,E), these three DA experiments obtain better analysis performance.However, the opposite happened at the Fuyun site.At the Fuyun site, these three DA experiments showed overestimated peak and melting snowpack as compared with the 1DVar experiment; the 1DVar analysis tended to be similar to the in situ SD time series trend as compared with other experiments in Figure 3C.The overestimation is caused by the large positive innovation − , which may be related to the predicted SCF observation underestimated by the simplified observation operator [20].It is expected to alleviate the overestimation by improving the simplified observation operator in a method similar to the SCF parameterization scheme [10,61].At the Jimunai site, the DEnVar-Beta0.5 experiment obtained the best analysis performance, but it tends to melt the snowpack too early and rapidly in Figure 3D.This is due to the presence of cloud cover, which accounts for the missing SCF observations in snow DA during this period [20].To compare the overall difference of the four DA experiments, these results were evaluated with the indices of bias, RMSE, and temporal correlation coefficient (Correlation) in Table 3.The DEnVar-Beta 0.5 experiment with the lowest bias and RMSE and the highest correlation coefficient outperformed other DA experiments at the Buerjin, Jimunai, and Qinghe sites.Among these three sites, Qinghe tended to have the lowest errors (the bias of 0.28 cm and the RMSE of 3.91 cm) and the highest correlation coefficient of 0.95 in the DEnVar-Beta0.5 experiment.The error of the analyzed SD in the DEnVar-Beta 0.5 experiment was reduced by 57.82% (NERP of 57.82%, see Figure 4) as compared with the CoLM simulations at the Qinghe site, and the errors in other DA experiments were also reduced by 18.12%, 51.13%, and 41.42%, respectively, in Figure 4.At the Aletai site, the DEnKF experiment with a bias of −0.71 cm obtained a slightly better agreement with in situ SD observations against the DEnVar-Beta 0.5 experiment with a bias of 1.85 cm in Table 3.The improved magnitude of the DEnKF and DEnVar-Beta 0.5 experiments at the Aletai site is extended upwards to 67.93%and 66.29% respectively in Figure 4.At the Fuyun site, the correlation coefficient between the 1DVar analyzed SD and in situ SD reached 0.98, with a bias of −1.41 cm and a RMSE of 3.36 cm in Table 3.These values imply that the 1DVar experiment obtains a better agreement and smaller differences with in situ SD observations as compared with other experiments.The large positive bias in Table 3 and negative NERP values in Figure 4 show overestimation of SD during the snow accumulation and melting periods in the experiments of DEnKF, EnVar-Beta 0.5, and DEnVar-Beta 0.5, especially in the DEnVar-Beta 0.5 experiments, with the largest positive bias of 8.97 cm and negative NERP of −83.36% at the Fuyun site in Figure 4.
The bias, RMSE, correlation coefficient, and NERP values significantly improved after the assimilation of available MODIS SCF observations.On average, the four DA experiments can improve SD simulations.In particular, the DEnVar-Beta 0.5 experiment obtains the best analysis performance while the 1DVar experiment obtains the worst analysis performance.The DEnKF analysis performs as well as the analysis of DEnVar-Beta 0.5 experiment and better than that of the EnVar-Beta 0.5 experiment.Since the DEnVar-Beta 0.5 experiment integrates the analysis of DEnKF and 1DVar, the DEnVar-Beta 0.5 experiment produces more snow accumulations than the DEnKF and 1DVar experiments.This is to say, the DEnVar-Beta 0.5 experiment makes significant improvements on SD simulations, especially in the case of underestimated SD simulations.The SWE state is updated with the error covariance between the SD and SWE states using a similar algorithm of the SD update.
In Figure 5, the four DA experiments significantly improve the SWE simulations at all sites, and all DA experiments show similar SWE time series trends.However, the analyzed SWEs in these DA experiments are still underestimated during the snow accumulation and melting periods in Figure 5. From November to late February, when perhaps the snow is dry, the AMSR-E-retrieved SWE has fewer errors introduced by liquid water content [62].These snow DA experiments show better agreement with the AMSR-E SWE observations during this period than during the melting period of March in Figure 5.
On the whole, the DEnVar-Beta 0.5 experiment obtained better analysis performance than other DA experiments at the Aletai, Fuyun, Jimunai, and Qinghe sites in Figure 5, with the lowest bias and RMSE and highest correlation coefficient values in Table 4.At the Aletai site, the bias between the analyzed SWE in the DEnVar-Beta 0.5 experiment and the AMSR-E SWE was reduced from −60.69 mm (the bias of the CoLM simulations) to −30.74 mm, and the RMSE was reduced from 72.17 mm to 44.78 mm.The DEnVar-Beta 0.5 experiment tended to have the largest improved magnitude of SWE at the Aletai site and the improved magnitude reached 37.95% in Figure 6.As compared with the EnVar-Beta 0.5 experiment, the DEnKF experiment showed analysis results similar to that of the DEnVar-Beta 0.5 experiment, and achieved a slightly better agreement with the AMSR-E SWE observations.The EnVar-Beta 0.5 experiment with a NERP of 21% outperformed the 1DVar experiment with a NERP of 8.45% in Figure 6, conversely, at the Buerjin site.At the Buerjin site, the 1DVar experiment offered better improvement than other DA experiments.As compared with the CoLM simulations, the SWE simulations of the 1DVar has obtained an 11.6% improved magnitude, while the largest improved magnitude among the DEnKF, EnVar-Beta0.5, and DEnVar-Beta 0.5 experiments was only 2.92%, in the EnVar-Beta 0.5 experiment in Figure 6.Although the analyzed SWE values during the snow accumulation and melting periods were improved, a large discrepancy still exists between the analyzed SWE and the AMSR-E SWE.That resulted partially from the imprecise covariance between SD and SWE states.An improved extent of the SWE state is determined by innovation = − and Kalman gain concerning the covariance between SD and SWE states.The covariance is underestimated with an inaccurate covariance.The underestimated covariance may decrease the Kalman gain, and further lead to a sub-optimal assimilation.

Sensitivity to Weighting Coefficients
Here we performed the sensitivity test of the weighting coefficient β to the proposed DEnVar method.The weighting coefficient β, as given in Equation (11), controls the respective weights of the static NMC and ensemble-based analysis error covariances.The weighting coefficient is set to 0.8 in the DEnVar.This explains that 20% of the covariance comes from the static NMC estimate and 80% from the ensemble error statistic, as described in the DEnVar-Beta 0.8 experiment.The weighting coefficient was set to 1.This explains that the hybrid error covariance only contains the ensemble-based analysis error covariance (as in the DEnVar-Beta 1.0 experiment).
Figure 7 shows the time series of the CoLM snow states from the DEnVar experiment, with β of 0.5 in the DEnVar-Beta 0.5 experiment, 0.8 in the DEnVar-Beta 0.8 experiment, and 1.0 in the DEnVar-Beta 1.0 experiment at the Aletai and Jimunai sites.The time series graphics at other sites are not shown given their similar analysis performance.In Figure 7, the three DA experiments show similar time series trends; the DEnVar-Beta 1.0 experiment obtained a slightly closer fit to snow observations than the DEnVar-Beta 0.5 and DEnVar-Beta 0.8 experiments.For the SD state, the three DA experiments have similar NERP values in Figure 8A, while the DEnVar-Beta 1.0 experiment has lower NERP values at the Aletai, Buerjin, and Fuyun sites, and higher NERP values at the Jimunai and Qinghe sites than in the other DA experiments.In addition, for the SWE state, the DEnVar-Beta 1.0 experiment with slightly higher NERP values achieved a slightly better analysis performance than the DEnVar-Beta 0.5 and DEnVar-Beta 0.8 experiments at all sites (Figure 8B).Therefore, it is concluded that the proposed DEnVar method is not highly sensitive to the weighting coefficient β.However, it prefers to give greater weight to the ensemble-based analysis error covariance since the DEnVar experiment with a greater weighting coefficient produces more snow accumulations and achieves better correlation with snow observations, especially in the case of SD-and SWE-underestimated simulations.

Sensitivity to Observation Error
Given that MODIS SCF products have approximately 10% error [54,63], some previous studies have assigned a 10% SCF observation error to MODIS SCF DA [16,17].To identify the influence of this 10% SCF observation error on the proposed DEnVar DA method, a sensitivity test extending the DEnVar-Beta 1.0 experiment, a DEnVar-Beta 1.0-2 experiment, was performed.
In Figure 9, the DEnVar-Beta 1.0-2 experiment showed an overestimated peak snowpack as compared with the DEnVar-Beta 1.0 experiment, which produced early snow accumulations and late snowpack melt-offs at all sites except for Buerjin.The DEnVar-Beta 1.0-2 experiment tended to overestimate the SD simulations at all sites, leading to large bias and RMSE values in Figure 10A.The Fuyun site has the largest bias of 20.06 cm and the largest RMSE of 22.10 cm.The lowest bias of 1.33 cm and RMSE of 5.19 cm occurred at the Buerjin site (Figure 10A).Figure 11 shows that the analyzed SWE values in the peak and melting stages were greatly increased.This led to overestimation of the SWE simulations at the Aletai and Jimunai sites, with positive biases of 11.39 mm and 18.99 mm, respectively (Figure 10B).Although the analyzed SWE is still underestimated at the Fuyun and Qinghe sites (as seen in the negative bias in Figure 10B), the time series trend from the DEnVar-Beta 1.0-2 experiment was closer to that of the AMSR-E-retrieved SWE observations than that from the DEnVar-Beta 1.0 experiment in Figure 11.On the whole, large changes in the SD and SWE states were produced by the DEnVar-Beta 1.0-2 experiment, with a 10% SCF observation error.In sum, the proposed DEnVar method is highly sensitive to observation error.Moreover, a reasonable MODIS SCF observation error is quite important for snow DA.A triple collocation technique [64] may be used to estimate the MODIS SCF observation error under the assumption that the three snow-related datasets are independent of each other.

Conclusions
In this study, four different data assimilation (DA) experiments assimilating MODIS SCF observations into the CoLM were performed using ensemble-based, variational, and hybrid methods, DEnKF, 1DVar, EnVar (a hybrid ensemble-variational DA system), and DEnVar (DEnKF-1DVar).EnVar and DEnVar are coupled DA methods, which are comparable to the method of Zhang et al. [36] and couple 1DVar and DEnKF methods.Unlike EnVar, the DEnVar method needs a two-step process.That is, the analysis ensemble mean and the analysis error covariance of DEnKF are introduced into the 1DVar hybrid cost function, and the analysis mean of DEnKF is replaced by the 1DVar analysis.
DA experiments were performed at five sites in the Altay region of China from November 2008 to March 2009.The analyzed SD and SWE were validated against in situ SD and AMSR-E SWE observations.Our results indicate that some improvements in the four DA experiments occurred at most sites when the available MODIS SCF is assimilated.It was shown that, with the same observation errors, the 1DVar of static NMC-based background error covariance performed worse than other DA methods of ensemble-based background error covariance.This indicates that the static NMC-based background error covariance magnitude from 1DVar is less than the ensemble-based background error covariance from the DEnKF.A possible reason is that a larger analysis error covariance = − + 1 4 of the DEnKF method is calculated by the equations = and = − 1 2 .This analysis error covariance explains that the DEnKF can avoid underestimating the analysis error covariance [41].
As compared with the stochastic EnKF, an additional term 1 4 in the DEnKF is interpreted as an inherent covariance inflation, which can effectively maintain the ensemble spread and avoid filter divergence.Furthermore, the DEnKF without observation perturbations does not introduce extra observation perturbation errors, which may render the filter sub-optimal [40].The analysis error covariance and ensemble mean ̅ of DEnKF are used in the hybrid 1DVar cost function (DEnVar-Beta 0.5).The 1DVar analysis in DEnVar-Beta 0.5 experiment is better than that in the EnVar-Beta 0.5 experiment, which uses the forecast error covariance and ensemble mean ̅ .As a result, our hybrid DEnVar method obtains better analysis performance than the EnVar method of Zhang et al. [36].The analysis performance of the two-step DEnVar method is better than that of the DEnKF method.
Although our hybrid DA method can improve the simulations of SD and SWE to some extent, large errors can still be found in snow simulations of the CoLM.The large errors may be partly due to the scale mismatch between the CoLM model and input datasets.In practice, the atmospheric forcing data with a 0.1° spatial resolution is effectively implemented to drive the CoLM with a 0.01° spatial resolution in simulations and assimilations of SD and SWE.However, this atmospheric forcing data at a larger scale may introduce some errors leading to a bias in simulations of SD and SWE.It is due to this that atmospheric forcing data may not satisfy the requirements of simulations and assimilations on the CoLM's scale, as well as the assimilated MODIS SCF with a 500 m spatial resolution.Furthermore, to be consistent with the CoLM's scale, the AMSR-E SWE data is downscaled to a spatial resolution of 0.01° using a nearest neighbor resampling.The value of the scaled SWE data does not differ from the original SWE data value with a 25 km spatial resolution.Thus, the scaled SWE data used to estimate the observation errors and validate the assimilated results may also have introduced some errors into our DA experiments.
Lastly, we tested the sensitivity of the weighting coefficient of hybrid background error covariance and the SCF observation error to the DEnVar method, respectively.It is found that the proposed DEnVar method is slightly sensitive to the weighting coefficient and is highly sensitive to observation errors.The DEnVar-Beta 1.0 experiment with the ensemble-based analysis error covariance obtains better analysis performance than the DEnVar-Beta 0.5 and DEnVar-Beta 0.8 experiments with the hybrid error covariance, because the element values of background error covariance calculated with the hybrid error covariance are decreased in the 1DVar cost function.It is thus feasible to use only the ensemble-based error covariance from the DEnKF in the hybrid DA system.

Figure 1 .
Figure 1.A hybrid analysis method of DEnVar which couples the DEnKF and the 1DVar.

Figure 2 .
Figure 2. The geographic locations of five sites in the Altay region of Northern Xinjiang, China on the land cover map of approximately 1 km spatial resolution [20].

Figure 3 .
Figure 3.Comparison of SD time series from in situ observations and the assimilated results from different DA experiments at the (A) Aletai, (B) Buerjin, (C) Fuyun, (D) Jimunai, and (E) Qinghe sites.

Figure 5 .
Figure 5.Comparison of the SWE time series from AMSR-E retrieved SWE and the assimilated results from different DA experiments at the (A) Aletai, (B) Buerjin, (C) Fuyun, (D) Jimunai, and (E) Qinghe sites.

Figure 7 .
Figure 7.The time series comparison of CoLM snow states (Left: SD, Right: SWE) from snow observations and the assimilated results in DEnVar-Beta 0.5, DEnVar-Beta 0.8, and DEnVar-Beta 1.0 sensitivity tests at the Aletai (A,B) and Jimunai sites (C,D).

Figure 9 .
Figure 9.The SD time series of in situ SD observations and the analyzed SD in the DEnVar-Beta 1.0 and DEnVar-Beta 1.0-2 sensitivity tests at the (A) Aletai, (B) Buerjin, (C) Fuyun, (D) Jimunai, and (E) Qinghe sites.The observation standard errors (σobs) at five sites in the DEnVar-Beta 1.0 experiment are given in Table1.The DEnVar-Beta 1.0-2 experiment is the same as the DEnVar-Beta 1.0 experiment, but the MODIS SCF observation error is set to 10%.

Figure 10 .
Figure 10.The bias and RMSE of the analyzed SD (A) and SWE (B) in DEnVar-Beta 1.0-2 experiment with the MODIS SCF observation error of 10%.

Table 1 .
Characteristics of five sites in the Altay region of China.The mean SD is calculated from daily in situ SD observations from November 2008 to March 2009, and the mean SWE is calculated by daily AMSR-E retrieved SWE from November 2008 to March 2009.The term "Pairs" implies the valid pairs of 48-h and 24-h CoLM snow forecasts from November 2007 to March 2008.The number (Num) of nonzero snow observations was calculated during the period of January 2004 to October 2008.The observation standard error (σobs) is calculated by fitting MODIS SCF observations to AMSR-E SWE-based SCF observations.

Table 3 .
Summary of the analyzed SD simulation statistics in 1DVar, DEnKF, EnVar-Beta 0.5, and DEnVar-Beta 0.5 DA experiments at five sites from November 2008 to March 2009.

Table 4 .
Summary of the analyzed SWE statistics.Cf.Table3.