Improving the Informational Value of MODIS Fractional Snow Cover Area Using Fuzzy Logic Based Ensemble Smoother Data Assimilation Frameworks

Abstract

Remote sensing fractional snow cover area (fSCA) has been increasingly used to get an improved estimate of the spatiotemporal distribution of snow water equivalent (SWE) through reanalysis using different data assimilation (DA) schemes. Although the effective assimilation period of fSCA is well recognized in previous studies, little attention has been given to explicitly account for the relative significance of measurements in constraining model parameters and states. Timing of the more informative period varies both spatially and temporally in response to various climatic and physiographic factors. Here we use an automatic detection approach to locate the critical points in the time axis where the mean snow cover changes and where the melt-out period starts. The assimilation period was partitioned into three timing windows based on these critical points. A fuzzy coefficient was introduced in two ensemble-based DA schemes to take into account for the variability in informational value of fSCA observations with time. One of the DA schemes used in this study was the particle batch smoother (Pbs). The main challenge in Pbs and other Bayesian-based DA schemes is, that most of the weights are carried by few ensemble members. Thus, we also considered an alternative DA scheme based on the limits of acceptability concept (LoA) and certain hydrologic signatures and it has yielded an encouraging result. An improved estimate of SWE was also obtained in most of the analysis years as a result of introducing the fuzzy coefficients in both DA schemes. The most significant improvement was obtained in the correlation coefficient between the predicted and observed SWE values (site-averaged); with an increase by 8% and 16% after introducing the fuzzy coefficient in Pbs and LoA, respectively.

1. Introduction

Snow has important relevance for society spanning from environmental and economical to recreational and aesthetic values. Environmentally snow influences the climate and ecology of an area through its effect on solar energy absorption and by constraining land cover characteristics [1,2]. Snowmelt results in significant decrease of surface albedo and releases the water stored in seasonal snowpack, providing water for humans and the environment [3,4]. Water managers in hydropower and water supply sectors consider snow storage in mountain areas as a natural reservoir [5]. During the past decades, hydrological models have been used to assess the amount and spatiotemporal distribution of this important resource [6,7]. However, the modeling process in general and model calibration in particular is prone to various sources of uncertainty [8]. Further, the objective functions commonly used for parameter inference using inverse modeling are based on aggregate measures of model performance and usually yield to poorly constrained model parameters [9]. A combined approach of model calibration followed by data assimilation is one of the strategies commonly followed to deal with model uncertainty [10,11]. Data assimilation (DA) is the means to obtain an estimate of the true state through use of independent observations and the prior knowledge (model state) with appropriate uncertainty modeling [9,12].

Ensemble-based filtering and smoothing techniques are commonly used in solving DA problems [9]. The filtering techniques involve updating of state variables and/or parameter values at each observation time with subsequent production of error statistics that can be cycled to future analysis times; while in smoothing all observations are assimilated in a single step [4,13]. Some of the filtering- (sequential) based DA schemes commonly used in hydrological data assimilation include the Ensemble Kalman Filter [14] and the Particle filter [15], while some examples from the smoother-based DA schemes include the ensemble batch smoother [16] as well as the particle batch smoother (Pbs) [17]. The smoother-based DA schemes have the advantage of low computational cost coupled with the flexibility in running independent of the forward model [4,18,19].

In recent years, remote sensing data in general and optical fractional snow cover area (fSCA) in particular has been increasingly used in hydrologic modeling for constraining model parameters through inverse modeling and various DA schemes [20]. This can be attributed to the readily available observations of this data for large areas and at a relatively high temporal resolution [21,22]. In DA, fSCA has been used from a simple rule-based (direct-insertion) approach [23] to advanced filtering and smoothing techniques [24,25]. The collection of fSCA measurements and their temporal evolution during the ablation period reflect the peak snow water equivalent (SWE) better than a single fSCA measurement does. Hence, the batch assimilation techniques are preferable over filtering techniques when the assimilated observation is fSCA in a snow depletion curve model [17]. The Pbs and particle filter DA schemes have been gaining interest as fSCA assimilation schemes in snow models due to their distribution free likelihood and their capability to estimate the SWE state vector directly at a relatively low computational cost [17,26]. However, these schemes assume crisp and equal informational value for all observations of the assimilated data. Further, in these schemes, most of the weights are assigned to one or very few ensemble members and this may lead to degeneration of the statistical information in the ensembles [27,28]. An alternative approach that ensures fair distribution of weights among ensemble members is the limits of acceptability (LoA) [29]. This approach provides flexibility in choosing different shapes for representing the error within specified error bounds. It further allows the incorporation of different hydrological signatures in the evaluation criteria. As such, it is more process-based as compared to the simple residual-based objective functions since it transforms the observed and simulated data into hydrologically relevant signatures [30].

In this study, the informational value of fSCA measurements is considered to be fuzzy, carrying a considerable amount of uncertainty. In previous studies different approaches have been used to increase the amount of information that can be extracted from time series data during parameter inference using inverse modeling [30]. Some of these measures include: breaking up the time series data into different periods, for example, the different components of a hydrograph [31] as well as transforming the time series data to emphasise on some aspects of the system response [32]. In several snow related DA studies all measurements are assumed to carry the same amount of information in constraining the model states and parameters. However, the amount of information that can be extracted from remote sensing data varies depending on several factors including the presence of certain physical limitations such as clouds and lack of sufficient light [33]. For example, high latitudes are characterized by polar darkness and frequent cloud cover and errors emanating from such sources can reduce the informational value of MODIS fSCA product for snow assimilation [34,35].

The informational value of MODIS fSCA was divided into three timing windows; and an automatic change point (CP) detection scheme was employed to identify the spatially and temporally variable location of the CP in the time-axis. CP detection deals with the estimation of the point at which the statistical properties of a sequence of observations changes [36]. In hydrology, CP detection methods have been applied in various studies including the assessment of non-stationarity in time series data [37,38] as well as in the assessment of flow-solute relationship in water quality studies [39]. Nevertheless, their application has been limited when it comes to the analysis of remote sensing snow cover data. Although, snowmelt at high latitudes is a rapid process, the start date of an ablation period may vary by several weeks from year to year in response to such factors as the amount of accumulated SWE before the onset of snowmelt and the availability of melt energy [40]. Spatially, the onset of snow melt varies from one grid-cell to another depending on variability in physiographic characteristics of the area in addition to the above mentioned climatic factors. During the accumulation period, spatial variability of snow results from snow-canopy interaction in forested areas, snow redistribution by wind, and orographic effects [41]. During the ablation period, snow melts from preferential locations, yielding heterogeneous patterns [42].

The main goal of this study is to get an improved estimate of SWE during the maximum accumulation period using fuzzy logic-based ensemble fSCA assimilation schemes and by incorporating uncertainties in selected model forcing and parameters. The first objective is to assess whether the assumption of variable informational value of fSCA observations depending on their location in different timing windows can reduce the uncertainty in SWE reanalysis. We introduce a novel approach that accommodates this concept by incorporating a fuzzy coefficient in the ensemble-based SWE reanalysis schemes. The second objective is to evaluate the viability of statistical change point detection methods to locate the critical points that constitute the timing windows. The third objective is to compare effect of the assumed variability in informational value of fSCA observations on relative performance of the Pbs and the LoA-based DA schemes. As of our best knowledge, the LoA approach was not used before for DA purposes; and this study will assess the feasibility of this approach as an ensemble smoother DA scheme.

This paper is organized as follows. Section 2 presents the snow model followed by description of the study area as well as the static and time-series data used in this study. This section also presents setup of the different DA schemes and methodologies followed to detect the critical dates in the ablation period. Results from the DA schemes in comparison to the observed values as well as sensitivity of the evaluation metrics to location of the change points in the ablation period are presented in Section 3. Finally, certain points from the methodology and results sections are discussed in Section 4 in light of previous studies; and relevant conclusions are drawn in Section 5.

2. Methods and Materials

2.1. The Hydrological Model

The forward model used in this study, PT_GS_K, is a conceptual hydrological model (method stack) within the Statkraft Hydrological Forecasting Toolbox (Shyft, https://github.com/statkraft/shyft). This model and the modeling framework were described in previous studies [43,44] and this section presents main features of the model with focus on its snow method. PT_GS_K encompasses the Priestley–Taylor (PT) method [45] for estimating potential evaporation; a quasi-physical-based snow method (GS) for snow melt, sub-grid snow distribution and mass balance calculations; and a simple storage-discharge function (K) [46,47] for catchment response calculation.

PT_GS_K operates on a single snow layer and it can be forced by hourly or daily meteorological inputs of precipitation, temperature, radiation, relative humidity, and wind speed. The Model uses a Bayesian Kriging approach to distribute the point temperature data over the domain, while for the other forcing variables it uses an inverse distance weighting approach. The model generates streamflow, fSCA, and SWE as output variables.

The potential evaporation calculation in the PT method requires net radiation, the slope of saturated vapor pressure, the Priestley–Taylor parameter, the psychometric constant, and the latent heat of vaporization [48]. Actual evapotranspiration is assumed to take place only from snow free areas and it is estimated as a function of potential evapotranspiration and a scaling factor.

The energy balance calculation in GS method follows a similar approach as used by DeWalle and Rango [49] (Equation (1)). The precipitation falling in each grid-cell is classified as solid or liquid precipitation depending on a threshold temperature (tx) and on the local temperature values. The available net snow melt energy ($NE$) is the sum effect of different energy sources and sinks in the system. These include: incoming short wave radiation ($SW$), incoming ($ILR$) and outgoing ($OLR$) long wave radiation, the subsurface energy flux ($G$) as well as the turbulent sensible ($SH$) and latent ($LH$) energy fluxes emanating from rainfall, solar radiation, wind turbulence and other sources.

$$NE=SW\ast \left(1-albedo\right)+ILR+OLR+SH+LH+G$$

(1)

Among other factors, the energy contribution from short wave radiation depends on snow albedo. For a given time step (t), the snow albedo of each grid cell depends on the minimum (${\alpha }_{min}$) and maximum (${\alpha }_{max}$) albedo values as well as on air temperature ($T_a$) (Equation (2)). In this method the decay rates of albedo due to snow aging as a function of temperature, i.e., the fast (${\alpha }_{fdr}$) and slow (${\alpha }_{sdr}$) albedo decay rates corresponding to temperature conditions above and below 0 °C respectively, are parameterized [50].

$${\alpha }_t=\{\begin{array}{c}{\alpha }_{min}+\left({\alpha }_{t-1}-{\alpha }_{min}\right)\ast \left(\frac{1}{2^{1/{\alpha }_{fdr}}}\right),T_a>0 °\mathrm{C}\\ {\alpha }_{t-1}+\left({\alpha }_{max}-{\alpha }_{min}\right)\ast \left(\frac{1}{2\left({\alpha }_{sdr}\right)}\right),T_a\le 0 °\mathrm{C}\end{array}$$

(2)

The incoming and outgoing longwave radiations are estimated based on the Stephan–Boltzmann law. Turbulent heat contribution is the sum of latent and sensible heat. Wind turbulence is linearly related to wind speed using a wind constant and wind scale from the intercept and slope of the linear function, respectively [50]. The subsurface energy flux is a function of snow surface layer and the snow surface temperature.

The sub-grid snow distribution in each grid cell is described by a three parameter Gamma probability distribution snow depletion curve (SDC) [51,52]. This SDC is based on the assumption that certain proportion of a grid cell area ($y_0$) remains snow-free throughout the winter season such as in steep slopes and ridges due to wind erosion or avalanches [53]. As such, the two parameter Gamma distribution function, characterized by the average amount of snow $m$ (mm) and the sub-grid snow coefficient of variation ($C{V}_s$) at the onset of the melt season (Equation (3)), is applied only to the remaining portion of a grid cell area to estimate the fraction of the initially snow covered area where snow has disappeared ($y_t^{\prime }$). Following this formulation, the bare ground fraction at each time step ($y_t)$ is estimated using Equation (4).

$$y_t^{\prime }=\underset{0}{\overset{\lambda \left(t\right)}{{\displaystyle \int }}}f\left(x;m,C{V}_s\right)dx=\gamma \left(C{V}_s^{-2},\frac{\lambda \left(t\right)}{m·C{V}_s^2}\right)$$

(3)

$$y_t=y_0+\left(1-y_0\right)\ast y_t^{\prime }$$

(4)

where, $f$ denotes the Gamma probability density function and $\gamma$ is the incomplete Gamma function with shape $C{V}_s^{-2}$ and scale $m·C{V}_s^2$. The variables $x$ and $\lambda \left(t\right)$ respectively refer to point snow storage and the accumulated melt depth (mm) at time t since the onset of the melt season.

Although, the catchment response function (K) was used for estimating the regional model parameter values through calibration against observed streamflow data, this method was not directly employed in this study.

2.2. Study Area and Data

The study sites are located in Nea-catchment, Sør-Trøndelag County, Norway (Figure 1). The Nea-catchment covers a total area of 703 km² and its geographical location ranges from 11.67390° to 12.46273° E and from 62.77916° to 63.20405° N. Altitude of the study sites ranges from around 700 masl at site-1 to above 1200 masl at site-7. The fact that the model was previously setup for the study area plus the availability of long term SWE data, i.e., for nine years and in nine different sites, make the Nea-catchment an appealing study area to apply the DA schemes. The observed SWE data cover the hydrological years 2008–2016; and the datasets in these years were used in validating the DA schemes. The highest and lowest average catchment peak SWE values estimated from the nine sites during these hydrological years were 717 and 331 mm respectively in years 2012 and 2014. Mean annual precipitation for the nine hydrological years was 1279 mm. The highest and lowest average daily temperature values for this period were 20 and −23 °C, respectively.

The SWE reanalysis was conducted on 33 grid cells located at nine sites (

Table 1. Summary of site-averaged annual peak observed snow water equivalent (SWE) (m w.e.) as well as average site elevation (masl) and number of model grid-cells in each site of the study area.

Site	SWE (m w.e.)			No. Of Grid-Cells	Avg. Elev. (masl)
	Min.	Max.	Average
1	0.269	0.558	0.400	4	746
2	0.411	0.931	0.606	9	907
3	0.159	0.764	0.487	2	873
4	0.299	0.659	0.454	5	958
5	0.398	0.818	0.603	3	1064
6	0.34	0.668	0.430	2	988
7	0.256	0.391	0.340	3	1251
8	0.258	0.383	0.316	1	976
9	0.152	0.549	0.274	4	993

). The location of these grid-cells was constrained by the availability of snow course measurements in the catchment. Daily time series data of most meteorological forcing variables, i.e., temperature, precipitation, radiation, and wind speed were obtained from Statkraft [54] as point measurement, while daily gridded relative humidity data were retrieved from ERA-interim [55]. Different physiographic and observational datasets were also used to generate the simulation results and in subsequent reanalysis using the DA schemes as well as during evaluation of the results. The model was setup over each grid cells of 1 km²; requiring average elevation, grid cell total area, as well as the areal fraction of forest cover. Data for these physiographic variables were retrieved from two sources: the land cover data from Copernicus land monitoring service [56] and the 10m digital elevation model (10m DEM) from the Norwegian mapping authority (Kartverket.no). The land cover data show that, the catchment is mainly dominated by moors, heathland and some sparse vegetation; and only limited part of the catchment is forest covered (3%). The grid cells considered in the DA study are located on the open areas of the catchment. The median SWE data of each grid cell were derived from radar measurements provided by Statkraft [54]. The data were collected once a year in the month of April, where accumulated snow storage approximately attains its peak magnitude. The radar measurements roughly follow the same course each year in the nine representative sites of the catchment. Daily fraction of snow cover area (fSCA) was retrieved from NASA MODIS snow cover products [57]. Frequent cloud cover is one of the major challenges when using MODIS and other optical remote sensing data in Norway. A composite dataset was thus formed using data retrieved from the Aqua and Terra satellites, MYD10A1 and MOD10A1 products respectively in order to minimize the effect of obstructions and misclassification errors emanating from clouds and other sources. The days with available snow cover data were used without filtering based on further criteria such as the remaining clouds. A direct assessment of the bias associated with this product was not done in this study as independent snow cover observations were not available for the study domain. The estimated annual average bias of these MODIS fSCA products for Northern Hemisphere is approximately 8 % in the absence of cloud [58] and in forest dominated areas it may reach up to 15% [59]. This MODIS product was reported to have a relatively lower accuracy as compared to other more recent products, e.g., the MODSCAG algorithm [1] mainly due to its limited use of the available information in a given scene [60].

2.3. The Data Assimilation Strategy

This section presents the perturbed forcing and model parameters followed by description of the two main DA schemes employed in this study, i.e., Pbs and LoA as well as other two similar versions of these schemes that take into account for variability in informational value of fSCA. The ensemble-based DA schemes involve two steps, namely prediction and updating. The prediction step provides the prior estimates of SWE using the dynamic (meteorological) and static (physiographic) forcings. In the updating step, the DA schemes are applied to constrain the prior SWE estimates using the fSCA observations and thereby yielding a posterior SWE estimates. Each grid cell is updated independently.

2.3.1. Perturbed Forcing and Model Parameters

In this study, one meteorological forcing and two model parameters are randomly perturbed to produce an ensemble of 100 model realizations. Precipitation was perturbed since meteorological forcing, in general and precipitation in particular, are the major sources of uncertainty in simulating snow storage and depletion processes [61,62]. Considering uncertainty in precipitation provides the means to control the amount of accumulated SWE in the absence of a viable spatial snow redistribution mechanism [17]. The perturbed model parameters are sub-grid snow coefficient of variation (CV_s) and initial bare ground fraction ($y_0$). Snow coefficient of variation affects the SDC through its impact on the rate of snow depletion. In snow models with the sub-grid snow distribution component parameterized using statistical probability distribution function, the rate of snow depletion is inversely related to CV_s [41]. Thus, when the sub-grid snow cover is more variable, the ground gets exposed earlier than when it is uniformly distributed [20]. The initial bare ground fraction plays an important role in the SDC especially during the onset of snow melt when the fraction of initial snow covered area where snow has disappeared ($y_t^{\prime },\mathrm{Equation}\left(3\right)$) is at its lowest [52].

Both log-normal and logit-normal probability distributions are used to perturb the model forcing and parameters. A log-normal distribution was assumed for precipitation, while for CV_s and $y_0$ a logit-normal distribution was assumed. A similar approach to previous SWE reanalysis studies [17,28] was followed in perturbing precipitation and the model parameters. Instead of directly perturbing the precipitation input at each time step ($P_t$), a perturbation parameter ($b_j$) was randomly generated from a log-normal distribution and the generated, $b_j$ values are used as multiplicative biases for their respective ensemble member, j (Equation (5)). Since the correction for wind induced precipitation under-catch was previously done, it was assumed that the average ensemble, $b$, value is unbiased with a value of unity. The snow coefficient of variation (CV_s) estimated during calibration of the model against observed streamflow was adopted as the prior mean value. The mean value of $y_0$ was set to 0.04 based on previous studies in Norway [52]; and its variability was kept to a narrow range, since this parameter can easily compensate for differences between the observed and simulated fSCA values [52].

$$P_{j,t}=b_j{P}_t$$

(5)

A standard deviation (${\sigma }_{fSCA}$) of 15% was used to represent the measurement error in fSCA. This value is set based on previous studies focused on MODIS measurement errors [59].

Table 2. Main input variables used in the SWE reanalysis using the different data assimilation (DA) schemes. fSCA: fractional snow cover area.

Variable	Value	Description
bmean	1.0	mean precipitation bias
bCV	1.0	coefficient of variation of precipitation bias
CVs,mean	0.68	mean snow coefficient of variation (calibrated value)
CVs,min	0.0	minimum snow coefficient of variation
CVs,max	1.5	maximum snow coefficient of variation
CVs,std	0.55	standard deviation of the snow coefficient of variation
y0,mean	0.04	mean initial bare ground fraction
y0,min	0.0	minimum initial bare ground fraction
y0,max	0.1	maximum initial bare ground fraction
y0,std	0.05	standard deviation of initial bare ground fraction
${\mathit{\sigma }}_{\mathit{f}\mathit{S}\mathit{C}\mathit{A}}$	15%	measurement error standard deviation of fSCA
Ne	100	ensemble size
No	-	number of fSCA observations

presents the perturbation parameters and their assumed prior values.

2.3.2. fSCA Assimilation using Particle Batch Smoother

Previous studies indicate that the batch smoother is better suited than the filtering approach when dealing with SWE reanalysis through assimilation of fSCA observations [17,63]. A particle batch smoother (Pbs) uses an ensemble of independent randomly drawn Monte Carlo samples (particles) and estimates the posterior weight of each particle. The procedure involves generating ensemble of model realizations by running the model over the full seasonal cycle followed by the assimilation of all fSCA measurements at once. The Bayes theorem forms the basis for estimating the updated relative importance (weight, w_j) of each ensemble member (Equation (6)).

$$w_j=\frac{p\left(y|{\widehat{\mathit{X}}}_j\right)p\left({\widehat{\mathit{X}}}_j\right)}{{{\displaystyle \sum }}_{j=1}^{N_e}p\left(y|{\widehat{\mathit{X}}}_j\right)p\left({\widehat{\mathit{X}}}_j\right)}$$

(6)

where $p\left(y|{\widehat{\mathit{X}}}_j\right)$ refers to the likelihood of the observations (y) given the state value of the jth ensemble member (${\widehat{\mathit{X}}}_j$); and $p\left({\widehat{\mathit{X}}}_j\right)$ denotes prior values.

In Pbs, all ensemble members are implicitly assigned equal prior weights of 1/Ne [17]. When using the Gausian likelihood, this yields [17,28]:

$$w_j=\frac{exp\left[-0.5{\left(y-{\widehat{\mathit{Y}}}_j\right)}^T{\mathit{R}}^{-1}\left(y-{\widehat{\mathit{Y}}}_j\right)\right]}{{{\displaystyle \sum }}_{j=1}^{N_e}exp\left[-0.5{\left(y-{\widehat{\mathit{Y}}}_j\right)}^T{\mathit{R}}^{-1}\left(y-{\widehat{\mathit{Y}}}_j\right)\right]}$$

(7)

where $y$ and ${\widehat{\mathit{Y}}}_j$ respectively refer to $N_o\times 1$ vector of perturbed fSCA observations and predicted fSCA for the jth particle in the $N_o\times N_e$ matrix ($\widehat{\mathit{Y}}$). $\mathit{R}$ denotes $N_o\times N_o$ diagonal observation error covariance matrix.

Pbs assumes that both the prior and posterior states of the particles remain the same [17]. It updates the particle weights in such a way that particles with their predictions closer to the observations get higher weight than those with farther from the observations. The median and prediction quantiles are estimated from the cumulative distribution of the sorted SWE values and their associated weights.

2.3.3. fSCA Assimilation Using the LoA

We use two hydrologic signatures to constrain the perturbation parameters by the LoA scheme through weighting the relative importance of the particles with a fuzzy rule. In a fuzzy rule, the information contained in a fuzzy set is described by its membership function [64]. The first signature is based on the ability of each ensemble member to reproduce the observations within pre-defined error bounds of fSCA observations. The same value of measurement error as used in Pbs was also adopted here, i.e., 15% (Table 2). The deviation between the observed and simulated values was converted into a normalized criterion using a fuzzy rule-based scoring function. A triangular membership function was assumed with its support representing the uncertainty in MODIS fSCA observations and the pointed core representing a perfect match between the observed and predicted fSCA (Equation (8)). The weight of each ensemble member with respect to this signature (${\mathit{\text{w\_e}}}_j$) is calculated as the membership grade of the prediction error, summed over all observations (Equation (9)).

$${\mathsf{\mu }}_{fSCA}\left(e\right)=\{\begin{array}{c}0,e\le l\\ \frac{e-l}{m-l},l<e\le m\\ \frac{u-e}{u-m},m<e<u\\ 0,e\ge u\end{array}$$

(8)

$${\mathit{\text{w\_e}}}_j={\displaystyle \sum }_{i=1}^{N_o}{\mathsf{\mu }}_{fSCA}\left(e\right)$$

(9)

where ${\mathsf{\mu }}_{fSCA}\left(e\right)$ is the membership grade of each prediction error $\left(e\right)$ corresponding to the observed fSCA value $i$; $m$ is the point in the support with perfect match between the observed and predicted fSCA values. The variables $l$ and $u$ respectively refer to the lower and upper bounds of the fSCA observation error.

The second hydrologic signature is based on the degree of persistency of an ensemble member in reproducing the observations. This is estimated by the percentage of observations ($p$) where the predicted values fall within the MODIS fSCA error bounds. Ensemble members whose prediction failed to reproduce at least 50% of the observations within the error bounds are considered as non-persistent and assigned a membership grade of zero. The degree of membership linearly increases for those that satisfy the condition from 50% to 95% of the observations; and attains its maximum value after 95%. Thus, ensemble members that are able to reproduce 95% to 100% of the observations within the error bounds are considered as the most persistent and accordingly assigned a membership grade of 1.0. The weight of each ensemble member with respect to its persistence in reproducing the observations (${\mathit{\text{w\_p}}}_j$) is equal to the degree of membership (${\mathsf{\mu }}_N\left(p\right),$ Equation 10). The combined weight ($w_j$) resulting from the two membership functions is determined by calculating the product of ${\mathit{\text{w\_e}}}_j$ and ${\mathit{\text{w\_p}}}_j$ followed by normalizing the result in such a way that all particle weights would sum up to 1.

$${\mathsf{\mu }}_N\left(p\right)=\{\begin{array}{c}0,p<l\\ \frac{p-l}{m-l},l\le p\le m\\ 1,p\ge m\end{array}$$

(10)

where ${\mathsf{\mu }}_N\left(p\right)$ is the membership grade of each ensemble member with respect to its persistency (p) in reproducing the observations; $l$ represents the minimum bound of the support (p = 50%) and $m$ refers to the point in the support where the membership grade starts to attain its maximum value (p = 95%).

2.3.4. DA Schemes with an Account for Fuzzy Informational Value of fSCA

Here, the ablation period is divided into three timing windows with respect to the assumed relative informational value of fSCA observations (Figure 2). This assumption is based on the concept that information is gained from an observation only if there is uncertainty about it; and an event that occurs with a high probability conveys less information and vice versa [65]. The plot shows temporal distribution of fSCA for a sample grid-cell in the study area. The first timing window (W-1) is characterized by some snow melt with a subsequent decrease in SWE while the fSCA remains close to its maximum value. The fSCA observations from this window are characterized by low signal-to-noise ratio. In the second window (W-2), fSCA starts to drop at a faster rate and this window is characterized by high signal-to-noise ratio with a strong decreasing trend with time. The third window (W-3) corresponds to the complete melt-out period; and with the exception of few fSCA observations from sporadic snow events, fSCA is expected to remain null throughout this period. The fSCA observations in W-1 and W-3 are thus relatively certain, conveying less information as compared to observations in W-2.

The informational value of MODIS fSCA observation is considered to be fuzzy and varying with location of each observation in the time axis and with respect to the aforementioned conceptualization of the informational value of fSCA in the three timing windows. A fuzzy coefficient (α) was introduced in the original formulations of Pbs and the LoA schemes in order to account for the variability in informational value of fSCA on the DA results. This coefficient was described by a fuzzy degree of membership function whose value is greater than zero based on the assumption that all observations contribute certain information in constraining the perturbation parameters. An exponential trapezoidal membership function was adopted to describe α (Equation (11)). The informational value of the MODIS fSCA observations is assumed to exponentially increase in W-1 until it reaches end of this timing window at $\tau$. After start of the complete melt-out period ($c$), the informational value is assumed to exponentially decay with time in W-3. The fSCA observation is assumed to be most informative in the period between these two timing windows, i.e., in W-2. Hence, all observations in this period are assigned a maximum α value of unity.

$${\mathsf{\mu }}_{fSCA}\left(t\right)=\{\begin{array}{c}exp\left(-\frac{\tau -t}{d_1}\right),t<\tau \\ 1,\tau \le t\le c\\ exp\left(-\frac{t-c}{d_2}\right),t\ge c\end{array}$$

(11)

where ${\mathsf{\mu }}_{fSCA}\left(t\right)$ is the membership grade of each fSCA observation date index ($t$); $\tau$ and $c$ respectively refer to the change point and start of the melt-out period in the time-axis. $d_1$ and $d_2$ respectively denote the widths of W-1 and W-3 in the time-axis.

The α value of each fSCA observation is equal to the membership grade corresponding to that observation (${\mathsf{\mu }}_{fSCA}\left(t\right)$); and this fuzzy coefficient is incorporated in the likelihood measure of Pbs in a way that it rescales the level of uncertainty of the residuals corresponding to each observation (Equation (12)).

$$w_j=\frac{exp\left\{-0.5{\left[\left(y-{\widehat{\mathit{Y}}}_j\right)\alpha \right]}^T{\mathit{R}}^{-1}\left[\left(y-{\widehat{\mathit{Y}}}_j\right)\alpha \right]\right\}}{{{\displaystyle \sum }}_{j=1}^{N_e}exp\left\{-0.5{\left[\left(y-{\widehat{\mathit{Y}}}_j\right)\alpha \right]}^T{\mathit{R}}^{-1}\left[\left(y-{\widehat{\mathit{Y}}}_j\right)\alpha \right]\right\}}$$

(12)

For the LoA-based DA scheme, the membership grades of the fSCA prediction errors (${\mathsf{\mu }}_{fSCA}\left(e\right)$, Equation (8)) were rescaled through multiplication by their respective α value before taking sum of ${\mathsf{\mu }}_{fSCA}\left(e\right)$. The remaining steps are similar to the procedure followed in the LoA scheme described earlier, i.e. to get the product of ${\mathit{\text{w\_e}}}_j$ and ${\mathit{\text{w\_p}}}_j$ with subsequent rescaling of the result in such a way that the weights of all ensemble members would sum up to unity.

2.4. Detection of Critical Periods with Respect to the Informational Values of fSCA

The critical dates that form the timing windows and the fuzzy coefficient curve (Figure 2) were identified using a change point detection algorithm. The two critical dates, i.e., the point in the time axis where the change in mean snow cover of a given grid-cell occurs and start of the melt-out period in that grid cell, vary from year to year; and spatially from one grid-cell to another in response to various climatic and physiographic characteristics of an area. The use of a flexible change point detection scheme is a viable option to address the impact of such spatially and temporally variable factors on critical dates and thereby on spatiotemporal variability of SWE. The timing for start of the melt-out period in each grid cell was identified as the first incident with null fSCA observation. Detection of the CP location was realized by performing an offline change point analysis using the likelihood-based parametric approach. For the sake of comparison, change points were also estimated using the non-parametric approach. In order to minimize effect of the noisy data, the change point detection methods were applied over the transformed fSCA data (cumulative sum of fSCA). This section describes the methodology followed in detecting the change points using these two approaches.

The CP detection schemes are based on the idea that given n ordered sequence of data, $x_1:n=\left(x_1,\dots ,x_n\right)$, a change point is expected to occur within this set when there exists a time, $\tau ∊\left\{1,\dots ,n-1\right\}$, such that the statistical properties (e.g., mean or variance) of $\left\{x_1,\dots ,x_{\tau }\right\}$ and $\left\{x_{\tau +1},\dots ,x_n\right\}$ are different in some way. The parametric change point analysis scheme employed in this study was adopted from previous studies focused on implementation of the methodology based on the likelihood ratio [36,66,67]. Change point detection was presented as a hypothesis test where the null hypothesis corresponds to no change in mean; while under the alternative hypothesis, a change point is expected to occur at time $\tau$. The parametric method based on the likelihood ratio requires the calculation of the maximum log-likelihood under both null and alternative hypotheses [36]. Assuming the transformed fSCA data ($x$) as independent and normally distributed random variable, the likelihood ratio under the null ($L_{H0}$) and alternative ($L_{H1}$) hypotheses as well as the resulting log-likelihood ratio ($R_{\tau }$) can be described as follows.

$$L_{H0}=\frac{1}{\sqrt{2\pi {\sigma }^{2n}}}{\displaystyle \prod }_{i=1}^{n}exp\left[-\frac{{\left(x_i-{\mathsf{\mu }}_0\right)}^2}{2{\sigma }^2}\right]$$

(13)

$$L_{H1}=\frac{1}{\sqrt{2\pi {\sigma }^{2n}}}{\displaystyle \prod }_{i=1}^{\tau }exp\left[-\frac{{\left(x_i-{\mathsf{\mu }}_1\right)}^2}{2{\sigma }^2}\right]{\displaystyle \prod }_{j=\tau +1}^{n}exp\left[-\frac{{\left(x_j-{\mathsf{\mu }}_2\right)}^2}{2{\sigma }^2}\right]$$

(14)

$$R_{\tau }=\log \left(\frac{L_{H1}}{L_{H0}}\right)=-\frac{1}{2{\sigma }^2}[{\displaystyle \sum }_{i=1}^{\tau }{\left(x_i-{\mathsf{\mu }}_1\right)}^2+{\displaystyle \sum }_{j=\tau +1}^n{\left(x_j-{\mathsf{\mu }}_2\right)}^2-{\displaystyle \sum }_{k=1}^n{\left(x_k-{\mathsf{\mu }}_0\right)}^2]$$

(15)

where ${\mathsf{\mu }}_0$, ${\mathsf{\mu }}_1$ and ${\mathsf{\mu }}_2$ respectively denote mean value of the transformed fSCA for the whole observation period, before the change point and after the change point. ${\sigma }^2$ refers to variance of the transformed fSCA data.

The generalized log-likelihood ratio $\left(G\right)$ is estimated as:

$$G={}_{1\le t\le n}^{ max}{R}_{\tau }$$

(16)

The null hypothesis is rejected if $2G>\lambda$, where $\lambda$ is a predefined critical value; and the value of $\tau$ that maximizes $R_{\tau }$ is taken as the change point position. In this study the Bayesian Information Criterion (BIC) was employed to define this critical value, i.e., $\lambda =k·log\left(n\right)$, where $k$ is the number of extra parameters used to define the change point (here $k$ = 1).

When performing the change point analysis using a nonparametric approach, the method of Taylor [68] was followed. The procedure involves calculation of the cumulative sums ($S_i$) of the differences between the observations and their mean value (Equation (17)). The index in the time axis corresponding to the maximum cumulative sum (Equation (18)) is adopted as location of the change point $\tau$.

$$S_i=S_{i-1}+\left(x_i-\overline{x}\right)$$

(17)

$$\left|S_m\right|={}_{i=0,\dots ,n}^{ max}\left|S_i\right|$$

(18)

where $x_i$ and $\overline{x}$ respectively refer to the transformed fSCA observations and their mean value. For the first observation $S_0$= 0 and the cumulative sum ends at $S_n$= 0.

The confidence level for the change point was determined by performing bootstrap analysis. First, the magnitude of the change is estimated as the difference between the maximum and minimum $S_i$ ($S_{diff}$). Then, N = 1000 bootstraps of n observations are sampled through random reshuffling of the original data. For each bootstrap sample, the $S_i$ values corresponding to each observation as well as the difference between the maximum and minimum $S_i$ values ($S^0{}_{diff}$) are calculated. Accordingly, the percentage confidence level ($Cl$) is estimated using Equation (19).

$$Cl=100\frac{{{\displaystyle \sum }}_{i=1}^{t}score}{N}$$

(19)

where, $score=\{\begin{array}{c}1,S^0{}_{diff}<S_{diff}\\ 0,S^0{}_{diff}\ge S_{diff}\end{array}$

2.5. Evaluation Metrics

The performance of each DA scheme was evaluated through comparison of the median value of the ensemble predictions ($sim$) against the observed ($obs$) values. The following three evaluation metrics, i.e., mean absolute bias ($MAB$), root mean of squared errors ($RMSE$) as well as correlation coefficient ($R$) are employed during evaluation of estimated values against fSCA and SWE observations.

$$MAB=\frac{1}{N}{{\displaystyle \sum }}_{k=1}^N\left|ob{s}_k-si{m}_k\right|$$

(20)

$$RMSE=\sqrt{\frac{1}{N}{{\displaystyle \sum }}_{k=1}^N{\left(ob{s}_k-si{m}_k\right)}^2}$$

(21)

$$R=\frac{cov\left(obs,sim\right)}{{\sigma }_{obs}·{\sigma }_{sim}}$$

(22)

where $\sigma$, $cov$ and N respectively refer to standard deviation, covariance, and the number of observations.

3. Results

3.1. Comparison of Change Points Detected Using the Parametric and Non-Parametric Approaches

In this study, location of the change point detected using the parametric approach was also compared against those identified by the nonparametric approach. The parametric approach was based on the maximum likelihood estimate assuming that the observations are drawn from independent normal random distribution. The test for normality using chi-square on the raw and transformed fSCA data shows that although the percentage of grid cells with data that satisfy this criterion is generally low for both datasets, the result is higher for the transformed than the raw data. For example, the highest percentage of grid cells with normally distributed data at a significance level of 1% was 85% when using the transformed data as compared to 21% using the raw data. On the other hand, the non-parametric approach followed in this study is distribution free, i.e., it can be applied to any time series data regardless of the underlying probability distribution the data follows.

The change points detected using the parametric and non-parametric approaches yielded close results. Figure 3 shows the change points detected using both approaches for all grid-cells in two sample years (2012 and 2014). As mentioned in Section 2.2, these years are characterized by contrasting mean catchment SWE values during the peak accumulation period. The effect of the relatively high SWE value observed in year 2012 was reflected in the snow depletion curve of that year with the earliest grid-cells melt out day occurring after day 75 and in most of the grid-cells around day 100 (a). In year 2014, on the other hand, an early melt out date (around day 60) was observed in most of the grid-cells (b). The change point detection schemes were able to reveal this phenomenon. For example, in most of the grid-cells the $\tau$ value detected using the parametric approach was around day 25 in year 2014 as compared to day 50 in year 2012. Within a given year, the fSCA dynamics over time varied spatially from one grid cell to another. As a result, the critical dates in the ablation period, i.e., $\tau$ and start of the melt-out period for many of the grid-cells occurred at different points in the time axis. This variability in the critical dates in turn affects width of the timing windows ($d_1$ and $d_2$, Equation (11)) and thereby the relative informational value of an fSCA observation from a particular date.

Table 3. Annual statistical summary of change points ($\tau$) in observed fSCA of the grid-cells detected using the parametric (Par) and non-parametric (Npar) approaches as well as the melt-out points (c) in the time axis: Critical point statistics conducted on days since April 1.

Year	Variable					Statistics								Avg. No. Obs.
		Min		Max		Mean		Median		Variance		Skewness
		Par	Npar	Par	Npar	Par	Npar	Par	Npar	Par	Npar	Par	Npar
2008	τ	22.0	23.0	41.0	49.0	28.1	33.0	25.0	34.0	26.1	56.8	0.8	0.4	42.5
	c	57.0		94.0		65.6		64.0		67.6		2.6
2009	τ	22.0	25.0	35.0	46.0	27.9	31.2	29.0	31.0	6.1	16.0	−0.1	1.5	42.7
	c	51.0		85.0		68.4		60.0		165.6		0.3
2010	τ	33.0	36.0	40.0	41.0	36.2	37.8	36.0	38.0	2.1	2.9	0.3	0.4	35.4
	c	69.0		92.0		76.6		70.0		102.7		0.9
2011	τ	22.0	23.0	35.0	38.0	27.5	29.6	28.0	29.0	8.4	14.0	0.0	0.3	33.7
	c	50.0		96.0		68.4		69.0		116.9		0.5
2012	τ	30.0	35.0	55.0	55.0	45.2	46.0	50.0	50.0	49.9	40.6	−0.5	−0.3	32.2
	c	78.0		133.0		99.8		101.0		96.2		1.2
2013	τ	15.0	15.0	39.0	43.0	30.8	34.3	34.0	37.0	65.8	57.3	−0.8	−1.1	36.0
	c	47.0		96.0		56.5		55.0		103.0		2.0
2014	τ	24.0	25.0	45.0	47.0	28.2	37.6	26.0	40.0	22.1	54.5	1.9	−0.3	47.1
	c	56.0		94.0		60.6		58.0		70.7		2.7
2015	τ	27.0	34.0	47.0	50.0	37.2	43.0	38.0	41.0	41.5	21.6	−0.1	−0.2	35.1
	c	80.0		122.0		90.3		88.0		78.0		1.9
2016	τ	24.0	33.0	40.0	45.0	34.7	37.0	37.0	37.0	12.4	6.0	−1.2	0.5	34.6
	c	60.0		93.0		66.6		61.0		108.5		1.6
Avg.	τ	24.3	27.7	41.9	46.0	32.9	36.6	33.7	37.4	26.0	30.0	0.0	0.1	37.7
	c	60.9		100.6		72.5		69.6		101.0		1.5

and

Table 4. Annual statistical summary of change points ($\tau$) in observed fSCA of the grid-cells detected using the parametric (Par) and non-parametric (Npar) approaches as well as the melt-out points (c) in the time axis: Critical point statistics conducted on time indices.

Year	Variable					Statistics								Avg. No. Obs.
		Min		Max		Mean		Median		Variance		Skewness
		Par	Npar	Par	Npar	Par	Npar	Par	Npar	Par	Npar	Par	Npar
2008	τ	7.0	8.0	14.0	16.0	10.0	11.6	10.0	12.0	2.8	4.1	0.2	0.0	42.5
	c	11.0		24.0		16.5		17.0		10.2		0.1
2009	τ	8.0	9.0	12.0	15.0	10.5	12.4	11.0	13.0	1.6	2.3	−0.7	−0.6	42.7
	c	13.0		21.0		17.7		18.0		5.2		−0.7
2010	τ	6.0	7.0	12.0	14.0	9.0	10.1	9.0	10.0	3.2	3.9	0.0	0.1	35.4
	c	10.0		19.0		14.5		14.0		7.9		0.0
2011	τ	3.0	4.0	13.0	14.0	8.8	10.0	9.0	11.0	6.1	6.8	−0.8	−0.9	33.7
	c	7.0		21.0		15.1		16.0		14.5		−0.8
2012	τ	8.0	9.0	15.0	15.0	12.1	12.5	12.0	13.0	3.4	3.0	−0.5	−0.4	32.2
	c	13.0		25.0		19.8		20.0		9.2		−0.4
2013	τ	4.0	5.0	11.0	12.0	7.4	8.3	7.0	8.0	3.0	3.8	−0.1	−0.2	36.0
	c	7.0		17.0		11.7		12.0		6.7		0.1
2014	τ	5.0	6.0	9.0	12.0	7.5	9.1	7.0	9.0	1.2	2.0	−0.3	0.0	47.1
	c	9.0		15.0		12.1		12.0		3.0		0.0
2015	τ	5.0	7.0	12.0	13.0	8.7	9.8	8.0	10.0	2.8	3.0	0.2	0.3	35.1
	c	10.0		19.0		14.1		14.0		6.4		0.4
2016	τ	7.0	7.0	12.0	13.0	9.4	10.2	9.0	10.0	1.9	1.9	0.0	−0.2	34.6
	c	11.0		20.0		15.3		16.0		5.5		−0.1
Avg.	τ	5.9	6.9	12.2	13.8	9.3	10.4	9.1	10.7	2.9	3.4	−0.2	−0.2	37.7
	c	10.1		20.1		15.2		15.4		7.6		−0.2

shows statistical summary of $\tau$ and melt-out dates of each analysis year based on spatial domain (i.e., grid-cell values). The statistical summary was computed based on the occurrence of the critical points in days since April 1 of each year (days) and time indices that show relative position of the critical points in the time axis (index). It can be noticed that the $\tau$ detected using non-parametric approach occurred at a relatively latter point in the time axis as compared to those identified using parametric approach with an average difference of about three days (one time-step). The spatial variability is relatively higher for $\tau$ detected using the non-parametric approach with an overall average value of 30 days as compared to 26 days for $\tau$ detected using the parametric approach. Based on the bootstrap test, the non-parametric approach was able to detect $\tau$ with over 99% confidence level for all grid cells and analysis years. From this table it can also be observed that start of the complete melt-out period varied spatially from one grid-cell to another within a given year and temporally between the different years. For example, highest negative skew of −0.8 and positive skew of 0.4 were respectively observed in years 2011 and 2015 (based on the time indices). Similarly, the average number of fSCA observations per grid cell varied from year to year ranging from 32 to 47.

3.2. Evaluation of the DA Schemes Against Observed Data

We assess the effect of the different DA schemes through comparison of the median simulated against observed SWE values. For each hydrologic year, the median values resulting from each of the 100 iterations are averaged and compared against the observed values using the evaluation metrics. Figure 4 shows the evaluation result of the DA schemes in their original formulation (Pbs and LoA) as well as with fuzzy coefficients (Pbs_F and LoA_F). The DA schemes yielded an improved estimate of SWE as compared to the prior simulation in all years with the exception of years 2010 and 2015. The relative performance of Pbs and LoA as well as Pbs_F and LoA_F varied from year to year. A similar result was also observed between LoA and LoA_F. On the other hand, the incorporation of a fuzzy coefficient in Pbs resulted to an improved estimate of SWE in all years as compared to the original formulation. The fSCA-based evaluation using similar metrics also revealed that the DA schemes resulted to an improved estimate of fSCA as compared to the prior values. While an improved performance was observed when using LoA_F as compared to LoA in terms of the three evaluation metrics, a slight deterioration in MAB and RMSE was observed when using Pbs_F as compared to Pbs. In most of the analysis years, LoA has shown significant improvement both in SWE and fSCA estimation after the persistency in reproducing the observations was included as a hydrologic signature.

Figure 5 compares observed and median simulated SWE for all years under the four DA schemes and at four elevation ranges. It shows that the prior predictions tend to overestimate the low SWE values especially for grid-cells located at the highest elevation range, i.e., above 1200 masl, and underestimate the high values (a). This phenomenon was improved in the posterior SWE estimates using the DA schemes. The line fitted to the observed versus reanalyzed SWE values (red line) yielded a closer slope to the identity line (broken grey line) as compared to the prior (dark full line) (b, c, d and e). The reanalyzed values are also closer to the identity line under Pbs with a slope of 0.57 and R value of 0.59 (b) than the LoA with a slope of 0.34 and R value of 0.56 (c). However, the number of outliers in the reanalyzed values is generally lower under LoA as compared to Pbs. This is due to the fairly distributed weight between ensemble members in LoA as compared to Pbs resulting to less surprises under variable conditions. A better fit was also observed when using the DA schemes with fuzzy coefficient (Pbs_F and LoA_F) as compared to the schemes with the original formulation (Pbs and LoA). The improvement was more significant for the DA schemes based on LoA where the slope of the line fitted to the reanalyzed versus observed values has increased from 0.34 to 0.53 and the R value has improved from 0.56 to 0.65 (e). It can also be noticed that most of the simulated SWE values for grid-cells that lie in the two extreme elevation ranges, i.e., <800 m and >1200 m lie above the identity line implying that they are overestimated. This can be attributed to the errors emanating from interpolation of the point forcings in these elevation ranges. The observed low SWE values are overestimated while high values are underestimated both before and after reanalysis using all DA schemes. This phenomenon is relatively more pronounced when using LoA than Pbs.

The DA result and prior values were also analyzed by aggregating the grid cells into their respective sites (Figure 6). The average SWE value of a given site was estimated from the individual SWE values of the grid-cells in the site. Generally, a better fit between observed and simulated SWE values can be noticed at site level as compared to the analysis result at grid-cell level. For example, when using Pbs_F, the slope and R values of the line fitted to the reanalyzed versus observed SWE values have respectively increased from 0.59 to 0.80 and from 0.63 to 0.69 (i.e., an 8% increase in R). This can be attributed to the lower number of outliers at site level due to the smoothing effect of the aggregation. Similar to the analysis at grid-cell level, a better fit was observed when using Pbs as compared to LoA; and the effect of accounting for the informational value of fSCA observations was more significant in LoA_F than in Pbs_F. The value of R has increased by 16% after introducing the fuzzy coefficient in LoA. No significant spatial pattern in effect of the DA was observed with the exception that the SWE values of site 6 were less constrained by the DA schemes and the estimated SWE values of sites 1 and 4 were relatively underestimated when using the LoA approach.

Figure 7 presents a sample annual inter-comparison of the prior and reanalyzed SWE estimates against observed values. The comparison was based on site-averaged observed and simulated SWE values. The relative performance of the DA schemes varied both spatially from one site to another and temporally from year to year. For example, while the SWE values reanalyzed using Pbs were closer to the observed values at sites 7 and 9 than those from LoA in year 2008 (a), the latter was better at reproducing the observed SWE values at sites 4 and 5 for the same year. A similar example can be mentioned from a temporal perspective; Pbs showed better performance in years 2009 (b) and 2010 (c) at site 1, while LoA yielded better result in years 2011 (d) and 2012 (e) for the same site. An improvement was also observed in uncertainty of the SWE estimates after reanalysis using the DA schemes, although the level of uncertainty for a given DA scheme varied spatially from one site to another and temporally from year to year. Relatively low level of uncertainty was observed on SWE values reanalyzed using Pbs and Pbs_F as compared to LoA and LoA_F. The inclusion of fuzzy coefficient in LoA yielded to an improved level of uncertainty. On the other hand, the prior and LoA resulted in SWE values with high level of uncertainty. In Figure 7, the upper quantile (95%) values from the prior and LoA are scaled down to 10% of their respective values in order to maintain proportional look with results from the other DA schemes displayed in the same plot.

The difference in the allocation of weight between the ensemble members under the four DA schemes and the prior partly explains the difference in the width of the resulting uncertainty bound of their SWE predictions. Some ensemble members may yield very high or low SWE prediction depending on the combination of the perturbed parameter values. The effect of these predictions on the resulting quantile values, however, depends on the weight assigned to such ensemble members. In prior estimate, all ensemble members are assigned equal weight of 1/Ne. LoA also tends to maintain fair distribution of the total weight among the ensemble members as compared to Pbs, while most of the weight is carried by few ensemble members when using the latter. For example, the maximum weight carried by a single ensemble member when using Pbs was 54.6% as compared to only 1.6% in LoA-based on analysis carried out in year 2008 and a sample grid-cell (result not shown). As such, when using LoA the quantile values fairly represent the SWE predictions of several ensemble members; and this yields to wide uncertainty bound. Whereas, Pbs yields to narrow uncertainty bounds of predicted SWE, although not necessarily bracketing the observed SWE values. For example, the uncertainty bounds fail to bracket the observed SWE values in sites 4, 5, and 8 when using Pbs in year 2008 (Figure 7a). However, this phenomenon was not consistent throughout all years since this DA scheme was able to bracket the observed SWE of the same sites in other years, for example in years 2011 and 2012 (Figure 7d,e).

Performance of the different DA schemes was also analyzed in terms of their capability in reproducing the snow depletion curve (Figure 8). The fSCA values of all grid-cells are aggregated to estimate a catchment scale daily mean value. The result shows that, the complete melt-out period occurred at a later time step when using the prior estimates as compared to the reanalyzed values using the DA schemes. As such, the median prior values tend to overestimate the fSCA observations during this period. This phenomenon was observed in all years with the exception of year 2011; where the prior yielded relatively better estimate than Pbs. Neither the prior nor the posterior median estimates were able to adequately reproduce the fSCA observations in the month of April. Visual inspection of these plots also shows that the uncertainty associated with fSCA estimates was highest under LoA as compared to the other DA schemes. On the other hand, the proportion of fSCA observations bracketed by the 90% uncertainty bound was highest under this scheme. This implies that fSCA observations during anomaly years might be better reproduced by the LoA as compared to the other DA schemes. The uncertainty bound of each DA scheme also varied from year to year.

3.3. Sensitivity of the DA Evaluation Metrics to Change in Location of the Critical Points

In this study the fSCA observations that fall in different timing windows are assumed to have varying degree of informational value in constraining model states and parameters; with observations between the critical points, i.e., the detected $\tau$ and the first incidence of null fSCA, carrying the highest informational value. The validity of this assumption was examined by moving the window with highest value of fuzzy coefficient (W-2) forward and backward in the time axis, followed by subsequent evaluation of its effect on the DA result using Pbs and LoA schemes (Figure 9). The percentage change in $\tau$ was calculated using the total number of observation time indices between April 1 and the detected $\tau$ as a reference. Since, the number of observations is generally low during this period, different percentage changes in $\tau$ might point to the same nearest time index; yielding similar value in the evaluation metrics. The result shows that the evaluation metrics started to deteriorate when the location of $\tau$ was shifted forward by about 20% (i.e., −20%) and backward by about 50% in the time axis. The sensitivity was also higher when the location of $\tau$ was moved forward as compared to backward in the time axis. This can be attributed to the relatively low signal to noise ratio of the observations in W-1 as compared to those in W-3. Further, the sensitivity to the change in location of $\tau$ was higher under Pbs than the LoA. The sensitivity curves of the three DA performance metrics have shown similar trend when applied on SWE (a, c and e) and fSCA (b, d and f) in response to the displacement of $\tau$ in the time axis.

4. Discussion

In this study four ensemble-based batch smoother schemes were used to retrospectively estimate SWE during the peak accumulation period through assimilation of fSCA into a snow hydrological model. Effect of the DA schemes was analyzed for nine years and with due consideration to different elevation zones as well as two spatial scales, i.e., at grid-cell and site levels. One of these schemes was based on a Bayesian approach (Pbs), while the other was based on the limits of acceptability concept (LoA). In both approaches, model predictions were scaled proportional to weight of each ensemble member. Due to nature of the likelihood function used, most of the weights in Pbs were carried by few ensemble members, while in the LoA framework, the weights were fairly shared by all ensemble members.

Although reanalysis using these DA schemes has yielded a significant improvement of the correlation coefficient in all analysis years, only a slight improvement was observed in most of these years when accuracy of the reanalyzed SWE was compared against the prior values based on the MAB and RMSE statistical criteria. The result was generally encouraging given the high level of uncertainty associated with remote sensing data in high latitude areas due to frequent clouds and polar darkness. Other similar studies focused on assimilation of MODIS fSCA into hydrological models have also reported only a modest improvement in SWE and streamflow estimates after assimilation; and it was mainly attributed to the very short transition period between the full snow cover and the complete melt-out dates [24,20]. On the other hand, studies based on the assimilation of multiple-source remote sensing data have reported significant improvement in these evaluation metrics [28].

The relative performance of Pbs and LoA varied in response to various factors. LoA was generally more resistant to uncertainties in the assimilated data. When using low quality fSCA measurements with several missing observations, the LoA may yield better result than Pbs. For example, LoA was able to yield relatively better result as compared to the prior estimates in years 2010 and 2015 where there was an overall low performance of the DA schemes in general and Pbs in particular (Figure 4). The sensitivity test shown in Figure 9 also reveals that the Pbs performance gets easily deteriorated as compared to the LoA due to the displacement of the timing window with highest coefficient (W-2) following the changes in location of $\tau$. A further advantage of the LoA is that different intuitive evaluation criteria including hydrologic signatures can easily be accommodated in constraining model states and parameters. For example, the overall grade and persistency of an ensemble member in reproducing the observed fSCA within the assumed measurement error bounds are included as evaluation criteria in this study. However, in the presence of several correlated observations such as the fSCA observations from the complete melt-out period, LoA may suffer from identifiability problem to distinguish the best performing ensemble members and this leads to low performance as compared to Pbs. On the other hand, the likelihood with the product of squared residuals renders Pbs a better capability in identifying the optimal ensemble members under such circumstances.

The overall low performance of the DA schemes in years 2010 and 2015 can be attributed to different factors. Following year 2014, these years have lowest observed mean catchment SWE value of 419 and 458 mm respectively as compared to the other years. However, these values are not too far from the average annual value of 525 mm to explain the relatively poor performance of the DA schemes in these years. The total number of fSCA observations in these years and distribution of the observations in the time axis were also assessed since they are expected to have an effect on the DA result. Similarly, no significant difference on the number and distribution of the observations was noticed as compared to the other years. Future studies may thus assess quality of the assimilated MODIS fSCA observations and climatic inputs in relation to the other years since such factors are also expected to play an important role in the DA result.

A fuzzy coefficient was introduced into the original formulations of the two DA schemes in order to take into account for the difference in informational value of fSCA measurements from different timing windows in the ablation period. An fSCA observation in a given time index can assume different informational value depending on its location in the time axis in relation to the critical points. In Pbs and other schemes with square of errors-based likelihood function, value of the efficiency criteria is biased towards data with highest uncertainty, with the bias proportional to the uncertainty squared [30]. The fuzzy coefficient was thus introduced to rescale the level of this uncertainty with due consideration to the assumed informational value of each observation. The fSCA observations from the more informative period were promoted by allocating higher importance weight, while observations from the other periods were punished in accordance to their location within their respective timing windows. A further assumption in this study was that all observations carry some informational value; albeit to varying degree.

The informational value of fSCA observations with time was assumed to follow an exponential-trapezoidal membership function that broadly divides the ablation period into three timing windows. In W-1, the fSCA measurements close to the onset of snowmelt are generally characterized by low signal to noise ratio. Similarly, during the complete melt-out period (W-3) the measurements are homogenous; and the informational value generally decreases with correlated data that leads to less power in constraining model states and parameters. The critical points of the timing windows, i.e., $\tau$ and start of the melt-out period, vary both spatially and temporally in response to various climatic and physiographic factors. This has highlighted the need for a flexible and automatic detection of the critical points. Both the parametric and nonparametric change point detection schemes used in this study yielded similar results with well characterization of the changes. A change point analysis conducted on the transformed data (running total) was found to be more powerful and robust to outliers than the one conducted on the raw data. This can be attributed to the smoothing effect of the noises in the former case. Studies dealing with critical point detection such as the timing of snowmelt onset are also crucial to a better estimate of the surface energy balance in high latitudes, since the fraction of available solar radiation that is actually absorbed almost quadruple during the onset of snowmelt [3].

The results obtained in this study through assimilation of fSCA into a hydrologic model could also be extended to land surface models (LSMs). Inaccuracies in snow estimates within LSMs can lead to substantial errors in simulation results since snow cover has a significant impact on the surface radiation budget as well as turbulent energy fluxes to the atmosphere [69]. An assimilation of fSCA into LSMs was reported to improve the simulation result as compared to an open-loop model without assimilation [69,70]. Thus, the fuzzy logic-based DA concept presented in this study can also be applied in the assimilation of fSCA, or other variables that display fuzzy informational value with time, into LSMs in order to get an improved result. Alternatively, the use of these schemes could indirectly yield an improvement in the performances of LSMs through improvement in the hydrologic component when LSMs are coupled with the latter. In recent years, there has been growing interest in incorporating a hydrologic component into land surface models (LSMs) in order to improve weather and extreme hydrologic condition forecasts [71,72].

Generally an improved estimate of SWE was observed as a result of reanalysis using the DA schemes with a fuzzy coefficient, although, a slight deterioration in some of the evaluation metrics of fSCA was observed when using Pbs_F. This can be attributed to the relatively lower weight assigned to observations from start of the melting period coupled with the fact that most of the weights are carried only by few ensemble members. During the onset of snowmelt, sensitivity of the SDC to the initial bare ground fraction is highest and updates during this period are more likely to be dominated by this parameter than the snow coefficient of variation or precipitation [52]. The increase in fSCA error was, however, very small as compared to the improvement attained in SWE estimate. Results from this study suggest that, the concept of variable informational value of fSCA observations is consistent with the notion that many of the variables that we usually consider to be crisp quantity and deterministic are actually fuzzy that carry considerable amount of uncertainty [64].

5. Conclusions

Two ensemble-based data assimilation schemes, i.e. particle batch smoother (Pbs) and the limits of acceptability (LoA) as well as newly introduced versions of these schemes that account for the variability in informational value of the assimilated observation (fraction of snow cover area, fSCA), i.e., Pbs_F and LoA_F, were used to reanalyze the model results. Using the LoA approach as a DA scheme yielded a promising result. All DA schemes resulted in a posterior snow water equivalent (SWE) estimate that is better than the prior estimate in terms of accuracy as measured using one or more of the efficiency criteria used in this study. The result was generally encouraging given the high level of uncertainty associated with optical remote sensing data in general, and MODIS fSCA product in particular in high latitude areas.

Although all fSCA observations in the ablation period were important in constraining the perturbation parameters, some observations were more important than others depending on their location in the time axis with respect to certain critical points in the melt season, i.e., the points where the mean snow cover changes ($\tau$) and the start of a melt-out period. These critical points varied from year to year and spatially between the grid cells. The parametric (likelihood-based) and non-parametric change point detection schemes employed to locate $\tau$ in each grid cell and year yielded reasonable results. A fuzzy coefficient was introduced in the original formulations of Pbs and LoA that renders highest importance weight to the assumed most informative timing window, i.e., to observations between these critical points. A sensitivity test conducted by moving the fuzzy coefficient value assigned to this timing window forward and backward in the time axis revealed that, although somewhat insensitive up to certain distance from its original location, moving this window farther lead to deterioration of the DA results.

The reanalyzed SWE using Pbs_F was superior to Pbs throughout all years, although the degree of improvement in the evaluation metrics varied from year to year. Similarly, an improved estimate of fSCA was obtained for all years when using LoA_F as compared to LoA. On the other hand, the improvement in SWE estimate when using LoA_F in comparison to LoA as measured using a given efficiency criterion was not consistent throughout the analyses years. In many of the analysis years, the relative performance of the DA schemes was also variable from one site to another. Future studies can examine the physiographic and climatic factors affecting the relative performance of these DA schemes in order to get a better insight into the relationship between the underlying concept of the DA schemes and these factors.

Generally, incorporating the variable informational value of the remote sensing data in the DA schemes was a viable option for an improved estimates of the perturbation parameters and thereby the reanalyzed SWE values. The assumed informational values of the fSCA observations with respect to their location in the time axis was based on the general concept that observations with high uncertainty are more likely to have more informational value than certain observations. Future studies can be focused on assessment of the informational values using the well-established techniques based on the information theory and statistical measures. Although this study was focused on fSCA assimilation, the fuzzy logic-based ensemble frameworks presented in this study can be applied to assimilation of other measurements that display variable informational value with time.