Feasibility of improving a priori regional climate model estimates of Greenland ice sheet surface mass loss through assimilation of measured ice surface temperatures

The Greenland ice sheet (GrIS) has been the focus of climate studies due to its considerable impact on sea level rise. Accurate estimates of surface mass fluxes would contribute to understanding the cause of its recent changes and would help to better estimate the past, current and future contribution of the GrIS to sea level rise. Though the estimates of the GrIS surface mass fluxes have improved significantly over the last decade, there is still considerable disparity between the results from different methodologies (e.g., Rae et al., 2012; Vernon et al., 2013). The data assimilation approach can merge information from different methodologies in a consistent way to improve the GrIS surface mass fluxes. In this study, an ensemble batch smoother data assimilation approach was developed to assess the feasibility of generating a reanalysis estimate of the GrIS surface mass fluxes via integrating remotely sensed ice surface temperature measurements with a regional climate model (a priori) estimate. The performance of the proposed methodology for generating an improved posterior estimate was investigated within an observing system simulation experiment (OSSE) framework using synthetically generated ice surface temperature measurements. The results showed that assimilation of ice surface temperature time series were able to overcome uncertainties in near-surface meteorological forcing variables that drive the GrIS surface processes. Our findings show that the proposed methodology is able to generate posterior reanalysis estimates of the surface mass fluxes that are in good agreement with the synthetic true estimates. The results also showed that the proposed data assimilation framework improves the root-mean-square error of the posterior estimates of runoff, sublimation/evaporation, surface condensation, and surface mass loss fluxes by 61, 64, 76, and 62 %, respectively, over the nominal a priori climate model estimates.

runoff will be the dominant mass loss process in the future due to the retreat of the tidewater glaciers above sea level; a recent study showing that the dynamic mass loss was reduced from 58 % before 2005 to 32 % for the period between 2009 and 2012 (Enderlin et al., 2014).
Many studies (e.g., Van de Wal et al., 2012) have taken advantage of in situ measurements to provide a direct pointscale estimate of the surface mass balance (SMB, i.e., the difference between accumulation and ablation).However, with these limited in situ measurements alone, large-scale mapping of the GrIS surface mass fluxes (i.e., precipitation, evaporation, sublimation, condensation, and runoff) is impossible.The availability of remote sensing data and/or products has taken GrIS from a remote "data poor" region that is reliant mostly on sparse in situ measurements to a potentially "data rich" environment.In this regard, a key research objective is to better understand how such data can be optimally leveraged for quantitatively estimating the SMB and its associated fluxes.
Surface remote sensing data and products (i.e., surface or skin temperature, multi-frequency brightness temperature, and albedo) have been used to characterize various aspects of SMB such as snow melt, melt extent, melt duration, new snow, and extreme melt events (e.g., Abdalati and Steffen, 1995;Tedesco et al., 2011;Box et al., 2012;Hall et al., 2013).However, the relationship between surface remote sensing data/products and surface mass fluxes are most often indirect and implicit.For example, ice surface temperature (IST) can be indicative of melt, but it fails to quantitatively estimate the volume of meltwater produced.More importantly, other surface mass fluxes such as evaporation, condensation, sublimation, and runoff cannot be directly quantified via remote sensing.This makes the possibility of quantitatively characterizing the surface mass fluxes from remote sensing retrieval algorithms difficult if not impossible.It can therefore be argued that the information content of remotely sensed data remains underutilized due to indirect and implicit links between the various data streams and surface mass fluxes.
Given the limitations of the observation-based methods, numerical models offer an alternative mechanism to quantify the GrIS surface mass fluxes.Several model-based approaches have been used to characterize the spatiotemporal variability of the GrIS surface mass fluxes in both historical and future contexts (e.g., Hanna et al., 2011Hanna et al., , 2013;;Box et al., 2006;Fettweis et al., 2011;Ettema et al., 2009;Lewis and Smith, 2009;Vernon et al., 2013;Franco et al., 2013).Although the aforementioned methodologies have provided the ability to estimate the GrIS SMB and related fluxes, their estimates vary considerably, mainly due to the different physics parameterizations in the models and simplifying assumptions, the inherent uncertainty of each method, error in model and input data, and the length of data records (e.g., Rignot et al., 2011;Vernon et al., 2013;Smith et al., 2015).Therefore, it is imperative to design techniques that bridge the gap between different methods by merging rele-vant data streams with a physical model with the aim of better spatiotemporal characterization of the GrIS surface mass fluxes.In this study, we provide an example of taking advantage of information in the relevant data streams to provide a better spatiotemporal characterization of the model outputs (i.e., the GrIS surface mass fluxes).This can be done using a data assimilation approach which attempts to merge model estimates with measurements in an optimal way (Evensen, 2009).

Motivation and science questions
To date, to the best of the authors' knowledge, there have been no attempts at merging surface remote sensing data with models using a data assimilation (DA) framework to fully resolve and quantify estimates of the GrIS surface mass fluxes.Data assimilation techniques have been heavily used in hydrology to estimate soil moisture (e.g., Reichle et al., 2002;Margulis et al., 2002;Al-Yaari et al., 2014), predict snow water equivalent (e.g., Durand et al., 2008;De Lannoy et al., 2012;Girotto et al., 2014a;Zhang et al., 2014), estimate runoff (e.g., Crow and Ryu, 2009;Franz et al., 2014), improve estimates of radiative fluxes (e.g., Forman and Margulis, 2010;Xu et al., 2011), and characterize snowpack properties and freeze-thaw state of the underlying soil (Bateni et al., 2013(Bateni et al., , 2015)).DA so far has been underutilized in applications aimed at characterizing GrIS dynamics.Recently, Goldberg and Heimbach (2013) and Morlighem et al. (2013) used variational DA methods to characterize the interior and basal properties of ice sheets and ice shelves.Larour et al. (2014) assimilated surface altimetry data into the reconstructions of transient ice flow dynamics to infer basal friction and surface mass balance of the northeast Greenland ice stream.However, the use of DA for estimating GrIS SMB terms remains relatively unexplored.Assessing the feasibility of such approaches in providing a mechanism for improving quantitative estimates of SMB is the key motivation of this work.
This study utilizes an observing system simulation experiment (OSSE) framework to assess the feasibility of the proposed DA system.The OSSE framework uses synthetically generated IST measurements consistent with a "true" realization of SMB evolution.This study addresses the following science questions: (1) can assimilation of IST measurements overcome errors and uncertainties in the near-surface meteorological forcing variables for snow/ice modeling?(2) Can a DA framework be used to reduce the uncertainty and/or correct biases in a priori estimates of surface mass fluxes from a regional climate model?
This paper is arranged as follows: Sect. 3 contains the description of the models and methods used in this work.The experimental design is given in Sect. 4. The results and evaluation of the proposed methodology are discussed in Sect. 5.
Finally, key conclusions and future research directions are reported in Sect.6.

