1. Introduction
Data assimilation has been widely used in meteorology and oceanography. It has proven to be an effective method for improving knowledge of the past, present and future state of the atmosphere and ocean. Assimilation of observations from very different sources improves operational weather forecasts and allows reanalysis of the state of the atmosphere and oceans over the past century (e.g. Reference CourtierCourtier and others, 1998; Reference UppalaUppala and others, 2005).
Recently, data assimilation has begun to be used in glaciology. For instance, Reference Arthern and GudmundssonArthern and Gudmundsson (2010) show that the basal friction parameter of ice sheets can be inferred using an inverse method, based on measurements of surface velocity. Several authors have also studied the assimilation of remotely sensed data in the scope of hydrological modeling (e.g. Reference Andreadis and LettenmaierAndreadis and Lettenmaier, 2006; Reference Slater and ClarkSlater and Clark, 2006; Reference Durand, Molotch and MargulisDurand and others, 2008; Reference De Lannoy, Reichle, Houser, Arsenault, Ver- hoest and PauwelsDe Lannoy and others, 2010; Reference ToureToure and others, 2011). These studies have demonstrated the positive impact of assimilating remotely sensed data on snowpack modeling and hydrological forecasting.
Data assimilation is a technique that allows a priori information on the state of a system to be combined with newer information, observations in our case (Fig. 1). The a priori information is often, and is usually referred to as, a guessed field. An a posteriori state is then developed from mixing the guess and the observed fields, i.e. the analysis. In this study, the guess field is known from model simulations. The mixing of the two sources of information is done according to the error statistics (model and observations) and to a direct relationship, which can be nonlinear, between the observed variables and the vector that describes the model state.
Fig. 1. Simplified schematic view of the methodology used to reconstruct the spatial distribution of the glacier mass balance.
Accurate knowledge of the spatially distributed mass balance of temperate glaciers is crucial for a better understanding of the physical processes controlling the mass balance, water-resource monitoring and reanalysis of past surface energy balances (SEBs). This knowledge is presently limited by a lack of observations. The temperate glacier SEB is mainly explained by the solar radiation budget during the ablation season (Reference Sicart, Hock and SixSicart and others, 2008). In addition, variations of the solar radiation budget (shortwave radiation balance) are governed by variations of surface albedo (Reference Sicart, Hock and SixSicart and others, 2008). Note that glacier surface albedo has a high spatial and temporal variability, with values extending from 0.15 for dirty ice to 0.9 for new snow (Reference Oerlemans and KnapOerlemans and Knap, 1998). Albedo can be observed using satellite sensors (Reference Dozier, Green, Nolin and PainterDozier and others, 2009). Consequently, albedo is a good candidate for assimilation in a snow model to improve simulation of glacier mass balances.
Over recent years, several studies have focused on the simulation of spatially distributed glacier mass balances (e.g. Reference HockHock, 2005; Reference AndersenAndersen and others, 2010; Reference Giesen and OerlemansGiesen and Oerlemans, 2010). Mass-balance models vary from simple temperature-index models to so-called energy- balance models (Reference HockHock, 2005), that evaluate the surface energy fluxes in detail. The latter approach is used in this study, since it allows detailed evaluation of each flux involved in the SEB and direct estimation of mass balance.
Modeling snowpack variations with time is generally a difficult step. Existing snow models have been reviewed and their performance compared by Reference Essery, Martin, Douville, Fernandez and BrunEssery and others (1999), Reference Essery and EtcheversEssery and Etchevers (2004) and Reference EtcheversEtchevers and others (2004). All these authors have pointed out the major influence of the albedo parameterization on the final mass-balance estimates. Reference Klok and OerlemansKlok and Oerlemans (2002) and Reference HockHock and Holmgren (2005) used models to simulate spatially distributed mass balances with two different albedo parameterizations and, once again, emphasized the crucial role of this variable.
Reference Gerbaux, Genthon, Etchevers, Vincent and DedieuGerbaux and others (2005), Reference LejeuneLejeune and others (2007) and Reference LejeuneLejeune (2009) developed a spatially distributed version of the CROCUS snow model (Reference Brun, Martin, Simon, Gendre and ColeouBrun and others, 1989, Reference Brun, David, Sudul and Brunot1992). Reference Gerbaux, Genthon, Etchevers, Vincent and DedieuGerbaux and others (2005) tested the model on Glacier de Saint-Sorlin, French Alps, the test site in this study. As input, they used meteorological variables from the SAFRAN meteorological analysis system (Reference Durand, Brun, Merindol, Guyomarc’h, Lesaffre and MartinDurand and others, 1993). SAFRAN disaggregates large-scale meteorological analysis and observations in the French Alps. Reference LejeuneLejeune and others (2007) and Reference LejeuneLejeune (2009) tested the model on the Glaciar Zongo (Bolivia) basin, carefully parameterizing the albedo in the case of shallow snow cover and using input meteorological data from automatic weather stations (AWSs).
All these studies point out the importance of the albedo scheme, the difficulties encountered in its parameterization and the complexity of comprehensive albedo field measurements, due to high spatio-temporal variability (Reference Pedersen and WintherPedersen and Winther, 2005; Reference Gardner and SharpGardner and Sharp, 2010). Some processes, such as those determining the snowline position, the albedo reduction at the end of the ablation season due to dust and algae accumulation, and the elevation of the rain/snow limit, are especially difficult to model but have a strong impact on the final mass-balance estimates. For this reason, we attempt to assimilate, in the SEB model, observed glacier albedo from remote-sensing data to better account for all these processes and thereby improve the mass-balance estimates.
1.1. Intent of the study
In this study, the spatially distributed version of CROCUS adapted to glaciers by Reference Gerbaux, Genthon, Etchevers, Vincent and DedieuGerbaux and others (2005) is tested for assimilation of the albedo observations from spatial (Moderate Resolution Imaging Spectroradiometer, MODIS) and terrestrial (digital photographs) remote-sensing measurements using a variational data assimilation scheme (Reference Bouttier and CourtierBouttier and Courtier, 2002).
The following section describes the study site, the albedo maps and the meteorological data. The methodology is then explained in detail. Section 4 introduces the results and provides details about the variational method at one point of the Col de Porte study site. The impact of the assimilation on the simulation of Glacier de Saint-Sorlin spatially distributed mass balance over five hydrological years selected within the 2000-10 decade is then presented and the results discussed.
2. Study Site and Data
2.1. Col de Porte
Col de Porte is a test site of Meteo-France located at 1326ma.s.l. in the Chartreuse area, near Grenoble, France (Reference Brun, Martin, Simon, Gendre and ColeouBrun and others, 1989, Reference Brun, David, Sudul and Brunot1992). This site was chosen for its comprehensive set of meteorological and snow-profile (weekly) observations, which are used for preliminary tests and validation of the method developed in this study.
2.2. Glacier de Saint-Sorlin
Glacier de Saint-Sorlin (Fig. 2) is located in the French Alps (Grandes Rousses area) at 45.10° N, 6.10° E. It is a small temperate glacier covering an area of ˜3km2. The tongue is at 2700 m a.s.l. and the top reaches the Pic de l’Etendard (3463 m a.s.l.). The glacier mass balance has been measured since 1956 (Reference Vincent, Vallon, Reynaud and Le MeurVincent and others, 2000; Reference VincentVincent, 2002). Numerous dynamic modeling and climate studies have been conducted on this glacier (e.g. Reference Le Meur and VincentLe Meur and Vincent, 2003; Reference Le Meur, Gerbaux, Schafer and VincentLe Meur and others, 2007). Detailed information on this site can be found at http://www-lgge.ujf-grenoble.fr/ServiceObs.
Fig. 2. Digital elevation model of Glacier de Saint-Sorlin. The outline of the glacier in 2003 is indicated by black stars. The locations of temporary and permanent AWSs are indicated by white crosses. Black dots, crosses, circles and diamonds represent stakes for mass-balance measurements for the four sections mentioned at the top. The SEB measurement station during summer 2006, SEB2006, was located near AWSabla, 2008
Two types of measurements available on this glacier are used in this study:
Automatic weather stations
A permanent AWS (AWSmoraine) has operated since 2005 on the moraine near the hut (Fig. 2). It provides half-hourly measurements of air temperature, relative humidity, wind speed and direction, along with downward and upward, longwave and shortwave radiations. Additionally, a complete SEB measurement station (SEB2006) was set up on the glacier during summer 2006 (Reference Sicart, Hock and SixSicart and others, 2008; Reference Six, Wagnon, Sicart and VincentSix and others, 2009). The SEB2006 station measured the turbulent fluxes of sensible and latent heat using an eddy covariance method. A Campbell CSAT-3 sonic anemometer and a LI-COR LI- 7500 infrared gas analyzer were installed 2 m above the glacier surface. A Campbell CR1000 data logger recorded the wind, temperature and water-vapor data at 20 Hz frequency. The eddy covariance data were processed with the EdiRe software. During summer 2008 and 2009, temporary radiation measurement stations (AWSabla, 2008 and AWSaccu, 2008 during summer 2008 and AWSabla200g during summer 2009) measured surface albedo using Kipp & Zonen CNR1 and CMP3 devices in the accumulation and ablation zones of the glacier. The locations of these AWSs are shown in Figure 2.
Mass-balance measurements
Stakes are used to measure the mass balance at different points of the glacier using the direct glaciological method (Reference Cogley and AdamsCogley and Adams, 1998). Locations of the 33 stakes for summer 2008 are indicated as an example in Figure 2. The global estimated accuracy of the annual mass- balance measurement with stakes is ±20 cm w.e. (Reference Cogley and AdamsCogley and Adams, 1998). For discussion purposes, the stakes are divided into four networks: the ablation area (diamonds), the intermediate zone (dots), the dense network (circles) characterized by rugged topography, and the accumulation zone (crosses) (Fig. 2). During the ablation season, intermediate mass-balance measurements are made approximately every 15 days.
2.3. Albedo maps derived from remote-sensing data
The albedo maps used in this study are derived from two data sources: terrestrial photographs (Reference Dumont, Sirguey, Arnaud and SixDumont and others, 2011) and MODIS data, based on the work of Reference Sirguey, Mathieu and ArnaudSirguey and others (2009). Two terrestrial digital cameras have automatically taken photographs of the glacier from the hut (Fig. 2) since summer 2008, allowing estimation of surface albedo at a spatial resolution of 10 m. The root-mean-square error (rmse) from comparison with field measurements made during two summers is <0.07 on the broadband albedo value (Reference Dumont, Sirguey, Arnaud and SixDumont and others, 2011). The method developed for the photographs is also applied to MODIS data at a spatial resolution of 250 m and the estimated broadband albedo rmse is 0.06 compared with the field measurements for the same period.
Table 1 provides an overview of the albedo maps used in the assimilation scheme for the five hydrological years selected in this study from the 2000-10 decade. The studied hydrological years were selected for the following reasons. The years 2000/01 and 2002/03 were selected because they represent extreme events for the glacier during the decade (positive annual mass balance and strongly negative mass balance, respectively). The year 2005/06 was chosen because the complete SEB measurement station was set up on the glacier during summer 2006 and the two last years were selected because terrestrial photographs of the glacier were available for this period.
Table 1. Overview of the albedo maps available for assimilation for each hydrological year. Data are chosen during the ablation period (May to October). The albedo maps are derived from MODIS data (250 m spatial resolution) or from terrestrial photographs (10 m spatial resolution) (Reference Dumont, Sirguey, Arnaud and SixDumont and others, 2011)
2.4. Meteorological data (SAFRAN)
The objective of this study is to develop a method that can be used on a large set of glaciers. The selected input data are therefore taken from disaggregation of a large-scale analysis and observations, as explained below.
SAFRAN is an original meteorological analysis system, designed to provide input meteorological data to snow models used in mountainous areas. It estimates the hourly values of temperature, wind speed, relative humidity, precipitation amount and phase, direct and diffuse solar radiation, nebulosity (global cloudiness) and longwave downward flux (Reference Durand, Brun, Merindol, Guyomarc’h, Lesaffre and MartinDurand and others, 1993). The meteorological variables are given at 300 m elevation steps and for the seven possible orientations north, east, southeast, northeast, southwest, west and flat (i.e. no orientation). They are representative of different mountainous regions in France, with areas varying from 400 to 1000 km2. For this study, data for the Massif des Grandes Rousses were selected and linearly interpolated on a 100 m resolution grid of the digital elevation map of Glacier de Saint-Sorlin, according to the method described by Reference Gerbaux, Genthon, Etchevers, Vincent and DedieuGerbaux and others (2005). Multiple reflections from adjacent mountains were not taken into account for either shortwave or longwave radiation. Shading was taken into account, but only for shortwave radiation.
SAFRAN estimates have been previously compared with measurements at AWSmoraine for three hydrological years (2005/06, 2007/08 and 2008/09). This comparison shows that temperature, relative humidity and shortwave incident radiation estimates are in good agreement with AWS measurements. The determination coefficients are, respectively, 0.98, 0.77 and 0.79 for hourly values. The comparison is less convincing for wind speed and longwave downward radiation. The AWS is located on the top of a hill, which implies a difference between the measured and estimated wind speed. This difference is due to local orographic features which are not taken into account in SAFRAN estimates. The longwave downward radiation is strongly influenced by local nebulosity. We found that SAFRAN overestimates the longwave radiation in cases of low nebulosity (positive bias of 17 W m-2 for the three years) and we therefore implemented a simple correction of SAFRAN longwave downward radiation as a function of nebulosity.
The use of SAFRAN precipitation fields leads to underestimation of the winter snow accumulation (difference of scale and topography effects). Precipitation amounts are consequently adjusted. This adjustment is based on the work of Reference Gerbaux, Genthon, Etchevers, Vincent and DedieuGerbaux and others (2005), i.e. using measurements of winter accumulation from stake readings. Fourteen years of winter mass-balance measurements (1995-2009) were used to compute the multiplication factors that allow reconstruction of the spatial distribution of winter accumulation over all these years. The multiplication factors are interpolated linearly with respect to latitude and longitude over the whole glacier. The impact of the interpolation method was tested and the mean difference between linear and cubic interpolation was found to be <1 mm w.e. for 1 year over the whole glacier (maximum value 8cm w.e.). We therefore used the linear interpolation method.
For the five hydrological years, the corrected precipitation field provides an unbiased estimation of winter accumulation. The rmse is 0.20 m w.e. with respect to 132 stake measurements. The mean multiplication factor is found to be 1.64 over the whole glacier. For comparison, Reference Gerbaux, Genthon, Etchevers, Vincent and DedieuGerbaux and others (2005) found a factor of 1.5 for the whole Reference Gottardiglacier and Gottardi (2009) found a value of 1.6 for solid precipitation in this massif. The precipitation correction factor takes into account the difference of scale (glacier vs massif) and the high-altitude effects on the precipitation field estimates.
3. Methodology
This section describes the methodology used to assimilate albedo data in the simulated snowpack. The method developed is based on the spatially distributed version of CROCUS adapted to Glacier de Saint-Sorlin (Reference Gerbaux, Genthon, Etchevers, Vincent and DedieuGerbaux and others, 2005). Figure 1 gives a schematic overview of the method.
3.1. CROCUS, a brief overview
The CROCUS snow model was developed by Reference Brun, Martin, Simon, Gendre and ColeouBrun and others (1989, Reference Brun, David, Sudul and Brunot1992) for operational avalanche forecasting. CROCUS is a one-dimensional snowpack model driven by meteorological variables alone and models each layer of the snowpack. Each layer is characterized by its temperature, density, liquid water content, depth, a historical variable for grain metamorphism and grain variables including size, sphericity and dendricity (Reference Brun, David, Sudul and BrunotBrun and others, 1992). Dendricity varies from 1 to 0 and describes the fraction of original precipitated crystal forms still present in the snowpack. Sphericity varies from 0 for completely faceted grains to 1 for rounded particles. Given the objective of the present study, we detail only the representation of turbulent fluxes and albedo in the model.
3.1.1. Turbulent flux representation
CROCUS uses a bulk formulation for turbulent fluxes (Reference Martin and LejeuneMartin and Lejeune, 1998). All roughness lengths, i.e. moisture and temperature in the CROCUS model, are assumed to be equal to an effective roughness length, z 0 (Reference Martin and LejeuneMartin and Lejeune, 1998). This assumption is not invalidated by Reference AndreasAndreas (1987), who shows that the moisture and temperature roughness lengths are only slightly different over all Reynolds numbers. In this study, z 0 is 3 mm for snow and 6 mm for ice, which is consistent with measurements over glaciers given by Reference MartinMartin (1975), Reference Greuell and SmeetsGreuell and Smeets (2001) and Reference Brock, Willis and SharpBrock and others (2006). These roughness lengths are not considered as accurate physical values, but as effective ones, i.e. they are used to slightly tune the mass-balance model while keeping a realistic value.
3.1.2. Albedo scheme
The snow albedo is parameterized as a function of the surface layer age, which is an implicit way to take into account the effect of impurities and of the optical grain size (Reference Brun, David, Sudul and BrunotBrun and others, 1992). The optical grain size is computed from grain variables (Reference WillemetWillemet, 2008). The albedo is split into three spectral bands, αi : [0.3-0.8], [0.8-1.5] and [1.5-2.8] |µm, for energy-balance purposes. In the operational version, the visible albedo, α1 , over the first spectral band stands between 0.7 and 1. The visible albedo decreases as a function of the optical diameter and the age of the layer up to a limit of 90 days, and is bounded to a minimum value of 0.7 for snow (Reference Brun, David, Sudul and BrunotBrun and others, 1992). In this study, the minimum albedo value in the first spectral band is set to 0.6, to be able to represent firn and very old snow with a high light-absorbing impurity content, in the late summer in the upper zone of the glacier. Reference PatersonPaterson (1994) cited a similar mean broadband firn albedo of 0.53.
The albedo of ice is simulated as a constant. In the operational version of CROCUS, α1,ice = 0.45, α2,ice = 0.30 and α3,ice = 0.1 (Reference Brun, David, Sudul and BrunotBrun and others, 1992). On Glacier de Saint-Sorlin, Reference Gerbaux, Genthon, Etchevers, Vincent and DedieuGerbaux and others (2005) used [0.23, 0.16, 0.05], based on measurements performed on the glacier. The latter values were used in this study because they give, without albedo assimilation, the best agreement between the simulated and measured mass balance.
In addition, the parameterization to take into account ice albedo for shallow snowpacks is adapted from the work of Reference Wagnon, Lafaysse, Lejeune, Maisincho, Rojas and ChazarinWagnon and others (2009).
3.2. Variational assimilation of snow albedo data
The observed albedo is assimilated in CROCUS using a variational data assimilation scheme (Reference Bouttier and CourtierBouttier and Courtier, 2002). This scheme modifies the simulated surface layer according to both the albedo observations and the first simulated (‘guess’) surface layer grains that are not linearly linked. It takes into account modeling and observation errors. This means that the significance given in the analysis to the observations and to the guess state is a function of the observation errror statistics and of the model error statistics, respectively. In other words, the smaller the observation error variances, the closer the analysis to the observations. The governing equations of the data assimilation scheme are given in the Appendix.
3.2.1. Observation operator, H
The data assimilation scheme relies on the capacity to simulate observations from the current state variables of the model. In our case, H is the observation operator, which estimates the albedo from simulated snow/ice surface layer characteristics (snow grain size, dendricity and sphericity, age, etc.). The albedo value is deduced from the optical grain radius, which is itself computed from the grain shape variables. Our observation operator is therefore composed of two successive elementary operators. They are built from the CROCUS routines.
Because CROCUS is a thermodynamic model, many thresholds are used in its implementation. Consequently we first need to implement a C∞-observation operator, i.e. an observation operator that has a derivative at all orders. This is done by replacing the thresholds by a hyperbolic tangent function. The differences between the original version of H and the class C∞ -version are smaller than 10-6 albedo units, regardless of whether the dendritic or non-dendritic case is chosen. Then a derivation of the code provides the linearized observation operator. From this linearized observation operator, the adjoint operator is implemented. This linear adjoint operator converts a small variation of the albedo values in terms of surface grain characteristics. A gradient test (Eqn (A2)) and verification of scalar product equality (Eqn (A3)) were performed.
3.2.2. Error covariance matrix estimation
A delicate problem in the implementation of a data assimilation scheme is the estimation of error covariance matrices (Reference Desroziers, Brousseau and ChapnikDesroziers and others, 2005). The values of these matrices are of crucial importance, since they determine the respective weights of the observations and the guess in the analysis. The guess error variances, i.e. the diagonal of the error covariance matrix, were computed using the difference between observed snow surface characteristics and CROCUSno assim results during one winter season at Col de Porte (2007/08). CROCUSno assim was then run using disturbed SAFRAN variables as input. The disturbed SAFRAN variables are SAFRAN estimates without assimilation of the observations from AWSs in the mountainous areas. The comparison between the CROCUS output in the case of disturbed meteorological variables and reference SAFRAN variables gives an estimate of the guess error covariances. The guess error covariance matrix is then well defined.
There is more uncertainty in the computation of the observation error covariance matrix, R. This matrix is estimated using the evaluation of accuracy of the albedo retrieval methods presented by Reference Dumont, Sirguey, Arnaud and SixDumont and others (2011). This matrix also needs to take into account the modeling errors of H, especially for the impurity effect. Since only broadband field measurements were available to estimate the error of the remotely sensed albedo, we were only able to determine accurately the broadband observation error variance. The effect of changing the magnitude of R has been investigated. When multiplying R by 10-2, i.e. decreasing the observation error variance to an unrealistically low value, no significant changes were seen in the analysis. By contrast, when R was multiplied by 104, i.e. increasing the observation error variances to an unrealistically high value, the analysis was closer to the guess than in the reference case.
3.2.3. Optimization
The surface layer grains are characterized by a state vector x = |(x 1, x 2, x 3), where x 1 is the dendricity, x 2 the sphericity and x 3 the grain size. This state vector is bounded (Reference Brun, David, Sudul and BrunotBrun and others, 1992, Fig. 1). In order to find the optimal surface grain state vector, x a, i.e. the vector that minimizes the distance both to the observed albedo and to the guess surface grain variables (CROCUS output), the cost function, J(x), defined in Eqn (A4) is minimized. To the classical formulation of the cost function we add a smooth constraint term J c(x) (Eqn (A4)),
where 
. These functions are chosen so as to maintain x inside its variation range.
We use an iterative algorithm to find the optimal state vector. The first gradient descent is performed using ηi = 0, ∀i ∊ 1,2,3, i.e. without any constraint. Then, if the resulting state vector is not in the variation range, ηi are updated increasingly, i.e. the constraints become stronger until the optimal state vector enters its variation range.
Once the optimal state vector has been found, the surface layer grain variables are replaced by the values of x a.
3.3. Strategy for composite snow/ice cases
The assimilation scheme described above is only usable if the albedo and the surface layer characteristics are linked through a non-constant relationship. In the CROCUS model, the ice is simulated with a constant albedo. Consequently, in cases where the observations or the guess surface layer can be classified as ice, the variational scheme is no longer usable.
Thus we had to implement different assimilation strategies, corresponding to four different cases. These cases are defined by comparing remotely sensed albedo observations and modeled surface layer grain variables:
-
Case 1: The observed albedo is larger than the ice albedo and the modeled surface layer is snow. In this case a variational assimilation method is used to modify surface grain variables.
-
Case 2: The observed albedo is larger than the ice albedo and the surface layer is ice. In this case a forcing strategy is used. A layer of new snow is added over the ice surface layer with density and temperature corresponding to the meteorological variables.
-
Case 3: The observed albedo is smaller than the ice albedo and the surface layer is ice. In this case a forcing strategy is used and the ice albedo is replaced by the observed values.
-
Case 4: The observed albedo is smaller than the ice albedo and the surface layer is snow. In this case a forcing strategy is used. The snow layers are removed from the simulated snowpack and the ice albedo is replaced by the observed value.
The ice vs snow albedo threshold for the observation is set to 0.35 for band 1, since it corresponds to the maximum ice albedo value measured on the glacier during summer campaigns 2008 and 2009 (relatively clean ice). The distinction between snow and ice for the simulated surface layer is based on the density threshold of 850 kg m-3 given by CROCUS.
In the following, the methodology described in this section will be referred as CROCUSassim while the reference run, without albedo data assimilation, will be referred to as CROCUSno assim.
4. Results
4.1. Test at one point
The performance of the data assimilation scheme was first evaluated at one point on the Col de Porte site (Reference Brun, Martin, Simon, Gendre and ColeouBrun and others, 1989).
Two types of data assimilation experiments were carried out. First, a CROCUSno assim simulation was used as a reference and the observed albedos were calculated from weekly observed snow characteristics at Col de Porte over two winter seasons (2006/07 and 2007/08). Secondly, a CROCUSno assim simulation with disturbed meteorological input (as for the computation of the guess error covariance matrix) was taken as the reference state and a CROCUSno assim simulated snowpack with reference SAFRAN input gave the observed state vector. This second virtual (i.e. not based on real observations) experiment was carried out for one winter season, 2007/08. Both experiments aimed at a formal validation of the adjoint technique and a better understanding of the behavior of the adjoint model with respect to the grain characteristics.
The analysis error variances for the two experiments described above were evaluated at each assimilation step using the inverse of the Hessian ofthe cost function (Eqn (A8)) (Reference Bouttier and CourtierBouttier and Courtier, 2002). The analysis error variances were systematically smaller than guess error variances for both experiments.
Figure 3 illustrates the results obtained for the first experiment for winter 2007/08. It presents the three components of the state vector (grain variables) and the optical diameter over the winter season 2007/08 (11 November 2007 to 26 April 2008) for CROCUSno assim (solid line) and CROCUSassim (dotted line). Observation and analysis state vectors are also indicated for each snow observation. The average time during winter 2007/08 over which CROCUSno assim and CROCUSassim paths were distinct, was 4.5 days, which gives us an idea of the spin-up time of the model. Figure 3 also shows that the analysed optical diameter is closer to the observations than the guess diameter, expect for 22 January. The analysed grain variables are closer to the observations than the guess ones, except for 1 7, 22, 30 January.
Fig. 3. Temporal evolution of the three components of the state vector for the surface layer (grain variables) of the snowpack during winter 2007/08 at Col de Porte. The solid line represents the CROCUSno assim simulation and the gray dotted line represents the CROCUSassim simulation. The crosses are observations of the snowpack used in the data assimilation scheme. The dark circles are the analysis state vector after each assimilation. The lower chart represents the evolution of the optical diameter which is used in CROCUS to compute the albedo.
We also compared the simulated snow water equivalent (SWE) and snow depth with in situ measurements made over winter season 2007/08. The mean snow-depth rmse with respect to daily measured values is 0.089 m for CROCUSassim and 0.104m for CROCUSnoassim. The mean observed value during the winter is 0.588 m. For SWE, the mean observed value is 255 mm and the rmse is 52 mm w.e. for CROCUSassim and 53 mm w.e. for CROCUSnoassim. The SWE rmse were computed using daily measurements over the winter.
4.2. Glacier de Saint-Sorlin spatially distributed mass balance
We now attempt to evaluate whether it is possible to use CROCUSassim with SAFRAN input and the observed albedo maps to simulate the glacier spatially distributed mass balance on a 100 m × 100 m grid. This section presents the simulation of energy fluxes, albedo and spatially distributed mass balance with CROCUSno assim and CROCUSassim compared with field measurements.
4.2.1. SEB simulation at one point
Figure 4 and Table 2 compare the simulated energy fluxes with the energy fluxes estimated from measurements at the SEB2006 location over 50 days in summer 2006. SEB2006 is located close to AWSabla,2008 in Figure 2. The simulations were done without assimilation of albedo data (CROCUSno assim, black curve) and with assimilation of MODIS data (CROCUSassim, gray curve). The shortwave radiation budget, SWnet, is the main source of error and error variability for the SEB simulation over this period. SWnet explains most of the difference between measurement and simulation up to event 1 in Figure 4 (27 July) and also for event 2 (8 August). Considering the simulated and measured albedo values (not presented here), we note that bias on the albedo is the main contributor to the difference in the shortwave radiation budget.
Fig. 4. Comparison between the measured and simulated SEBs at the location of the SEB measurement station during summer 2006 (close to AWSabla,2008 in Fig. 2). (a) Net shortwave radiation, (b) net longwave radiation, (c) sensible heat flux, (d) latent heat flux and (e) sum of the atmospheric fluxes. The dashed curve shows the measured flux. The black curve is the flux simulated with CROCUSno assim and the gray curve is the flux simulated with CROCUSassim using MODIS data. The vertical bars indicate measurement uncertainties. Events 1, 2 and 3 are indicated for discussion purposes. The measurements are presented by Reference Six, Wagnon, Sicart and VincentSix and others (2009).
Table 2. Comparison between measured and simulated surface daily energy fluxes at location SEB2006 (close to AWSabla2008 in Fig. 2) for 50 days (Fig. 4). SWnet is the shortwave radiation budget, LWnet the longwave radiation budget, H the sensible heat flux and LE the latent heat flux. ΔQ is the SEB expressed as the sum of the four previous fluxes. µaws and σ aws are, respectively, the mean and standard deviation of the daily value over the whole measurement period. m is the mean daily bias between measured and simulated fluxes in case of CROCUSno assim (subscript ‘na’) and CROCUSassim (subscript ‘a’). r is the root-mean-square deviation between daily measured and simulated fluxes. The simulations are done using SAFRAN/CROCUSno assim (without albedo data assimilation) and SAFRAN/CROCUSassim (with assimilation of MODIS data)
The longwave radiation budget does not show any major bias, but the error variability is significant (Table 2). For turbulent fluxes, the difference between measured and simulated fluxes explains most of event 3 on 20 August. During the considered period, the bias on H seems to compensate for the bias on LE. The sum of H and LE is positive for measurements and simulation, except for 3 days. Consequently, they are mainly an energy source for the glacier surface. For CROCUSassim the bias and rmse on the SWnet and SEB are smaller than for CROCUSno assim, but are still significant.
4.2.2. Albedo estimates
Figure 5 shows the evolution of the measured and simulated broadband albedo during summer 2009 at AWSabla,2009. The simulations were run with CROCUSno assim (gray curve), CROCUSassim with MODIS maps (black curve) and CROCUSassim with terrestrial photograph maps (gray crossed curve). The albedo simulated with CROCUSassim is closer to the measured albedo than the albedo simulated with CROCUSno assim. Nevertheless, the agreement between the simulated and the measured albedo is far from perfect and the temporal variability is not well captured. The first day when the surbface is ice (albedo <0.35) is 2 July for the measurements, 12 July for CROCUSassim with photographs, 22 July for CROCUSassim with MODIS and 13 August for CROCUSno assim. After 22 July the albedo simulated with CROCUSassim is within the estimated accuracy of the measurements. During the first period in which the surface is snowy, the simulated albedo is systematically higher than the measurements.
Fig. 5. Measured vs simulated broadband albedo at AWSabla,2009 during summer 2009 (Fig. 2). The dotted curve plots the measurements obtained using the Kipp & Zonen CNR1 device. The vertical bars are an estimate of the uncertainties. The gray curve shows the results of the CROCUSno assim simulation. The black curve shows the results of the CROCUSassim simulation with assimilation of MODIS albedo maps. The gray crossed curve shows the results of the CROCUSassim simulation with assimilation of albedo maps derived from terrestrial photographs. The ice albedo is set at [0.23, 0.16, 0.05] at the beginning of each simulation.
The simulated and observed albedo maps of the glacier were also compared. Figure 6 gives an example of these maps and of the observed albedo for 14 June 2009. Old snow albedo (the lower part of the glacier in the figure) is overestimated by at least 0.05 in the simulation. The albedo map simulated with assimilation of the terrestrial photographs shows more variability than the albedo maps simulated with assimilation of MODIS data.
Fig. 6. Glacier albedo maps: (a) CROCUSno assim; (b) CROCUSassim with assimilation of MODIS and photographs albedo; (c) CROCUSassim from MODIS; and (d) CROCUSassim from photographs. The albedo values are for 14 June 2009. The CROCUS-derived maps were taken just after the assimilation of the albedo observations.
4.2.3. Spatially distributed mass balance
Figures 7 and 8 and Table 3 compare simulated and measured mass balance for each of the five hydrological years at each stake location on the glacier. Over the whole period, 145 annual mass-balance measurements were available on the glacier and were compared with the simulated values. CROCUSnoassim-simulated mass balance has a positive bias of 0.53 m w.e. for these years. The results for CROCUSassim with MODIS data assimilation are not significantly biased and the rmse is 0.48 m w.e. compared with the field measurements. The largest difference is found for years 2002/03 and 2007/08. In addition, the assimilation of the terrestrial photographs improves the results for 2007/08 but has hardly any effect on 2008/09. Lastly, Figure 8 shows that the stakes located in the dense network have the largest rmse.
Fig. 7. Comparison between the measured and simulated mass balance for the five hydrological years without assimilation. The mass balance is cumulated over each hydrological year from the first measurement at the beginning of the ablation period to the last measurement at its end. The simulations were done using CROCUSno assim. The different symbols indicate the different zones of the glacier as shown in Figure 2. The linear regression analysis gives 0.85x + 0.39 with r2 = 0.88
Fig. 8. Comparison between the measured and simulated mass balance for the five hydrological years. The mass balance is cumulated over each hydrological year from the first measurement at the beginning of the ablation period to the last measurement at its end. The simulations were done using CROCUSassim with assimilation of MODIS albedo maps. The different symbols indicate the different zones of the glacier, as shown in Figure 2. The linear regression analysis gives 0.97x + 0.11 with r2 = 0.93.
Table 3. Comparison between measured and simulated annual mass balance for the five hydrological years. µs(σs ) is the mean annual mass balance and its standard deviation (measured values) for all the stakes. mC and rmseC are the bias and rmse when comparing measurements with annual mass balances simulated by CROCUSno assim for all stakes. mC,a and rmseC,a are the bias and rmse when comparing measurements with annual mass balances simulated by CROCUSassim for all stakes. All these results are given for assimilation of MODIS albedo maps, except for values in italics for years 2007/08 and 2008/09, which are for assimilation of albedo maps derived from terrestrial photographs. ‘Stakes’ gives the number of stakes used for each year
5. Discussion
5.1. General performance of the method
The results described for the data assimilation scheme show that albedo assimilation improves the snowpack simulation. The analysis error variances are systematically smaller than the guess error variances. During winter 2007/08 at Col de Porte (Fig. 3), there are three cases for which the analysis state vector is not closer to the observations than the guess state vector (17, 22 and 30 January). It snowed on 22 January, so the modifications induced by the assimilation were immediately erased by new snow. On 17 and 30 January the analysed optical diameter is closer to the observations than the guess, but this is not the case for the state vector. This illustrates the fact that we are trying to advect information that is contained in one variable (albedo or optical diameter) to three grain variables. In general, the information seems to be correctly reported to the grain variables but it is not systematic. In both cases, the grain size decreased but there was no effect on the dendricity. This emphasizes the need for a more physical albedo parameterization to further improve the impact of the assimilation. The SWE and snow depth estimates during winter 2007/08 are closer to field measurements for CROCUSassim than for CROCUSn0 assim. The assimilation impact is small since CROCUSno assim was already very close to measurements. In addition, Figure 3 shows that the spin-up time of the model is a few days. This means that the period of the assimilation has to be short (no more than a few days) to have a significant impact on the simulation (at least during winter at Col de Porte). In this study, the albedo observations are available for a much longer period (Table 1). Nevertheless, using a mixed assimilation/forcing algorithm improves the performance of the algorithm despite the large temporal resolution of available observations.
Comparison of the CROCUSno assim simulated and measured SEB at one point shows that the main error can be attributed to the simulated value of albedo and that errors in the other fluxes can be significant. Improving the shortwave radiation budget does not allow complete correction of the whole SEB and can unbalance some errors which would have otherwise been compensated for by the present features and performance characteristics of the CROCUS model. Error in turbulent fluxes can be attributed to the parameterization itself, in particular the input wind speed or the selected roughness length, which may not be representative of the surface at the beginning of July. During this period, the measurement station was located not far from the snowline. In the neighborhood of the snowline, snow/ice patches and meltwater channels appear and the roughness lengths can vary significantly.
Nevertheless, the comparison between the simulated and measured albedo shows that assimilation improves the simulation. In Figure 5, first ice appearance is more consistent for CROCUSassim than for CROCUSno assim. This underestimated melt rate (or overestimated accumulation) at the beginning of the ablation season for CROCUSno assim may be attributed to the use of a mean correction factor for precipitation. The simulated precipitation field may indeed induce an over-accumulation at the beginning of the ablation season for this year. Assimilating the albedo data corrects this effect without the need for in situ precipitation measurements. Old-melting-snow albedo is nonetheless overestimated in the CROCUSassim simulation. This may be due to the setting of a minimum value in the albedo parameterization. This parameterization has been developed for the winter seasonal snowpack and may not be able to represent old dirty wet snow and firn. In Figure 5, the temporal variability of the measured albedo is poorly reproduced on the simulated albedo. This might be explained by the fact that there is only about one available albedo map per week, which is not temporally sufficient information considering the spin-up time of the model. The remaining difference between the observed and simulated snow albedo in June 2009 cannot be attributed to a too-small value of the observation covariance matrix, as explained above (Section 3.2.2).
The assimilation of albedo observations in the model largely improved the simulation of the spatially distributed mass balance for the five studied hydrological years. It provides an unbiased estimation of spatially distributed mass balance with a rmse <0.5 m w.e., whereas the simulation without assimilation shows a positive bias of 0.53 mw.e. This positive bias can be attributed to an albedo positive feedback mechanism. At least for years 2000/01, 2007/08 and 2008/09, the albedo is overestimated at the beginning of the ablation season, leading to underestimation of the amount of absorbed shortwave radiation in the snowpack. The snowmelt is therefore also underestimated, which, in turn, amplifies the overestimation of the albedo. This effect is corrected when assimilating albedo observations for years 2000/01 and 2008/09. For 2007/08, MODIS data spatial resolution seems to be insufficient and only assimilation of terrestrial photographs allows correct simulation of the spatially distributed mass balance, through a more precise implicit knowledge of the snowline position captured by the assimilation scheme. Part of this positive bias and albedo positive feedback might be due to an insufficient description of the effect of impurities on the albedo parameterization implemented in the H operator included in the CROCUS snow model. Assimilating albedo observation helps to reduce the effect of this shortcoming. However, we believe that using a more accurate albedo parameterization (e.g. as described by Reference Gardner and SharpGardner and Sharp, 2010) would improve the simulated mass balance but would require additional information about the impurity content and its evolution. This impurity knowledge is implicitly contained in the albedo observations and would be transferred to the simulated snowpack using such an assimilation technique.
Year 2002/03 is an exception, since assimilating albedo slightly decreases the accuracy of the simulation. A comparison between the meteorological variables for each summer of the five years, shows that the wind speed and temperature values were maximum and nebulosity minimum for 2002/03. These facts can explain the overestimated melt. First, this meteorological situation implies that the turbulent fluxes have a greater effect. For example, during a similar period (hot and sunny) in summer 2006, turbulent fluxes derived from measurements reached 21% of the sum of the atmospheric fluxes at SEB2006 (Reference Six, Wagnon, Sicart and VincentSix and others, 2009). Consequently, errors in the estimation of these fluxes (Table 2) may have a greater weight on the SEB. Secondly, as discussed above, longwave radiation estimates for SAFRAN are overestimated in the case of low nebulosity. Although a correction has been applied, it may not be sufficient for this case and can also lead to melt overestimation.
The results shown in Figure 8 indicate that mass-balance estimates on the dense network have the largest rmse. The dense network is characterized by rugged topography and stakes are located ˜40m apart. The model grid is 100 m × 100 m. Thus, in this area of the glacier, the mass-balance spatial variability depicted by the stakes is smaller than that which can be reached by the model.
Furthermore, using albedo maps with 10m instead of 250 m resolution improves the simulation for year 2007/08 but has hardly any effect on year 2008/09. The simulation for year 2008/09 was already good using MODIS data (bias -0.05 m w.e.) and we believe that the accuracy limit of the model had already been reached. Adding additional spatial information about the albedo value does not improve the simulation in such a case.
5.2. Precipitation correction
A mean precipitation correction factor may be used for the whole glacier, instead of a detailed map, when applying the method to a large set of glaciers. We ran the mass- balance simulation for the five hydrological years using a mean correction factor of 1.64 instead of the map. The results show that the estimation of the annual mass balance is still unbiased (0.06 m w.e). The rmse with respect to the 145 measured mass-balance values is now 0.54m w.e., instead of 0.48 m w.e. with the map. Using the mean factor slightly reduced the accuracy of the estimates but the simulation of the spatially distributed mass balance was still reasonable. This effect must also be investigated on larger glaciers where the mass-balance variations over the glacier are higher. In that case, it may be necessary to use different multiplication factors over the glacier.
In addition, we used the same multiplication factor for winter and summer precipitation, ignoring the phase of the precipitation. The correction factor could be different and closer to that for liquid precipitation (Reference GottardiGottardi, 2009). Nonetheless, the rain heat flux is known to be of minor importance for the total SEB on glaciers in the Alps (Reference HockHock, 2005).
5.3. Roughness lengths
Finally, two constant roughness lengths were used in this study, one for ice and one for snow. As discussed in Section 3, these values cannot be fully considered as physical values but as effective ones. The spatial and temporal variability of roughness lengths is important at the glacier surface. Measured values vary from 0.1 mm forfresh snow (Reference Brock, Willis and SharpBrock and others, 2006) to 30 mm for penitents (Reference Wagnon, Ribstein, Francou and PouyaudWagnon and others, 1999). The surface macro-roughness (snow/ice patches, meltwater channels and crevices) has the largest influence on the value of the roughness length. In this study, we evaluate the influence of roughness length variations on the simulated melt rate. An illustration is given in Figure 9, which presents the impact of the variation in snow and ice roughness lengths on the simulated mass balance at one point. The results show that changing the snow roughness length has more effect on the simulated mass balance than changing the ice roughness length. Indeed when the snow roughness length varies from 1 to 10 mm, the simulated mass balance at the end of the summer decreases 20%, whereas the same variation in the ice roughness length only causes a decrease of 9% in the simulated mass balance. The snow roughness length affects the melting rate of the winter snowpack and, consequently, determines snowline position. It could therefore be useful in further work to use different roughness lengths, depending on the snow age and depth.
Fig. 9. Effect of variation in snow and ice roughness lengths on the simulated mass balance. These simulations were run with CROCUSno assim at the location of SEB2006 (close to AWSabla,2008 in Fig. 2) for 2005/06. The solid and dotted curves show the effect of snow roughness length variations. The dashed curves present the effect of ice roughness length variations.
6. Conclusions and Perspectives
Mesoscale meteorological variables used with the CROCUS snow model and assimilation of albedo observations from MODIS and/or terrestrial photographs provide a reasonable estimation of the spatially distributed mass balance of a glacier in the French Alps. The estimates are unbiased, with a rmse <0.5 m w.e. over the five hydrological years. This study is the first attempt to assimilate albedo observations in order to numerically simulate the high-resolution distribution of the SEB of an alpine glacier. It demonstrates the potential contribution of remote-sensing data to the improvement of snowpack and glacier mass-balance simulation.
The method described here can be used on a large set of glaciers to estimate the spatial distribution of temperate glacier mass balance simply by combining remotely sensed albedo observations with a numerical snow energy-balance model. It is applicable to the seasonal snow cover on larger areas. The present limitations of the method are the accuracy of the knowledge of precipitation fields and the tuning of the roughness lengths. Finer precipitation fields will soon be routinely available (Reference SeitySeity and others, 2011). Concerning the snow roughness lengths, research is still in progress, but this study provides a first guess for their determination. For temperate glaciers, we believe the roughness length can be tuned using only one value of the specific annual mass balance. For the seasonal snow cover, realistic profiles of the snowpack are nowadays provided by CROCUS, using a unique roughness length over the Alps which can be used as a first guess (Reference Durand, Giraud, Brun, Merindol and MartinDurand and others, 1999). Another limitation is that the albedo parameterization used in this study has no explicit formulation for the snow impurity content and uses non-continuous grain variables. Work is in progress to implement a more physical parameterization of the albedo in CROCUS, and we believe such parameterization would further improve the impact of the assimilation scheme.
Using an energy-balance approach gives better quantification of the different fluxes involved in the snowpack or glacier mass balance. In addition, such a variational approach is consistent with theoretical sensitivity analysis to determine, for example, which parameters are the most important in the model. A possible application might be the determination of the most important solar spectral band in the model formulation.
Furthermore, this variational method could be used to take into account effects that are not yet accurately parameterized, such as snow impurity content and snowdrifts. Impurity effects on snow albedo are mainly found in the first spectral band (Reference WarrenWarren, 1982). Assimilating this spectral band introduces information about the impurity content of the observed snowpack. Consequently, impurity content can be considered as a variable in the snowpack simulation model. Similarly, it is known that snowdrifting by wind changes snow albedo (Reference Domine, Taillandier, Cabanes, Douglas and SturmDomine and others, 2009). Changes in albedo due to wind can be detected through remote-sensing data and therefore could be included in the snowpack simulation model (Reference CorripioCorripio, 2004). More generally, such methods have strong potential in ice and snow modeling, in that they add various sources of information to the simulation.
Lastly, albedo was a first choice for testing the effect of variational data assimilation in a snowpack model. Other remote-sensing variables that provide information on the surface or on the whole snowpack, such as snow water equivalent, depth, surface temperature and specific surface area, could be assimilated in the snow model using a similar scheme and further improve the results.
7. Acknowledgements
We thank the LEFE-ASSIM program (ADAC project) for funding the main part of this study, and the GLACIOCLIM observatory and ANR TAG for providing the field data. We also thank Yves Lejeune, Jean-Marie Willemet and Bernard Lesaffre for their help at Col de Porte and with CROCUS, and Patrick Wagnon, Jean-Emmanuel Sicart, Christian Vincent and Marc Bocquet for helpful comments. We are grateful to Eric Brun, Samuel Morin, Peter Kuipers Munneke and an anonymous referee for their helpful and relevant comments on the first draft.
Appendix: Variational Data Assimilation Method, Definitions
This Appendix gives definitions that are essential to understand the methodology (Reference Bouttier and CourtierBouttier and Courtier, 2002).
Let x = (x 1, x 2, x 3) be the surface grain variables vector and y = (α1, α2, α3) be the vector composed of the three spectral bands of albedo.
Let β be the model space (grain variables) and O the observation space (albedo).
A.1. Observation operators
We define H, the observation operator, as y = H[x]. In the absence of any modeling error, this operator would generate the albedo value, y, if both observations and state vectors were perfect.
The linearized observation operator at x 0, H, is defined as
Consequently, H is the true linearized operator of H if
The adjoint observation operator, H× , is defined as
A.2. Errors and error covariances
If ε represents the error vector, then 
is the error covariance matrix. x
b is the guess state vector given as a priori information by the model, x
t is the true state vector, y is the vector of observations and x
a is the analysis model state (after assimilation).
∊b = x b - x t is the guess error with covariances B. It is the estimation error of the guess state and ∊0 = y - H[xt] is the observation error with covariances R. It contains the error in the observation process, errors due to the scale difference between observations and the model and due to the design of H and discretization errors that prevent x t from being a true state
∊a = x a - x t is the analysis error with covariances A. It is the estimation error of the analysis state which has to be minimized.
A.3. Optimization
Using the variational method implies that we define a cost function, J, in the model space. The analysis state vector, x a, is then defined as the state vector that minimizes this cost function:
Jb(x) = ||x — x b||2β is the guess term, J o(x) = ||y — H[x]2O the observation term and J c(x) the constraint term that prevents x a from escaping from its validity domain.
From Eqn (A4), we can calculate ∇J and ∇2J:
To estimate x a, we use a popular gradient descent algorithm:
Note thatthe analysis error covariance matrix, A, is directly linked to the Hessian of J:
