<figure><img src="image_1.png" /></figure>

HAL
open science

# Variational assimilation of albedo in a snowpack model and reconstruction of the spatial mass balance distribution of an alpine glacier

M. Dumont, Yves Durand, Y. Arnaud, D. Six

► To cite this version:

M. Dumont, Yves Durand, Y. Arnaud, D. Six. Variational assimilation of albedo in a snowpack model and reconstruction of the spatial mass balance distribution of an alpine glacier. Journal of Glaciology, 2011, 58 (207), pp.151-164. 10.3189/2012JoG11J163. insu-00648256

HAL Id: insu-00648256

https://insu.hal.science/insu-00648256v1

Submitted on 12 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Variational assimilation of albedo in a snowpack model and reconstruction of the spatial mass-balance distribution of an alpine glacier

Marie DUMONT,¹,² Yves DURAND,² Yves ARNAUD,¹,³ Delphine SIX¹

¹Université Joseph Fourier Grenoble-CNRS, LGGE UMR 5183, Grenoble, France
E-mail: marie.dumont@meteo.fr

²Météo-France-CNRS, CNRM-GAME URA 1357, CEN, Grenoble, France

³Université Joseph Fourier Grenoble-CNRS-IRD-Grenoble INP, LTHE UMR 5564, Grenoble, France

**ABSTRACT.** Accurate knowledge of the spatial distribution of the mass balance of temperate glaciers is essential for a better understanding of the physical processes controlling the mass balance and for the monitoring of water resources. In relation to albedo variations, the shortwave radiation budget is a controlling variable of the surface energy balance of glaciers. Remotely sensed albedo observations are here assimilated in a snowpack model to improve the modeling of the spatial distribution of the glacier mass balance. The albedo observations are integrated in the snowpack simulation using a variational data assimilation scheme that modifies the surface grain conditions. The study shows that mesoscale meteorological variables and MODIS-derived albedo maps can be used to obtain a good reconstruction of the annual mass balance on Glacier de Saint-Sorlin, French Alps, on a 100 m × 100 m grid. Five hydrological years within the 2000–10 decade are tested. The accuracy of the method is estimated from comparison with field measurements. Sensitivity to roughness lengths and winter precipitation fields is investigated. Results demonstrate the potential contribution of remote-sensing data and variational data assimilation to further improve the understanding and monitoring of the mass balance of snowpacks and temperate glaciers.

## 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. Courtier and others, 1998; Uppala and others, 2005).

Recently, data assimilation has begun to be used in glaciology. For instance, Arthern 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. Andreadis and Lettenmaier, 2006; Slater and Clark, 2006; Durand and others, 2008; De Lannoy and others, 2010; Toure 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.

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 (Sicart and others, 2008). In addition, variations of the solar radiation budget (shortwave radiation balance) are governed by variations of surface albedo (Sicart 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 (Oerlemans and Knap, 1998). Albedo can be observed using satellite sensors (Dozier 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. Hock, 2005; Andersen and others, 2010; Giesen and Oerlemans, 2010). Mass-balance models vary from simple temperature-index models to so-called energy-balance models (Hock, 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 Essery and others (1999), Essery and Etchevers (2004) and Etchevers and others (2004). All these authors have pointed out the major influence of the albedo parameterization on the final mass-balance estimates. Klok and Oerlemans (2002) and Hock and Holmgren (2005)
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Simplified schematic view of the methodology used to reconstruct the spatial distribution of the glacier mass balance.</figcaption></figure>

used models to simulate spatially distributed mass balances with two different albedo parameterizations and, once again, emphasized the crucial role of this variable.

Gerbaux and others (2005), Lejeune and others (2007) and Lejeune (2009) developed a spatially distributed version of the CROCUS snow model (Brun and others, 1989, 1992). Gerbaux 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 (Durand and others, 1993). SAFRAN disaggregates large-scale meteorological analysis and observations in the French Alps. Lejeune and others (2007) and Lejeune (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 (Pedersen and Winther, 2005; Gardner 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 Gerbaux 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 (Bouttier 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 Météo-France located at 1326 m a.s.l. in the Chartreuse area, near Grenoble, France (Brun and others, 1989, 1992). 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 ~3 km². 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 (Vincent and others, 2000; Vincent, 2002). Numerous dynamic modeling and climate studies have been conducted on this glacier (e.g. Le Meur and Vincent, 2003; Le Meur and others, 2007). Detailed information on this site can be found at http://www-lgge.ujf-grenoble.fr/ServiceObs.

Two types of measurements available on this glacier are used in this study:

**Automatic weather stations**

A permanent AWS (AWS<sub>moraine</sub>) 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 (SEB<sub>2006</sub>) was set up on the glacier during summer 2006 (Sicart and others, 2008; Six and others, 2009). The SEB<sub>2006</sub> 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 (AWS<sub>abla,2008</sub> and AWS<sub>accu,2008</sub> during summer 2008 and AWS<sub>abla,2009</sub> 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
<figure><img src="image_4.png" /><figcaption>**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, SEB<sub>2006</sub>, was located near AWS<sub>abla,2008</sub></figcaption></figure>

(Cogley 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. (Cogley 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 (Dumont and others, 2011) and MODIS data, based on the work of Sirguey 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 (Dumont 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.

## 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,

<table><caption>**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) (Dumont and others, 2011)</caption><thead><tr><td>Hydrological year</td><td>MODIS albedo maps</td><td>Photograph albedo maps</td></tr></thead><tbody><tr><td>2000/01</td><td>20</td><td></td></tr><tr><td>2002/03</td><td>17</td><td></td></tr><tr><td>2005/06</td><td>15</td><td></td></tr><tr><td>2007/08</td><td>16</td><td>11</td></tr><tr><td>2008/09</td><td>20</td><td>19</td></tr></tbody></table>
nebulosity (global cloudiness) and longwave downward flux (Durand 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 km². 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 Gerbaux 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 AWS<sub>moraine</sub> 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<sup>-2</sup> 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 Gerbaux 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 8 cm 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, Gerbaux and others (2005) found a factor of 1.5 for the whole glacier 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 (Gerbaux 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 Brun and others (1989, 1992) 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 (Brun 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 (Martin 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$ (Martin and Lejeune, 1998). This assumption is not invalidated by Andreas (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 Martin (1975), Greuell and Smeets (2001) and Brock 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 (Brun and others, 1992). The optical grain size is computed from grain variables (Willemet, 2008). The albedo is split into three spectral bands, $\alpha_1$: [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, $\alpha_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 (Brun 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. Paterson (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, $\alpha_{1,ice} = 0.45$, $\alpha_{2,ice} = 0.30$ and $\alpha_{3,ice} = 0.1$ (Brun and others, 1992). On Glacier de Saint-Sorlin, Gerbaux 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 Wagnon and others (2009).

## 3.2. Variational assimilation of snow albedo data

The observed albedo is assimilated in CROCUS using a variational data assimilation scheme (Bouttier 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 error 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, $\mathcal{H}$

The data assimilation scheme relies on the capacity to simulate observations from the current state variables of the model. In our case, $\mathcal{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^\infty$-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 $\mathcal{H}$ and the class $C^\infty$-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 (Desroziers 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 CROCUS$_{no assim}$ results during one winter season at Col de Porte (2007/08). CROCUS$_{no 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, $\mathbf{R}$. This matrix is estimated using the evaluation of accuracy of the albedo retrieval methods presented by Dumont and others (2011). This matrix also needs to take into account the modeling errors of $\mathcal{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 $\mathbf{R}$ has been investigated. When multiplying $\mathbf{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 $\mathbf{R}$ was multiplied by $10^4$, 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 $\mathbf{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 (Brun and others, 1992, fig. 1). In order to find the optimal surface grain state vector, $\mathbf{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(\mathbf{x})$, defined in Eqn (A4) is minimized. To the classical formulation of the cost function we add a smooth constraint term $J_c(\mathbf{x})$ (Eqn (A4)),

$$J_c(\mathbf{x}) = \eta_1 J_1(x_1) + \eta_2 J_2(x_2) + \eta_3 J_3(x_3) \quad (1)$$

where $J_1(x_1) = \frac{1}{2}\{\tanh[K_1(x_1 + 0.1)] + 1\}x_1^2$, $J_2(x_2) = \frac{1}{2}(\frac{x_2-50}{50})^8$ and $J_3(x_3) = \frac{1}{2}(\frac{x_3-0.00245}{0.00245})^8$. These functions are chosen so as to maintain $\mathbf{x}$ inside its variation range.

We use an iterative algorithm to find the optimal state vector. The first gradient descent is performed using $\eta_i = 0, \forall i \in \{1, 2, 3\}$, i.e. without any constraint. Then, if the resulting state vector is not in the variation range, $\eta_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 $\mathbf{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.
<figure><img src="image_5.png" /><figcaption>**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 CROCUS<sub>no assim</sub> simulation and the gray dotted line represents the CROCUS<sub>assim</sub> 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.</figcaption></figure>

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⁻³ given by CROCUS.

In the following, the methodology described in this section will be referred as CROCUS<sub>assim</sub> while the reference run, without albedo data assimilation, will be referred to as CROCUS<sub>no assim</sub>.

# 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 (Brun and others, 1989).

Two types of data assimilation experiments were carried out. First, a CROCUS<sub>no assim</sub> 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 CROCUS<sub>no assim</sub> simulation with disturbed meteorological input (as for the computation of the guess error covariance matrix) was taken as the reference state and a CROCUS<sub>no assim</sub> 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 of the cost function (Eqn (A8)) (Bouttier 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 CROCUS<sub>no assim</sub> (solid line) and CROCUS<sub>assim</sub> (dotted line). Observation and analysis state vectors are also indicated for each snow observation. The average time during winter 2007/08 over which CROCUS<sub>no assim</sub> and CROCUS<sub>assim</sub> 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 17, 22, 30 January.
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Comparison between the measured and simulated SEBs at the location of the SEB measurement station during summer 2006 (close to AWS<sub>abla,2008</sub> 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 CROCUS<sub>no assim</sub> and the gray curve is the flux simulated with CROCUS<sub>assim</sub> using MODIS data. The vertical bars indicate measurement uncertainties. Events 1, 2 and 3 are indicated for discussion purposes. The measurements are presented by Six and others (2009).</figcaption></figure>

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 CROCUS<sub>assim</sub> and 0.104 m for CROCUS<sub>no assim</sub>. 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 CROCUS<sub>assim</sub> and 53 mm w.e. for CROCUS<sub>no assim</sub>. 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 CROCUS<sub>assim</sub> 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 CROCUS<sub>no assim</sub> and CROCUS<sub>assim</sub> 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 SEB<sub>2006</sub> location over 50 days in summer 2006. SEB<sub>2006</sub> is located close to AWS<sub>abla,2008</sub> in Figure 2. The simulations were done without assimilation of albedo data (CROCUS<sub>no assim</sub>, black curve) and with assimilation of MODIS data (CROCUS<sub>assim</sub>, gray curve). The shortwave radiation budget, SW<sub>net</sub>, is the main source of error and error variability for the SEB simulation over this period. SW<sub>net</sub> 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.

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
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Measured vs simulated broadband albedo at AWS<sub>abla,2009</sub> 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 CROCUS<sub>no assim</sub> simulation. The black curve shows the results of the CROCUS<sub>assim</sub> simulation with assimilation of MODIS albedo maps. The gray crossed curve shows the results of the CROCUS<sub>assim</sub> 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.</figcaption></figure>

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 CROCUS<sub>assim</sub> the bias and rmse on the SW<sub>net</sub> and SEB are smaller than for CROCUS<sub>no assim</sub>, 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 AWS<sub>abla,2009</sub>. The simulations were run with CROCUS<sub>no assim</sub> (gray curve), CROCUS<sub>assim</sub> with MODIS maps (black curve) and CROCUS<sub>assim</sub> with terrestrial photograph maps (gray crossed curve). The albedo simulated with CROCUS<sub>assim</sub> is closer

to the measured albedo than the albedo simulated with CROCUS<sub>no assim</sub>. 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 surface is ice (albedo <0.35) is 2 July for the measurements, 12 July for CROCUS<sub>assim</sub> with photographs, 22 July for CROCUS<sub>assim</sub> with MODIS and 13 August for CROCUS<sub>no assim</sub>. After 22 July the albedo simulated with CROCUS<sub>assim</sub> 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.

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.

<table><caption>**Table 2.** Comparison between measured and simulated surface daily energy fluxes at location SEB<sub>2006</sub> (close to AWS<sub>abla,2008</sub> in Fig. 2) for 50 days (Fig. 4). SW<sub>net</sub> is the shortwave radiation budget, LW<sub>net</sub> 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. μ<sub>AWS</sub> and σ<sub>AWS</sub> 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 CROCUS<sub>no assim</sub> (subscript ‘na’) and CROCUS<sub>assim</sub> (subscript ‘a’). *r* is the root-mean-square deviation between daily measured and simulated fluxes. The simulations are done using SAFRAN/CROCUS<sub>no assim</sub> (without albedo data assimilation) and SAFRAN/CROCUS<sub>assim</sub> (with assimilation of MODIS data)</caption><thead><tr><th>Flux</th><th>μ<sub>AWS</sub><br>W m<sup>-2</sup></th><th>σ<sub>AWS</sub><br>W m<sup>-2</sup></th><th>m<sub>na</sub><br>W m<sup>-2</sup></th><th>r<sub>na</sub><br>W m<sup>-2</sup></th><th>m<sub>a</sub><br>W m<sup>-2</sup></th><th>r<sub>a</sub><br>W m<sup>-2</sup></th></tr></thead><tbody><tr><td>SW<sub>net</sub></td><td>161.5</td><td>59.3</td><td>-30.6</td><td>59.3</td><td>-9.5</td><td>50.6</td></tr><tr><td>LW<sub>net</sub></td><td>-32.9</td><td>21.0</td><td>3.2</td><td>26.4</td><td>3.2</td><td>26.4</td></tr><tr><td>H</td><td>33.0</td><td>29.14</td><td>-5.5</td><td>28.9</td><td>-5.1</td><td>28.7</td></tr><tr><td>LE</td><td>-7.9</td><td>19.7</td><td>20.2</td><td>26.6</td><td>18.3</td><td>25.1</td></tr><tr><td>ΔQ</td><td>154</td><td>64.8</td><td>-12.7</td><td>48.1</td><td>6.9</td><td>46.1</td></tr></tbody></table>

### 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. CROCUS<sub>no assim</sub>-simulated mass balance has a positive bias of 0.53 m w.e. for these years. The results for CROCUS<sub>assim</sub> 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.
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Glacier albedo maps: (a) CROCUS<sub>no assim</sub>; (b) CROCUS<sub>assim</sub> with assimilation of MODIS and photographs albedo; (c) CROCUS<sub>assim</sub> from MODIS; and (d) CROCUS<sub>assim</sub> from photographs. The albedo values are for 14 June 2009. The CROCUS-derived maps were taken just after the assimilation of the albedo observations.</figcaption></figure>

# 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 advance 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 CROCUS<sub>assim</sub> than for CROCUS<sub>no assim</sub>. The assimilation impact is small since CROCUS<sub>no assim</sub> 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

<table><caption>**Table 3.** Comparison between measured and simulated annual mass balance for the five hydrological years. $\mu_s(\sigma_s)$ is the mean annual mass balance and its standard deviation (measured values) for all the stakes. $m_C$ and $rmse_C$ are the bias and rmse when comparing measurements with annual mass balances simulated by CROCUS<sub>no assim</sub> for all stakes. $m_{C,a}$ and $rmse_{C,a}$ are the bias and rmse when comparing measurements with annual mass balances simulated by CROCUS<sub>assim</sub> 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</caption><thead><tr><th>Year</th><th>μ<sub>s</sub>(σ<sub>s</sub>)<br>m w.e.</th><th>m<sub>C</sub><br>m w.e.</th><th>rmse<sub>C</sub><br>m w.e.</th><th>m<sub>C,a</sub><br>m w.e.</th><th>rmse<sub>C,a</sub><br>m w.e.</th><th>Stakes</th></tr></thead><tbody><tr><td>2000/01</td><td>0.05(0.82)</td><td>0.26</td><td>0.48</td><td>-0.00</td><td>0.44</td><td>20</td></tr><tr><td>2002/03</td><td>-3.37(0.69)</td><td>-0.42</td><td>0.68</td><td>-0.54</td><td>0.66</td><td>19</td></tr><tr><td>2005/06</td><td>-1.87(0.92)</td><td>0.16</td><td>0.46</td><td>-0.02</td><td>0.50</td><td>36</td></tr><tr><td>2007/08</td><td>-2.23(0.71)</td><td>1.25</td><td>1.14</td><td>0.40<br>0.13</td><td>0.75<br>0.69</td><td>33</td></tr><tr><td>2008/09</td><td>-3.53(0.68)</td><td>0.70</td><td>0.83</td><td>-0.05<br>0.04</td><td>0.41<br>0.42</td><td>35</td></tr><tr><td>Total</td><td></td><td>0.53</td><td>0.81</td><td>0.07</td><td>0.48</td><td>143</td></tr></tbody></table>
<figure><img src="image_9.png" /></figure>

**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 CROCUS<sub>no assim</sub>. 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 $r^2 = 0.88$.

<figure><img src="image_10.png" /></figure>

**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 CROCUS<sub>assim</sub> 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 $r^2 = 0.93$.

of the algorithm despite the large temporal resolution of available observations.

Comparison of the CROCUS<sub>no assim</sub> 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 CROCUS<sub>assim</sub> than for CROCUS<sub>no assim</sub>. This underestimated melt rate (or overestimated accumulation) at the beginning of the ablation season for CROCUS<sub>no assim</sub> 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 CROCUS<sub>assim</sub> 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 m w.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 Gardner 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 SEB<sub>2006</sub> (Six 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 ~40 m 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 10 m 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.54 m 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 (Gottardi, 2009).
Nonetheless, the rain heat flux is known to be of minor
importance for the total SEB on glaciers in the Alps
(Hock, 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 for fresh snow (Brock and
others, 2006) to 30 mm for penitents (Wagnon 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.

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 (Seity 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 (Durand 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.
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Effect of variation in snow and ice roughness lengths on the simulated mass balance. These simulations were run with CROCUS<sub>no assim</sub> at the location of SEB<sub>2006</sub> (close to AWS<sub>abla,2008</sub> 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.</figcaption></figure>

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 par-
ameterized, such as snow impurity content and snowdrifts.
Impurity effects on snow albedo are mainly found in the first
spectral band (Warren, 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 (Domine 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 (Corripio, 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.

REFERENCES

Andersen ML and 14 others (2010) Spatial and temporal melt variability at Helheim Glacier, East Greenland, and its effect on ice dynamics. *J. Geophys. Res.*, **115**(F4), F04041 (doi: 10.1029/2010JF001760)

Andreadis KM and Lettenmaier DP (2006) Assimilating remotely sensed snow observations into a macroscale hydrology model. *Adv. Water Resour.*, **29**(6), 872–886

Andreas EL (1987) A theory for the scalar roughness and the scalar transfer coefficients over snow and sea ice. *Bound.-Layer Meteorol.*, **38**(1-2), 159-184

Arthern RJ and Gudmundsson GH (2010) Initialization of ice-sheet forecasts viewed as an inverse Robin problem. *J. Glaciol.*, **56**(197), 527–533

Bouttier F and Courtier P (2002) *Data assimilation concepts and methods, March 1999*. European Centre for Medium-Range Weather Forecasts, Reading (Meteorological Training Course Lecture Notes)

Brock BW, Willis IC and Sharp MJ (2006) Measurement and parameterization of aerodynamic roughness length variations at Haut Glacier d'Arolla, Switzerland. *J. Glaciol.*, **52**(177), 281–297

Brun E, Martin E, Simon V, Gendre C and Coléou C (1989) An energy and mass model of snow cover suitable for operational avalanche forecasting. *J. Glaciol.*, **35**(121), 333–342
Brun E, David P, Sudul M and Brunot G (1992) A numerical model to simulate snow-cover stratigraphy for operational avalanche forecasting. *J. Glaciol.*, **38**(128), 13–22

Cogley JG and Adams WP (1998) Mass balance of glaciers other than the ice sheets. *J. Glaciol.*, **44**(147), 315–325

Corripio J (2004) Snow surface albedo estimation using terrestrial photography. *Int. J. Remote Sens.*, **25**(24), 5705–5729

Courtier P and others (1998) The ECMWF implementation of three- dimensional variational assimilation (3D-Var). I: Formulation. *Q. J. R. Meteorol. Soc.*, **124**(550), 1783–1807

De Lannoy GJM, Reichle RH, Houser PR, Arsenault KR, Verhoest NEC and Pauwels VRN (2010) Satellite-scale snow water equivalent assimilation into a high-resolution land surface model. *J. Hydromet.*, **11**(2), 352–369

Desroziers G, Brousseau P and Chapnik B (2005) Use of randomization to diagnose the impact of observations on analyses and forecasts. *Q. J. R. Meteorol. Soc.*, **131**(611), 2821–2837

Domine F, Taillandier A-S, Cabanes A, Douglas TA and Sturm M (2009) Three examples where the specific surface area of snow increased over time. *Cryosphere*, 3(1), 31–39

Dozier J, Green RO, Nolin AW and Painter TH (2009) Interpretation of snow properties from imaging spectrometry. *Remote Sens. Environ.*, **113**, Suppl. 1, S25–S37

Dumont M, Sirguey P, Arnaud Y and Six D (2011) Monitoring spatial and temporal variations of surface albedo on Saint Sorlin Glacier (French Alps) using terrestrial photography. *Cryosphere*, 5(3), 759–771

Durand Y, Brun E, Mérindol L, Guyomarc'h G, Lesaffre B and Martin E (1993) A meteorological estimation of relevant parameters for snow models. *Ann. Glaciol.*, **18**, 65–71

Durand Y, Giraud G, Brun E, Mérindol L and Martin E (1999) A computer-based system simulating snowpack structures as a tool for regional avalanche forecasting. *J. Glaciol.*, **45**(151), 469–484

Durand M, Molotch NP and Margulis SA (2008) A Bayesian approach to snow water equivalent reconstruction. *J. Geophys. Res.*, **113**(D20), D20117 (doi: 10.1029/2008JD009894)

Essery R and Etchevers P (2004) Parameter sensitivity in simulations of snowmelt. *J. Geophys. Res.*, **109**(D20), D20111 (doi: 10.1029/2004JD005036)

Essery R, Martin E, Douville H, Fernández A and Brun E (1999) A comparison of four snow models using observations from an alpine site. *Climate Dyn.*, **15**(8), 583–593

Etchevers P and 22 others (2004) Validation of the energy budget of an alpine snowpack simulated by several snow models (SnowMIP project). *Ann. Glaciol.*, **38**, 150–158

Gardner AS and Sharp MJ (2010) A review of snow and ice albedo and the development of a new physically based broadband albedo parameterization. *J. Geophys. Res.*, **115**(F1), F01009 (doi: 10.1029/2009JF001444)

Gerbaux M, Genthon C, Etchevers P, Vincent C and Dedieu JP (2005) Surface mass balance of glaciers in the French Alps: distributed modeling and sensitivity to climate change. *J. Glaciol.*, **51**(175), 561–572

Giesen RH and Oerlemans J (2010) Response of the ice cap Hardangerjøkulen in southern Norway to the 20th and 21st century climates. *Cryosphere*, **4**(2), 191–213

Gottardi F (2009) Estimation statistique et réanalyse des précipitations en montagne. Utilisation d'ébauches par types de temps et assimilation de données d'enneigement. Application aux grands massifs montagneux français. (PhD thesis, Institut National Polytechnique de Grenoble)

Greuell W and Smeets P (2001) Variations with elevation in the surface energy balance on the Pasterze (Austria). *J. Geophys. Res.*, **106**(D23), 31717–31727

Hock R (2005) Glacier melt: a review on processes and their modelling. *Progr. Phys. Geogr.*, **29**(3), 362–391

Hock R and Holmgren B (2005) A distributed surface energy-balance model for complex topography and its application to Storglaciären, Sweden. *J. Glaciol.*, **51**(172), 25–36

Klok EJ and Oerlemans J (2002) Model study of the spatial distribution of the energy and mass balance of Morteratschgletscher, Switzerland. *J. Glaciol.*, **48**(163), 505–518

Le Meur E and Vincent C (2003) A two-dimensional shallow ice-flow model of Glacier de Saint-Sorlin, France. *J. Glaciol.*, **49**(167), 527–538

Le Meur E, Gerbaux M, Schäfer M and Vincent C (2007) Disappearance of an Alpine glacier over the 21st Century simulated from modeling its future surface mass balance. *Earth Planet. Sci. Lett.*, **261**(3–4), 367–374

Lejeune Y (2009) Apports des modèles de neige CROCUS et de sol ISBA à l'étude du bilan glaciologique d'un glacier tropical et du bilan hydrologique de son bassin versant. (PhD thesis, Université Joseph Fourier)

Lejeune Y and 7 others (2007) Melting of snow cover in a tropical mountain environment in Bolivia: processes and modeling. *J. Hydromet.*, **8**(4), 922–937

Martin E and Lejeune Y (1998) Turbulent fluxes above the snow surface. *Ann. Glaciol.*, **26**, 179–183

Martin S (1975) Wind regimes and heat exchange on Glacier de Saint-Sorlin. *J. Glaciol.*, **14**(70), 91–105

Oerlemans J and Knap WH (1998) A 1 year record of global radiation and albedo in the ablation zone of Morteratschgletscher, Switzerland. *J. Glaciol.*, **44**(147), 231–238

Paterson WSB (1994) *The physics of glaciers*, 3rd edn. Elsevier, Oxford

Pedersen CA and Winther J-G (2005) Intercomparison and validation of snow albedo parameterization schemes in climate models. *Climate Dyn.*, **25**(4), 351–362

Seity Y and 7 others (2011) The AROME-France convective-scale operational model. *Mon. Weather Rev.*, **139**(3), 976–991

Sicart JE, Hock R and Six D (2008) Glacier melt, air temperature, and energy balance in different climates: the Bolivian Tropics, the French Alps, and northern Sweden. *J. Geophys. Res.*, **113**(D24), D24113 (doi: 10.1029/2008JD010406)

Sirguey P, Mathieu R and Arnaud Y (2009) Subpixel monitoring of the seasonal snow cover with MODIS at 250 m spatial resolution in the Southern Alps of New Zealand: methodology and accuracy assessment. *Remote Sens. Environ.*, **113**(1), 160–181

Six D, Wagnon P, Sicart JE and Vincent C (2009) Meteorological controls on snow and ice ablation for two contrasting months on Glacier de Saint-Sorlin, France. *Ann. Glaciol.*, **50**(50), 66–72

Slater AG and Clark MP (2006) Snow data assimilation via an ensemble Kalman filter. *J. Hydromet.*, **7**(3), 478–493

Toure AM and 6 others (2011) A case study of using a multilayered thermodynamical snow model for radiance assimilation. *IEEE Trans. Geosci. Remote Sens.*, **49**(8), 2828–2837

Uppala SM and 45 others (2005) The ERA-40 re-analysis. *Q. J. R. Meteorol. Soc.*, **131**(612), 2961–3212

Vincent C (2002) Influence of climate change over the 20th century on four French glacier mass balances. *J. Geophys. Res.*, **107**(D19), 4375 (doi: 10.1029/2001JD000832)

Vincent C, Vallon M, Reynaud L and Le Meur E (2000) Dynamic behaviour analysis of glacier de Saint Sorlin, France, from 40 years of observations, 1957–97. *J. Glaciol.*, **46**(154), 499–506

Wagnon P, Ribstein P, Francou B and Pouyaud B (1999) Annual cycle of energy balance of Zongo Glacier, Cordillera Real, Bolivia. *J. Geophys. Res.*, **104**(D4), 3907–3924

Wagnon P, Lafaysse M, Lejeune Y, Maisincho L, Rojas M and Chazarin JP (2009) Understanding and modeling the physical processes that govern the melting of snow cover in a tropical mountain environment in Ecuador. *J. Geophys. Res.*, **114**(D19), D19113 (doi: 10.1029/2009JD012292)

Warren SG (1982) Optical properties of snow. *Rev. Geophys.*, **20**(1), 67–89

Willemet J-M (2008) *The snow cover model CROCUS. User's guide version 2.4*. Météo France, Centre National de Recherches Météorologiques/Centre d'Etude de la Neige, Saint-Martin-d'Hères
# APPENDIX: VARIATIONAL DATA ASSIMILATION METHOD, DEFINITIONS

This Appendix gives definitions that are essential to understand the methodology (Bouttier and Courtier, 2002).

Let $\mathbf{x} = (x_1, x_2, x_3)$ be the surface grain variables vector and $\mathbf{y} = (\alpha_1, \alpha_2, \alpha_3)$ be the vector composed of the three spectral bands of albedo.

Let $\mathcal{B}$ be the model space (grain variables) and $\mathcal{O}$ the observation space (albedo).

## A.1. Observation operators

We define $\mathcal{H}$, the observation operator, as $\mathbf{y} = \mathcal{H}[\mathbf{x}]$. In the absence of any modeling error, this operator would generate the albedo value, $\mathbf{y}$, if both observations and state vectors were perfect.

The linearized observation operator at $\mathbf{x}_0$, $H$, is defined as

$$ \mathcal{H}[\mathbf{x}] = \mathcal{H}[\mathbf{x}_0] + H(x - x_0) + \circ (\|x - x_0\|_{\mathcal{O}}) \quad (A1) $$

Consequently, $H$ is the true linearized operator of $\mathcal{H}$ if

$$ \forall \mathbf{x}, \mathrm{d}\mathbf{x}, \lim_{\epsilon \to 0} \left| 1 - \frac{\mathcal{H}[\mathbf{x} + \epsilon \mathrm{d}\mathbf{x}] - \mathcal{H}[\mathbf{x}]}{H(\epsilon \mathrm{d}\mathbf{x})} \right| = 0 \quad (A2) $$

The adjoint observation operator, $H^\times$, is defined as

$$ \forall \mathbf{x}, \mathbf{y}, \langle H^\times \mathbf{y}, \mathbf{x} \rangle_B = \langle \mathbf{y}, H \mathbf{x} \rangle_O \quad (A3) $$

## A.2. Errors and error covariances

If $\epsilon$ represents the error vector, then $\mathbf{E} = (\epsilon - \bar{\epsilon})(\epsilon - \bar{\epsilon})^T$ is the error covariance matrix. $\mathbf{x}_b$ is the guess state vector given as a priori information by the model, $\mathbf{x}_t$ is the true state vector, $\mathbf{y}$ is the vector of observations and $\mathbf{x}_a$ is the analysis model state (after assimilation).

$\epsilon_b = \mathbf{x}_b - \mathbf{x}_t$ is the guess error with covariances $\mathbf{B}$. It is the estimation error of the guess state and

$\epsilon_0 = \mathbf{y} - \mathcal{H}[\mathbf{x}_t]$ is the observation error with covariances $\mathbf{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 $\mathcal{H}$ and discretization errors that prevent $\mathbf{x}_t$ from being a true state

$\epsilon_a = \mathbf{x}_a - \mathbf{x}_t$ is the analysis error with covariances $\mathbf{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, $\mathbf{x}_a$, is then defined as the state vector that minimizes this cost function:

$$ \mathbf{x}_a \text{ such that } J(\mathbf{x}_a) = \min_{\mathbf{x} \in \mathcal{B}} J(\mathbf{x}); \quad J(\mathbf{x}) = J_o(\mathbf{x}) + J_b(\mathbf{x}) + J_c(\mathbf{x}) \quad (A4) $$

$J_b(\mathbf{x}) = \|\mathbf{x} - \mathbf{x}_b\|_B^2$ is the guess term, $J_o(\mathbf{x}) = \|\mathbf{y} - \mathcal{H}[\mathbf{x}]\|_{\mathcal{O}}^2$ the observation term and $J_c(\mathbf{x})$ the constraint term that prevents $\mathbf{x}_a$ from escaping from its validity domain.

From Eqn (A4), we can calculate $\nabla J$ and $\nabla^2 J$:

$$ \nabla J(\mathbf{x}) = 2\mathbf{B}^{-1}(\mathbf{x} - \mathbf{x}_b) + 2H^\times \mathbf{R}^{-1}(\mathbf{y} - \mathcal{H}[\mathbf{x}]) + \nabla J_c(\mathbf{x}) \quad (A5) $$

$$ \nabla^2 J(\mathbf{x}) = 2\mathbf{B}^{-1} + 2H^\times \mathbf{R}^{-1} H + \nabla^2 J_c(\mathbf{x}) \quad (A6) $$

To estimate $\mathbf{x}_a$, we use a popular gradient descent algorithm:

$$ \mathbf{x}_a = \lim_{n \to \infty} \mathbf{x}^n, \quad \mathbf{x}^{n+1} = \mathbf{x}^n + (\nabla^2 J(\mathbf{x}^n))^{-1} \nabla J(\mathbf{x}^n). \quad (A7) $$

Note that the analysis error covariance matrix, $\mathbf{A}$, is directly linked to the Hessian of $J$:

$$ \mathbf{A} = \left( \frac{1}{2} \nabla^2 J(\mathbf{x}) \right)^{-1} \quad (A8) $$

MS received 1 August 2011 and accepted in revised form 14 November 2011