Study domain
The study domain covers the entire GrIS, which is discretized with a grid size of 25 km by 25 km to match the domain used in the regional atmospheric model described below.The focus is on fully snow/ice-covered pixels.Figure 1 shows the different GrIS mass balance zones based on a forward simulation for the year 2010.The ablation zone is defined as the region of the GrIS where the annual surface mass balance is negative.The dry snow zone is defined as the region where the mean annual temperature is less than −25 • C (Cuffey and Paterson, 2010) and melt generally does not occur.The area between the ablation zone and the dry snow zone is considered the percolation zone where surface meltwater percolates downward into the snow layers.It should be noted that the digital elevation model over the ice sheet originates from a high-resolution map generated by Bamber and Layberry (2001).The elevation of the ice sheet increases from almost 0 in the coastal regions up to about 3400 m at the summit.

Data
Surface temperature plays an important role in the coupled GrIS surface energy and surface mass budget.It is the key factor that regulates partitioning of net radiation into the subsurface snow/ice and sensible and latent heat fluxes.Surface temperature also influences the generation of runoff, the temperature profile evolution, and even basal melt (Hall et al., 2013).Space-borne instruments can provide estimates of IST.The retrieved IST is directly related to snow surface emissivity (Hook et al., 2007).The emissivity of the snow surface is a function of grain size and liquid water content, both of which are under the influence of surface processes (Hall et al., 2009).These facts support the idea that clearsky IST, of all remote sensing products available, may contain the most information about physical processes that drive the GrIS accumulation and mass loss.Therefore, this work focuses on testing the feasibility of using products such as Moderate Resolution Imaging Spectroradiometer (MODIS) IST as an extra source of information to enhance the utility of modeling techniques.The possibility of using additional remotely sensed data streams (e.g., passive microwave brightness temperature and albedo) will be investigated in future studies.
The Greenland ice surface temperature product (GrIS IST) is available from the MODIS Terra satellite (http:// modis-snow-ice.gsfc.nasa.gov/?c=greenland) and provides up to one (clear-sky) measurement per day at a native resolution of 1.5 km and an accuracy of ∼ 1-1.5

2012
).However, cloud contamination and occasional instrument outages play an important role in the availability of the MODIS IST measurements.These two factors along with some other technical and quality considerations can reduce the availability of the IST measurements to less than 10 highquality clear-sky measurements in some months (Hall et al., 2012).In the context of the OSSE used in this work, synthetic IST was generated based on the temporal resolution and acquisition time of the actual GrIS IST product by perturbing the modeled surface temperature with assumed measurement error described below.

Regional climate model (RCM)
The a priori (or prior) estimate used in the DA framework in this study is based on output from the regional climate model Modèle Atmosphérique Régional (MAR; Gallée and Schayes, 1994;Gallée and Duynkerke, 1997).The version of the model used here (i.e., MARv2) has been applied extensively over the GrIS and is described in more detail in previous studies (Lefebre et al., 2003;Fettweis et al., 2005).This version has also been used to generate future projections for the ICE2SEA European project (Fettweis et ulation (1979) and to force the atmospheric lateral boundaries as well as the oceanic conditions (surface temperature and sea ice extent) every 6 h over 1979-2010.MAR was not reinitialized every day by the ECMWF reanalysis and its results were not recalibrated after the simulation to better compare with observations as in other approaches (e.g., Box et al., 2004Box et al., , 2006)).The reader is referred to Fettweis et al. (2005Fettweis et al. ( , 2011) ) and Lefebre et al. (2003) for detailed information on the MAR setup used here.

Surface mass/energy balance and snow physical model
The key equations related to SMB are the water and energy balance of the near-surface ice sheet.The bulk surface mass balance for each model pixel (i.e., integrated over the top ∼ 10 m of the ice sheet) can be written as where P is the surface precipitation, E is the surface evaporation/sublimation, C includes both liquid and solid condensation, and R is the meltwater runoff from the snow/ice pack.Note that refreezing is implicitly included in the runoff term.Evaporation, sublimation, condensation, and runoff are the key variables that drive the surface mass loss (SML), while precipitation is the key meteorological driver for GrIS surface accumulation.
The temporal evolution of snow temperature in a vertical snow column is constrained by the conservation of energy equation, i.e., (Brun et al., 1989) where ρ is the snow density, c p is the snow heat capacity, T is the snow temperature at depth z and time t, κ is the snow heat conductivity, and q represents a sink (melt) and source (refreezing).It is worth noting that Eq. ( 2) is valid for T < 273.15 K; any energy inputs that would raise the temperature beyond freezing instead contribute directly to melt.Equation ( 2) is subject to the surface energy balance as a boundary condition, which is the key driver of the snowpack energy budget: where R ↓ s is the downward shortwave radiation, α is the (broadband) snow albedo, and R ↓ l and R ↑ l are the downward and upward longwave radiation (all terms are positive values).R n is the net radiation that is partitioned among the surface sensible (Q SH ), latent (Q LH ), and surface (Q G ) heat fluxes (into the snow).Q SH and Q LH are positive when directed toward the atmosphere and Q G is positive when directed toward the snow/ice surface.The sensible/latent heat fluxes represent the turbulent heat/vapor exchange between the surface and overlaying air due to the temperature/water vapor gradient between the surface and the reference level (i.e., meteorological forcing variables).The ground heat flux is driven by the temperature difference between the surface temperature and subsurface layers and hence can have a significant impact on the ice/snow melt and runoff.Based on Eq. ( 3), R ↓ s , R ↓ l , α and air temperature, specific humidity, and wind speed (embedded in Q SH and Q LH ) are the key meteorological variables controlling the downward energy into the snowpack (Q G ), which ultimately contributes to runoff (R).
The above coupled surface mass/energy balance represented by the CROCUS snow physical model was used in this study to provide a prior estimate of the GrIS surface mass fluxes that is consistent with the nominal forcings provided by MAR.CROCUS is a 1-D energy balance model consisting of a thermodynamic module, a water balance module taking into account the refreezing of meltwater, a turbulent module, a snow metamorphism module, a snow/ice discretization module, and an integrated surface albedo module.CROCUS derives the turbulent sensible and latent heat fluxes using a bulk method (Brun et al., 1989), which applies Monin-Obukhov similarity theory to estimate turbulent fluxes using the near-surface wind speed and the temperature and humidity differences between the surface and the temperature at ∼ 3 m, prescribed by MAR.CROCUS uses the bulk Richardson number to adapt the fluxes for stable and unstable atmospheric conditions.Note that a similar approach has been used by Van den Broeke et al., (2009a, b).CRO-CUS computes albedo and absorbed energy in each layer for three spectral bands (i.e., one visible and two near-infrared bands).The capability of the model to partition the incident solar radiation between the layers allows melt occurs on multiple depths.In CROCUS each snow layer in the snow column is treated as a reservoir with a maximum water holding capacity of 5 % of the pore volume.When the liquid water content (LWC) exceeds the threshold, excess water moves toward the layer below and the process continues until the water reaches the bottom layer and generates runoff.In addition, CROCUS takes into account changes in LWC due to snow melt, refreezing, and evaporation during a model time step.The physics of CROCUS and its validation are detailed in Brun et al. (1989Brun et al. ( , 1992)).
Assimilation of data into an RCM is another option for attempting to improve RCM fields (such as precipitation), but that is beyond the scope of this work.The focus of this work is to improve surface mass fluxes using RCM outputs and assimilation of a surface remote sensing data stream.Fur-thermore, the use of a fully coupled MAR-CROCUS system to generate an a priori ensemble estimate would be computationally prohibitive.To reduce the computational burden, an offline version of CROCUS was implemented (i.e., MAR was run over the whole modeling period, and then MAR outputs were used to force CROCUS over the same period).One can think of the DA framework outlined below as providing an update to an initial (prior) estimate of the surface mass fluxes from MAR (or any other regional climate model) using IST data as an additional constraint.
Of particular relevance to this study is the connection between CROCUS states and the measured variables used in the DA (i.e., IST).Surface temperature (synthetic IST) is an output of the forward model (CROCUS); therefore, it can directly be used as a prediction of the measurement in the DA system.One key aspect is that the raw measurements are available at higher spatial resolution than the model state (i.e., 1.5 km vs. 25 km).This was handled via an assumed change in the measurement error due to aggregation as described in more detail below.

Model adaptation
The CROCUS snow/ice model was originally developed for operational avalanche forecasting.Therefore, the model must be modified for SMB ice sheet applications.Following Fettweis (2006), the bottom boundary condition was modified for simulating approximately the top 10 m of the ice sheets.In this context, this represents the "surface" mass and energy balance via the vertically integrated states and fluxes within these top layers of the ice sheet.This method consists of the following rules.First, if during the model integration the sum of the snow and ice layer heights becomes less than 8 m, the bottom layer is extended for 2 m.Second, in the case that the sum of the snow and ice layer heights becomes larger than 15 m, the bottom layer is divided by 2. This is consistent with the methodology used in nominal MAR simulations.

Ensemble batch smoother (EnBS) framework
The EnBS is a technique that conditions a prior estimate of model states on measurements taken over an assimilation window to generate a posterior reanalysis estimate rather than a real-time (or sequential) estimate (Girotto et al., 2014a;Bateni et al., 2013Bateni et al., , 2015)).In the context of this paper, the assimilation window is a full annual cycle and measurements consist of IST data over this period.Using the generated forcing fields from MAR, the CROCUS model was run forward in time to provide an ensemble of a priori estimates of snow/ice state variables (e.g., surface temperature, snow/ice layer temperature, density, grain size) and different surface mass fluxes (e.g., evaporation, sublimation, runoff).The propagation of the CROCUS model forward in time can be shown in state-space form as where y j (t) is the vector of states for the j th realization at time t, f (.) represents the CROCUS model operator, y j (τ ) is the vector of states at previous times (τ ), u j (t) is the forcing fields for realization j , and β j is the model parameter vector for replicate j .Conventionally, the generated snow/ice states and surface mass fluxes by the forward propagation of CROCUS are called the open-loop (prior) estimates.Note that y j (τ = 0) represents the initial snow profile (IC is the initial condition).
The main source of uncertainty in a priori snow/ice states and surface mass fluxes is hypothesized to be most likely due to errors in the meteorological forcings (u j (t), see Eq. 4) generated by a parent model (in this case MAR): incoming shortwave and longwave radiation, air temperature (T a , which is implicit in the latent and sensible heat fluxes), precipitation, wind speed, relative humidity, and cloudiness.Herein, our focus is on the subset of key forcings that are the postulated main drivers of SMB (i.e., P , R l , R s , and T a ).It is hypothesized that the a priori uncertainty in forcings can be modeled via where ), and T a,MAR (x, t) are the nominal near-surface meteorological outputs from MAR; and γ P,j (x), γ s,j (x), γ l,j (x), and γ T ,j (x) are log-normally distributed multiplicative coefficients designed to capture uncertainty in the forcing inputs.The subscript j represents an individual ensemble member sampled from the postulated uncertainty distribution (j = 1, . . ., N e , where N e represents the ensemble size) and x shows the spatial index (i.e., implicitly represents an individual computational pixel in the domain).It should be noted that a multiplicative log-normal perturbation model (e.g., Margulis et al., 2002;Andreadis and Lettenmaier, 2006;Forman and Margulis, 2010) was used since all forcing variables (i.e., P , R l , R s , and T a ( • K)) are positive quantities and it provides a simple mechanism for capturing the expected uncertainty in the inputs.This type of perturbation model characterizes the ensemble using the first two moments (i.e., mean and coefficient of variation -CV) (Forman and Margulis, 2010).In this study, the mean, CV, and cross correlation between the forcing variables were obtained using the reported values in De Lannoy et al. (2010Lannoy et al. ( , 2012)).All of the parameters for each forcing are shown in Table 1.
Traditional DA applications are posed as state estimation problems where the vector of state variables (i.e., snow temperature, density, grain size, depth) is estimated via conditioning on measurements.In the current application, this can Table 1.Postulated parameters (coefficient of variation (CV) and cross-correlation) for multiplicative perturbations to hourly meteorological forcing inputs (the units for each forcing are P in mm h −1 , R s and R l in W m −2 , and T a in K).

Perturbation
CV Cross correlation 0.1 0.5 −0.3 1.0 0.6 Air temperature (T a ) 0.005 −0.1 0.3 0.6 1.0 become prohibitive since the state vector dimension is extremely large (i.e., each snow state profile involves 50 layers with several states per pixel and several thousand pixels over the domain).More importantly, updated states do not provide quantitative information about surface mass fluxes.
Hence, here we took a different approach.Rather than estimating the states directly, we treated the multiplicative coefficients γ i,j in Eq. ( 5) as the "states" to be estimated.In other words, the multiplicative coefficients have been used to transfer the nominal MAR forcing into probabilistic space (i.e., prior and posterior forcings).The DA algorithm uses IST measurements to condition the probability density function (pdf) of the prior multiplicative coefficients to compute the posterior pdf of the multiplicative coefficients.This strategy, which was also used specifically for precipitation in Durand et al. (2008) and Girotto et al. (2014a), is in direct recognition of the fact that the primary source of uncertainty in surface mass fluxes is due to error in the near-surface meteorological forcing inputs.The added benefit of this approach is that the size of the state vector is significantly reduced even in the case of time variant multiplicative states.Such a strategy derives a posterior estimate of the forcing variables directly (via the updated γ i,j ) and consequently allows for improved estimates of the surface mass fluxes via a posterior integration of CROCUS (with the posterior forcing inputs).The DA system theoretically allows the multiplicative states to vary on any arbitrary timescale.However, for simplicity, we implemented time-invariant perturbations (i.e., assumed γ i,j were unchanged over the annual modeling period) herein.In this way the update to the states was designed to allow for biases and/or low-frequency errors in individual realizations in the prior multiplicative states.
It would be ideal to characterize the uncertainties for all inputs from the information content in the assimilated data stream(s).However, in many cases available measurements are not relevant to some sources of uncertainty in the models.For instance, in this study, IST is less likely to have information about precipitation because there is no expected meaningful correlation between precipitation and IST.With regard to the fact that precipitation cannot be updated using the IST data the focus of this work has involved constrain-ing the GrIS surface mass loss (SML) components (i.e., sublimation/evaporation, condensation, and runoff), while still including the expected uncertainty in the accumulation term (precipitation).In other words, all forcing inputs were perturbed to take into account their respective postulated uncertainties, but only longwave, shortwave, and surface air temperature coefficients were updated as part of the assimilation system.
In the update step, the EnBS merges IST measurements with prior multiplicative states in order to generate a posterior estimate of those multiplicative states.In this study, we used an EnBS, which was implemented in a batch mode over a pre-defined window (i.e., applied over 1 year) with a single update.This feature of the EnBS (i.e., the batch mode update) allows running MAR and CROCUS in an offline mode that could be applied to the historical record.The open-loop (prior) estimate of the variables of interest (i.e., γ S , γ l , and γ T ) were collected into the state matrix − .Similarly, the vector of synthetically generated IST measurements was assembled into a vector: where v is the assumed additive white Gaussian error and T true is the synthetic truth (see Sect. 4.1).Finally, each ensemble member was updated individually via a Kalman-type update equation (Durand and Margulis, 2008;Bateni et al., 2013Bateni et al., , 2015)): where − j and + j represent the j th ensemble member before and after the update, respectively, and T predicted is the matrix of predicted measurements consisting of predicted IST.V is the measurement error that was synthetically produced and added to the measurements in order to avoid correlation among the replicates (Burgers et al., 1998), and K is the Kalman gain matrix which is given by where C V is the error covariance of the measurements, C T is the cross-covariance between the prior states and predicted measurements, and C T T is the covariance of the predicted measurements.In this framework, the state variables are related to the measurements in the batch through the covariance matrices that are obtained from the ensemble.The update in Eq. ( 10) can be seen as a projection of measurement-prediction misfits onto the states.The updated (posterior) multiplicative states were used in Eq. ( 5) to retrieve updated (posterior) forcing.The posterior forcings and initial snow profile (IC) were used as inputs in CROCUS to estimate the posterior surface mass fluxes.The proposed methodology can simply be extended to multiple years by applying the DA sequentially and independently for each year (e.g., Girotto et al., 2014b) or via applying the DA to a moving window (e.g., Dunne and Entekhabi, 2005).A schematic illustration of the methodology is presented in Fig. 2. The proposed methodology can be thought of as a post-processing (reanalysis) of MAR estimates by constraining the model using independent IST observations.

Experimental design
An OSSE or synthetic twin experiment offers a controlled setting in which the true forcing variables (i.e., γ S , γ l , and γ T ) are available.The goal of an OSSE is to evaluate the feasibility of the new methodology prior to assimilating real spaceborne measurements.In an OSSE, a synthetic true state and corresponding noisy measurements of the system are generated and used to evaluate the feasibility of the DA framework (e.g., Durand and Margulis, 2006;Crow and Ryu, 2009;De Lannoy et al., 2010).

True selection
The synthetic truth uses realistic input and measurement error characteristics in conjunction with the forward models to generate a realistic realization of the true system.In this study, the synthetic truth was selected as an outlier (defined below) from the generated ensemble due to the fact that errors in forcings can yield differences between a forward model (open-loop) estimate and the true surface mass fluxes.
In the OSSE system, traditionally the synthetic true ensemble is chosen from state space trajectory of the forward model (e.g., Crow and Van Loon, 2006;Durand and Margulis, 2006;Bateni et al., 2013).While an alternative approach could involve choosing the synthetic truth from the trajectory space of another well-developed RCM model, running multiple RCM models to generate a synthetic truth is prohibitive.
The ensemble of forcing data was generated via Eq.( 5) for the year 2010 and then the offline CROCUS implementation was run using the ensemble of forcing data to generate estimates of the GrIS surface mass fluxes in 2010.The year 2010 was chosen, at least in part, since it was characterized by an extreme melt rate (Tedesco et al., 2011).Considering the fact that runoff is the main component of the GrIS surface mass loss, the true ensemble (synthetic truth) was selected in a way that the integrated true runoff over the GrIS was an outlier relative to the median of the ensemble simulations.The forcing variables, states, and fluxes corresponding to the synthetic truth were also considered as the true forcings, the true states, and the true fluxes, respectively.It should be highlighted that in a synthetic DA experiment, any generated realization from the forward model (CROCUS) can be used as the synthetic truth, but one that is significantly different from the prior mean/median allows for a more robust assessment of the value of the assimilated measurements.In other words, in an OSSE the goal is to assess whether a DA framework can

Assimilated measurement characteristics
Surface temperature from the forward model can be considered as a close approximation of the remotely sensed IST.Here, the synthetic DA experiments were designed to mimic reality as much as possible.Hence, the DA system was run with a realistic representation of the temporal frequency of real space-borne IST measurements; e.g, the GrIS IST measurements from MODIS have a daily temporal resolution.However, in many instances daily observations are not available due to cloud contamination, instrument outage, and quality-related considerations.To take this issue into account, the number of available daily IST measurements (i.e., synthetic measurements) for assimilation in each month was derived from the spatial average seen in the actual Greenland IST product (e.g., Hall et al., 2012).The days with measurements were selected randomly so that the total number per month was consistent with the real number of available measurements.
Since the raw MODIS IST measurements are available at a much finer spatial resolution (i.e., ∼ 1.5 km) than the model scale (25 km), the measurements themselves and their error characteristics would require a pre-processing spatial aggregation to match the resolution of computational pixels (∼ 25 km).In the context of the OSSE in this study, the synwww.the-cryosphere.net/10/103/2016/The Cryosphere, 10, 103-120, 2016 thetic measurements and forward model both have the same spatial resolution therefore there is no need for spatial aggregation of the predicted measurement.However, specification of realistic measurement errors need to take into account the difference in spatial resolution between MODIS IST measurements and the model pixel scale.Measurement errors for MODIS IST at its raw resolution (i.e., 1.5 km) are expected to be ∼ 1-1.5 • K (e.g., Hall et al., 2012).Hence the measurement errors at the model scale (25 km) are expected to be less than or equal to this value depending on the level of correlation of the measurement errors at the sub-pixel scale.In the case of perfectly uncorrelated sub-pixel measurement errors, the aggregated measurement would be expected to have a measurement error equal to the fine-scale value divided by the number of sub-grid MODIS pixels.Assuming uncorrelated sub-grid errors are likely overly optimistic, we postulated that the measurement error standard deviation of IST at the 25 km scale is 1 K.

Implementation
The feasibility of the new DA system was evaluated via assimilation of IST as follows: a synthetically generated data stream was assimilated within an EnBS framework to assess the information content of the IST and explore whether it can overcome errors in forcing inputs.This was examined by comparing the open-loop and EnBS estimates of multiplicative states with the synthetic truth.Thereafter, the posterior meteorological forcings were fed into CROCUS to estimate the surface mass fluxes.The performance of the EnBS algorithm was further evaluated through the comparison of the posterior estimates with the prior estimates and the true estimate for all surface mass fluxes.It is worth noting that in the OSSE in this study the ensemble size was set to 100 replicates which has been shown to be adequate in previous relevant studies (e.g., Margulis et al., 2002;Huang et al., 2008;Evensen, 2009).

Performance of the EnBS via assimilation of IST
To provide an illustrative example of the methodology, Fig. 3a-c show the distribution of prior (open-loop) and posterior (obtained by assimilating IST) multiplicative state variables corresponding to the different forcings for a sample pixel (the red square in Fig. 1) in the ablation zone (latitude 67 • N, longitude 49.8 • W), which is the critical zone in terms of the GrIS surface mass loss.The prior distribution of multiplicative coefficients for each forcing variable is wide, representing the postulated uncertainty in the prior forcings.In contrast, Fig. 3a shows that the histogram of the posterior estimates of γ T is tightly distributed around the true estimate.
A narrow distribution around the true estimate means that the DA system uses the information contained in the IST se-quence and moves the ensemble members toward the true estimate while reducing the uncertainty of γ T .The reduction in uncertainty is evident by comparing the base of the posterior histogram with that from the prior estimates.The positive update by the DA system can be explained based on the fact that IST and air temperature are coupled and each one affects the other (Hall et al., 2008).Figure 3b illustrates that the median of the posterior estimate of γ l agrees well with the corresponding synthetic truth.Incoming longwave radiation is correlated with the effective (near-surface) air temperature and, as stated above, IST and surface air temperature are closely tied to each other.Prior to melt, solar radiation goes into heating the snow/ice surface; during the melt period, energy input drives sublimation or evaporation and melt (Box and Steffen, 2001).Therefore, it can be stated that IST is positively correlated with the incoming shortwave radiation.The EnBS system takes advantage of this correlation and provides improved estimates of the multiplicative state related to shortwave radiation (Fig. 3c). Figure 3d presents the time series of the IST for the prior, posterior, synthetic true, and assimilated measurements during a portion of the assimilation window.For the purpose of illustration, IST data for 10 days during the dry period (January) and beginning of the melt period (April) were selected to show the ability of the algorithm to estimate the true IST (Fig. 3d and e).It is evident in Fig. 3d-e that the EnBS captures the diurnal variability of IST and closely estimates the true IST both during the days and nights during the dry and melt periods.Moreover, Fig. 3d shows that the EnBS successfully estimates the true IST even when the temporal resolution of the IST measurements significantly decreases.This is important since the IST record shows that there are fewer measurements available during the months of December and January (Hall et al., 2012) where in some years the available measurements during these 2 months drop to fewer than 10 measurements per month.Comparing Fig. 3d with Fig. 3e also shows that during the month of January when there are fewer IST measurements the posterior estimates are in good agreement with the true IST, however, the uncertainty of the estimates is slightly larger.These results illustrate that information from IST measurements can be exploited to estimate the multiplicative states (i.e., γ s , γ l , and γ T ) and consequently the IST.
Results for the whole domain are presented in terms of relevant bulk metrics that capture the integrated impact of the forcings.Specifically, the pixel-wise cumulative incoming shortwave and incoming longwave radiation (in MJ m −2 yr −1 ) were used to represent the total energy input into the ice sheet and provide insight into the surface energy balance of the GrIS.For the air temperature, negative degreeday temperature (NDD) (i.e., cumulative mean daily air temperature for days in which the mean daily air temperature is below 0 • C) and the positive degree-day temperature (PDD) (i.e., cumulative mean daily air temperature for days in which the mean daily air temperature is above 0 • C) are two other Table 2.The spatial mean bias, the spatial RMSE, and improvement metric κ for the prior and posterior estimates of the forcing variables via assimilation of IST over the entire GrIS.metrics which are indicative of snow accumulation and melt periods, respectively.These bulk metrics were used to evaluate the performance of the DA algorithm over the entire ice sheet using root-mean-square error (RMSE) and an improvement metric.
The spatial mean bias and the spatial RMSE of the prior and posterior estimates of the integrated forcing variables over the GrIS were computed using the prior, posterior, and true cumulative longwave, shortwave, and air temperature (i.e., PDD and NDD).Table 2 summarizes the spatial mean bias and the spatial RMSE of the different forcing variables.As can be seen for the entire simulation pe-riod, the mean bias (RMSE) of cumulative shortwave, longwave, PDD, and NDD are, respectively, 84 % (70 %), 82 % (85 %), 94 % (71 %), and 65 % (86 %) less than the mean bias (RMSE) of the prior estimates.
An alternative method to evaluate the DA system is to determine the contribution of remote sensing (RS) data to the estimate explicitly.Following Durand and Margulis (2006) and Bateni et al. (2013) an improvement metric based on the prior and posterior error relative to the true was defined as follows: where the Y i (−) and Y i (+) represent the cumulative ensemble median of the prior and posterior estimates of the forcing i, respectively, and Y True i is the cumulative synthetic true for the forcing i.The improvement metric κ i can be used to interpret the contribution of the IST measurements to the posterior estimates of the forcing.This formulation suggests a value greater than 0 when the posterior error is less than the prior error (i.e., measurement improves the posterior estimates), a value equal to 0 when the prior and posterior errors are equal, and a value less than 0 when the error in the posterior estimates is greater than that in the prior estimates (the measurement degrades the posterior estimates).Table 2 shows that IST measurements make a large contribution to correct the forcing variables.IST contributed an integrated sum of 452, 375 (MJ m −2 yr −1 ), 14 and 257 ( • C day) to correct the shortwave, longwave, PDD, and NDD.The improvement metric of the PDD is much smaller than that of the NDD due the fact that there are many fewer days in which the mean daily near-surface air temperature is above the freezing point.
In order to further investigate the performance of the EnBS, the prior errors (i.e., prior − true) and the posterior errors (i.e., posterior − true) were computed for each forcing variable.Figure 4a-d show the histograms of the prior and posterior errors for cumulative R s , R l , PDD, and NDD over the spatial domain.The EnBS reduces the uncertainty of the posterior estimates for all forcing variables and effectively removes any of the prior biases.Therefore, using the improved surface energy terms to force CROCUS improves vertically integrated melt energy and enhances the estimates of the states and fluxes over the vertical snow/ice column.

Updating the SML terms
While updating the forcing variables is the mechanism by which the EnBS transfers information from IST into the posterior estimates, the main objective of the DA framework in this study is to assess the feasibility of providing better estimates of the GrIS SML and related fluxes using the improved forcings.To generate a benchmark for our analysis, CROCUS was run in open-loop mode using the prior forcings (explained above).The SML terms obtained from the prior (open-loop) simulation constitute a basis for evaluation of the methodology implemented in this study.Using the posterior forcing, CROCUS was executed for each grid cell to obtain posterior estimates of surface mass fluxes (i.e., runoff, sublimation/evaporation, and condensation) and consequently SML.
Runoff plays an important role in the GrIS net mass loss and is the main component of the GrIS SML.The GrIS meltwater runoff is heavily concentrated in the ablation zone along the ice sheet margin where the width of the ablation zone in the GrIS in some regions is very narrow and does not exceed tens of kilometers.The map of synthetic true runoff (Fig. 5a) shows that the west and southwest margins experience the highest rates of runoff that exceeds 6 m water equivalent per year.It is worth remembering that the true runoff is an outlier in the context of ensemble modeling as explained previously.Figure 5b-c show the runoff anomaly for the prior (i.e., prior − true) and the runoff anomaly for the posterior (i.e., posterior − true), respectively.The gray areas represent the percolation and dry snow zones, which do not generally contribute to surface runoff during the simulation period.It should be noted that in this area the snowmelt is not necessarily zero but refreezing can inhibit runoff.The prior anomaly map (Fig. 5b) shows that the open-loop simulation consistently underestimates the true runoff across the domain with a strong negative anomaly in the southwest margin (more than 1600 mm water equivalent below the true).Comparing the GrIS margin pixels in the prior and posterior maps (Fig. 5b-c) shows that the anomaly of the posterior estimates is significantly lower than that of the prior estimates.Reduced anomalies indicate that the EnBS successfully recovers the true estimates of the runoff in most pixels.However, the posterior results are not perfect and the algorithm slightly underestimates and overestimates runoff in some pixels.
Scatter plots of the runoff for the prior and posterior estimates versus the true estimates are illustrated in Fig. 5d-e.Each data point in Fig. 5d-e represents the ensemble median of the estimate (i.e., prior, posterior) versus the true estimate in a single pixel; the error bar illustrates the corresponding ensemble interquartile range of the estimates in the same pixel.The scatter plot of the prior runoff shows that almost all data points lie below the 1 : 1 line, indicating that the prior estimates were significantly biased (by construct in this OSSE).The posterior scatter plot (Fig. 5e) displays that the data points are narrowly distributed around the 1 : 1 line and the error bars are much smaller than that in the prior estimates, implying that the proposed algorithm significantly removes the bias and decreases the uncertainty of the estimates.
Sublimation and evaporation play an important role in the GrIS surface mass loss.However, it should be noted that MAR and CROCUS estimate surface sublimation which is considerably smaller than drifting snow sublimation.Lenaerts et al. (2012) reported for the period 1960-2011 on average surface sublimation is responsible for 40 % of total sublimation and drifting snow sublimation is responsible for another 60 %.Here, the discussion focuses on sublimation rather than evaporation due to the fact that sublimation is 1 order of magnitude larger than evaporation.The map of synthetic true sublimation (Fig. 6a) shows that the west and southwest of the GrIS in the ablation zone experience the largest sublimation rates.Box and Steffen (2001) explained that at the edge of the ice sheet, where slopes become steeper, the katabatic wind accelerates and tends to increase sublimation.Furthermore, the net radiation increases during the summertime, especially at lower latitudes, which in turn generates a vertical temperature gradient and increases the sublimation.Higher energy input also contributes to a positive albedo feedback (e.g., Tedesco et al., 2011) and further increases the sublimation rates.The prior anomaly map (Fig. 6b) illustrates that the open-loop model underestimates the sublimation at the ice sheet margin and slightly overestimates it in the ice sheet interior.The results demonstrate that posterior sublimation estimates from the assimilation of IST are much closer to the truth than are the prior estimates (Fig. 6c).Comparing the scatter plots of the posterior versus the true estimates with that of the prior versus the true estimates reveals that the methodology successfully overcomes the bias and significantly reduces the uncertainty of the sublimation estimates and increases the confidence of the results (see Fig. 6d-e).
Surface solid condensation (deposition) also influences surface mass fluxes of the GrIS by adding mass to the ice sheet.Similar to sublimation, wind and the vertical specific humidity gradient are two key factors that control the deposition.To be more precise, colder temperatures and lower winds enhance the deposition rates.In contrast with sublimation, deposition occurs at night and during winter, mainly due to radiative cooling (Box and Steffen, 2001).Figure 7a shows that the surface solid condensation (SSC) is greater in the ice sheet interior where winds are weak and there is sufficient moisture in the air column.The high elevation central regions, however, show less condensation due to distance from moisture sources.High speed winds in the ice sheet margins prevent condensation despite the availability of moisture.Figure 7b shows that the prior estimates for SSC is not in good agreement with the truth and that the prior simulation both underestimates and overestimates surface solid conden-  the IST measurements to eliminate the bias and reduce the uncertainties of the posterior estimates.
Herein, the SML is defined as the sum of the mass loss terms (i.e., runoff and sublimation/evaporation) and mass gain term (i.e., surface solid condensation) discussed above.Figure 8a shows that SML is greater in the west and southwest of the ice sheet where runoff is the dominant mass loss mechanism and is smaller in the ice sheet interior where mass loss mainly occurs through sublimation.Similar to runoff, the prior anomaly is largely concentrated in the ablation zone and, since runoff is roughly 2 orders of magnitude larger than sublimation and condensation, the anomaly due to these two fluxes is almost undetectable in the anomaly map (see Fig. 8b).Comparing the posterior anomaly map (Fig. 8c) with that of the prior clearly shows that the posterior SML is closely matched with the true estimates across the domain.Scatter plots (Fig. 8d-e) also confirm that the EnBS effectively removes the bias and increases the confidence level of SML estimates.
To provide an integrated picture over the full domain, Fig. 9a-d show the time series of the cumulative runoff, sublimation, surface solid condensation, and SML over the GrIS, respectively, in 2010.As illustrated in Fig. 9a, the true runoff starts in late April and increases rapidly during the melt season (to a cumulative value of 408 mm) until late August.The central tendency of the prior simulation (as indicated by the ensemble median) underestimates the runoff by about 35 % owing to errors in the forcing inputs.The posterior estimates show a cumulative runoff of 394 mm over the GrIS, which is in good agreement with the truth.Table 3 shows that the EnBS reduces the spatial mean bias (RMSE) of the prior estimates of runoff by 90 % (61 %) from −552 mm (646 mm) to −54 mm (250 mm).Note that runoff occurs in the ablation zone therefore the spatial mean bias and spatial RMSE for runoff were computed over the ablation zone.The spatial mean bias and spatial RMSE for sublimation, condensation, and SML were computed over the entire ice sheet.As evident in Fig. 9b, sublimation accelerates during the summer season owing to increased energy input to the snow/ice surface.The true estimate suggests that in total net sublimation (i.e., sublimation and evaporation) accounts for about 66 mm (∼ 15 %) mass loss over the GrIS.The median of the prior simulation shows a total sublimation loss of ∼ 56 mm which is 10 mm less than the truth.The EnBS significantly im- proves the results where the posterior median estimate shows a total sublimation of 65 mm.From Table 3 the spatial mean bias (RMSE) of the posterior estimate shows a 90 % (64 %) reduction relative to the prior.In general surface solid condensation accelerates during the winter and decelerates in the summer season (Fig. 9c).The true simulation suggests a cumulative SSC of 27 mm, and the median of the prior and posterior estimates is 25 and 27 mm, respectively.The 76 % reduction of the spatial RMSE of the posterior estimates and 80 % reduction of the spatial mean bias (Table 3) also supports the accuracy of the posterior estimates.Finally, the true SML estimate is 450 mm, and the prior and posterior median of SML are 295 and 435 mm, respectively.Clearly the posterior SML estimate is in better agreement with the truth.The IST measurements contribute an integrated sum of 140 mm to correct the posterior estimates of the GrIS SML and also reduce the spatial mean bias and the spatial RMSE of the estimates by 90 and 62 %, respectively (Table 3).
A probabilistic approach also provides information about the uncertainty of the estimates.Figure 9a-d

Sensitivity to the synthetic truth values
As in any OSSE, the synthetic measurements are, by construct, a function of the chosen true and therefore the posterior results could be impacted by the particular selection of the true realization.To address this concern, and show the robustness of the proposed algorithm, the simulation was repeated for two different true values: one smaller than the baseline simulation and the other larger.In the first case the synthetic true runoff was set to 330 mm, which is the average of the runoff estimates from the open-loop simulation (i.e., ∼ 260 mm) and the true runoff from the baseline simulation (i.e., ∼ 400 mm).In the second case the true runoff was set to 470 mm, which is 70 mm larger than the baseline simulation.Table 4 shows the RMSE of the surface mass fluxes for all simulation cases.The posterior RMSE of each mass flux for all simulation cases are very similar even when the prior RMSE of the estimates are significantly different.For example, the prior RMSE of the runoff (SML) for the second simulation case (true runoff equal to 470 mm) is 2.5 (2.6) times larger than the prior RMSE of the first simulation case (true runoff equal to 330 mm), but the posterior RMSE differs by only 4 % (10 %).Therefore, it can be stated that the DA algorithm robustly retrieves the true estimates of the surface mass fluxes and the performance of the algorithm is relatively insensitive to the selected truth.

Discussion and conclusions
A new data assimilation methodology for improving estimates of the GrIS surface mass loss fluxes has been tested and presented using an observing system simulation experiment framework.The prior estimates were derived from an offline surface module (CROCUS) forced by an ensemble of meteorological forcing fields that were based on a nominal regional climate model simulation (in this case MAR).A posterior estimate was generated by conditioning the forcings on the synthetically generated IST measurements using an ensemble batch smoother (EnBS) approach.Specifically, it was shown that using the EnBS with IST measurements, we were able to improve nominal estimates derived from MAR that result from erroneous forcing fields that drive surface mass and energy balance processes.The results illustrated that IST measurements have potential information on shortwave, longwave, and surface air temperature that allows for correction of errors in these terms.However, due to the lack of meaningful correlation between precipitation and IST measurements, the precipitation flux was not updated in this context (i.e., the prior and posterior precipitation is the same).Hence the assimilation of IST is primarily beneficial for estimating the surface mass loss terms and not the accumulation term.However, it should be noted that using MAR-CROCUS to generate the synthetic truth might lead to optimistic results since the truth is taken from the same model.Mitigation of this was attempted by using an outlier for the truth.An expensive alternative, but worth pursuing in future work, would be to use other RCM models to generate the synthetic truth.That said, it can be argued that using another model such as RACMO2 (Ettema et al., 2009) to generate the true realization will not significantly affect the results because the synthetic truth from RACMO2 is likely to fall within the ensemble spread of MAR-CROCUS trajectory.The main reasons for that are (1) the SMB fluxes from MAR and RACMO2 are highly correlated (Fettweis et al., 2013) and (2) the trends of SMB fluxes from two models are very similar (Vernon et al. 2013).Furthermore, sensitivity analysis shows that the proposed algorithm is able to retrieve the synthetic truth for the extreme cases where the real true stats fall beyond the chosen values.
The new methodology has several advantages over the traditional state-space data assimilation approaches.First, in this new application the multiplicative perturbation variables are considered as states to be updated.Reduction of the size of the state vector and consequently computational costs is the direct outcome of this approach.Second, mass loss terms cannot directly be sensed by the means of satellite sensors; using this methodology, the mass loss fluxes were estimated indirectly by reducing the error in forcing variables.Finally, the modularity of the proposed methodology would allow for incorporation of any regional climate model and additional remotely sensed observations in future applications.All of these advantages should make such data assimilation approaches attractive and complementary approaches to better resolve and diagnose the ice sheet surface mass fluxes.The improved mass loss estimates could also be used as input to net mass balance estimates and ultimately a sea level rise projection when applied to real data over the remote sensing record.
As a final note, it should be emphasized that the application presented in this study does not attempt to optimize or include uncertainty in any model parameters.Rather, the focus is on the uncertainty of time-varying model forcing inputs, which is expected to be the primary source of uncertainty in estimates of surface melt.We acknowledge that the model parameters are treated as certain and, therefore, any uncertainty/error in model parameters (e.g., water holding capacity that impacts the transformation of meltwater into runoff) would increase the expected error in posterior SML in an application with real data.A more general case where estimation of parameters is included in the data assimilation framework could be the basis of future work.
The next logical step is to apply the methodology with real IST measurements to further validate the robustness of the proposed approach.This future work will include the use of the MODIS IST product for estimating GrIS SML.The data assimilation framework is general and could also include the potential application of assimilation of passive microwave, albedo, and even Gravity Recovery and Climate Experiment (GRACE) data to further constrain GrIS SMB estimates.

Figure 1 .
Figure 1.The Greenland ice sheet mask (filled area), including the ablation zone (blue), the percolation zone (dark green), and the dry snow zone (bright green), based on an offline CROCUS simulation for the year 2010.The contour lines show the topography of the ice sheet with an interval of 500 m.The red square shows the location of a sample pixel in the ablation zone (see Sect. 5.1).

Figure 2 .
Figure 2. Schematic illustration of the proposed methodology.The posterior SMB/SML is effectively a post-processing (reanalysis) of regional climate model (in this case MAR) estimates conditioned on IST measurements.The term IC represents the initial snow profile.

Figure 3 .
Figure 3. Ensemble histogram of the prior (red bars) and the posterior (after assimilation of IST) multiplicative states (blue bars) for (a) surface air temperature, (b) longwave radiation, and (c) shortwave radiation for a sample pixel in the ablation zone.The prior (red line) and posterior (blue line) median values and truth (black line) are also shown for reference.Also shown are the time series of (d) the IST for the 10-day period during the dry season and (e) the IST for the 10-day period during the melt season.The red and blue shaded areas represent the prior and posterior uncertainty band (interquartile range) and the red, blue, and black lines represent the median of the prior, the median of the posterior, and the truth, respectively.The green circles represent the synthetically generated (noisy) IST measurements that are assimilated to generate the posterior estimates.

Figure 4 .
Figure 4.The histogram of the prior errors (red) and posterior (after assimilation of IST) errors (blue) for cumulative (a) shortwave radiation, (b) longwave radiation, (c) PDD, and (d) NDD over the full GrIS.

Figure 5 .
Figure 5.The (a) synthetic true runoff (mmWE yr −1 ) for the year 2010, (b) runoff anomaly (mmWE yr −1 ) for the prior (i.e., difference between the prior and true runoff), (c) runoff anomaly (mmWE yr −1 ) for the posterior, (d) scatter plot of the prior runoff estimates, (e) scatter plot of the posterior runoff estimates.Black dots are the ensemble median of the estimates and the error bars represent the corresponding ensemble interquartile range of the estimates.

Figure 6 .
Figure 6.The same as Fig. 5 but for sublimation and evaporation.

Figure 9 .
Figure 9.The time series of (a) cumulative runoff, (b) cumulative sublimation and evaporation, (c) cumulative surface solid condensation, and (d) cumulative mass loss over the GrIS (in millimeter of water equivalent).The truth is the black dashed line, the prior ensemble median is the red line, and the posterior ensemble median is the blue line.The red shaded area corresponds to the ensemble interquartile range (IQR) for the prior simulation and the blue shaded area corresponds to the ensemble IQR for the posterior estimates.
show that the prior estimates of all surface mass fluxes have a large ensemble spread, reflecting the propagation of a priori forcing uncertainties to SML terms.During the update process the EnBS significantly reduces the uncertainties of the posterior estimates of forcing variables and consequently the posterior estimates of the surface mass fluxes.Comparing the narrow blue shaded area with the wide red shaded area illustrates that the EnBS increases the confidence of the model predictions by decreasing the error and uncertainties of the posterior estimates relative to the prior estimates.

Table 3 .
The spatial mean bias and the spatial RMSE of runoff, sublimation/evaporation, surface solid condensation, and net mass loss estimates via assimilation of IST measurements.The spatial mean bias and the spatial RMSE for runoff were computed over the ablation zone and for the other surface mass fluxes were computed over the entire ice sheet.

Table 4 .
The spatial RMSE of runoff, sublimation/evaporation, surface solid condensation, and net mass loss estimates via assimilation of IST measurements for three different true values.