WATERRESOURCESRESEARCH, VOL.???, XXXX,DOI:10.1029/, Assimilation of Terrestrial Water Storage from 1 GRACE in a Snow-Dominated Basin 2 1,2 1 3 B. A. Forman , R. H. Reichle , and M. Rodell Index Terms: 1847 (Modeling); 1855 (Remote sensing); 1863 (Snow); 3315 (Data Assimila- 3 tion); GRACE; 4 B. A. Forman, Global Modeling and Assimilation Office, NASA Goddard Space Flight Center, Code 610.1, Greenbelt, MD 20771, USA. ([email protected]) R. H. Reichle, Global Modeling and Assimilation Office, NASA Goddard Space Flight Center, Code 610.1, Greenbelt, MD 20771, USA. M. Rodell, Hydrological Sciences Branch, NASA Goddard Space Flight Center, Code 614.3, Greenbelt, MD 20771, USA. 1Global Modeling and Assimilation Office, NASA Goddard Space Flight Center 2Oak Ridge Associated Universities 3Hydrological Sciences Branch, NASA Goddard Space Flight Center D R A F T October 21, 2011, 1:30pm D R A F T X - 2 FORMANETAL.: GRACE-DAINASNOW-DOMINATEDBASIN Abstract. Terrestrial water storage (TWS) information derived from 5 Gravity Recovery and Climate Experiment (GRACE) measurements is as- 6 similated into a land surface model over the Mackenzie River basin lo- 7 cated in northwest Canada. Assimilation is conducted using an ensemble 8 Kalman smoother (EnKS). Model estimates with and without assimilation 9 are compared against independent observational data sets of snow water 10 equivalent (SWE) and runoff. For SWE, modest improvements in mean 11 difference (MD) and root mean squared difference (RMSD) are achieved as 12 a result of the assimilation. No significant differences in temporal correla- 13 tions of SWE resulted. Runoff statistics of MD remain relatively unchanged 14 while RMSD statistics, in general, are improved in most of the sub-basins. 15 Temporal correlations are degraded within the most upstream sub-basin, 16 but are, in general, improved at the downstream locations, which are more 17 representative of an integrated basin response. GRACE assimilation using 18 an EnKS offers improvements in hydrologic state/flux estimation, though 19 comparisons with observed runoff would be enhanced by the use of river 20 routing and lake storage routines within the prognostic land surface model. 21 Further, GRACE hydrology products would benefit from the inclusion of 22 better constrained models of post-glacial rebound, which significantly af- 23 24 fects estimationGRACE estimates of interannual hydrologic variability in the Mackenzie River basin. 25 D R A F T October 21, 2011, 1:30pm D R A F T FORMANETAL.: GRACE-DAINASNOW-DOMINATEDBASIN X - 3 1. Introduction Snow is an important component of the hydrologic cycle that accounts for a large 26 fraction of the available freshwater resources in many parts of the northern hemisphere 27 [Barnett et al.,2005]. Accurateestimationofsnowmass, orsnowwaterequivalent(SWE), 28 acrossspaceandtimeusingpoint-scale, ground-basedtechniquesisadifficulttask. There- 29 fore, in an effort to better quantify this potential freshwater supply, many researchers have 30 turned to remote sensing estimates derived from space-based instrumentation used in con- 31 junction with land surface models. 32 Despite recent popularity in the utilization of passive microwave and visible spectrum 33 imagery for the purpose of snow pack estimation (e.g., Andreadis and Lettenmaier [2006]; 34 Durand and Margulis [2006];Dong et al.[2007];Su et al.[2008]), satellite-derivedmeasure- 35 menttechniquespossesssignificantlimitations. Passivemicrowaveestimates, forexample, 36 are prone to large errors for snow packs that are either wet, deep (> 1 m), or contain ice 37 and/or depth hoar layers [Clifford, 2010]. Similarly, visible imagery often provides little 38 information outside of the initial accumulation and final ablation periods of the snow 39 season [Clark et al., 2006]. 40 An alternative to passive microwave and visible spectrum-based SWE estimation is the 41 use of gravimetry. Gravimetric techniques focus on the measurement of gravitational 42 anomalies associated with the accumulation (or loss) of mass near the Earth’s surface. 43 In the context of snow, changes in the Earth’s gravitational field are associated with the 44 accumulation of snow during the snow season and the subsequent ablation and runoff of 45 the snow mass during the melt season. Gravimetry is capable of capturing snow mass 46 D R A F T October 21, 2011, 1:30pm D R A F T X - 4 FORMANETAL.: GRACE-DAINASNOW-DOMINATEDBASIN throughouttheaccumulationseason, includingpeakaccumulationwhenSWEinformation 47 is most valuable to water resource managers. Unfortunately, the drawback of space-based 48 gravimetry is its coarse spatial (∼hundreds of km) and temporal (∼monthly) resolution 49 that limits its applicability for smaller domains. When satellite gravimetric measurements 50 are combined with a land surface model as part of a data assimilation (DA) framework, 51 however, there is the potential to effectively downscale gravimetric estimates in time and 52 space while simultaneously providing useful information content when passive microwave 53 and visible spectrum measurements cannot. 54 2. Background One such satellite gravimetry mission is the Gravity Recovery and Climate Experiment 55 (GRACE). GRACE provides approximately monthly estimates of variations in terrestrial 56 water storage (TWS), which includes snow, ice, surface water, soil moisture, and ground- 57 water. The mission is a major step towards understanding regional TWS dynamics [Tang 58 et al., 2010] and offers significant insight into regional- and continental-scale hydrologic 59 processes [Syed et al., 2009; Rodell et al., 2009]. 60 Relatively few studies have been conducted that utilize GRACE measurements within a 61 DA framework. The first study by Zaitchik et al. [2008] assimilated GRACE information 62 into a land surface model of the Mississippi River basin. When compared against in- 63 situ groundwater observations, Zaitchik et al. [2008] found reduced errors and increased 64 temporal correlations as a result of the assimilation. Further, the results suggested the 65 potential to downscale the coarse-scale GRACE measurements via use of a relatively fine- 66 scale land surface model. However, due to the fact that snow contributes little to TWS in 67 D R A F T October 21, 2011, 1:30pm D R A F T FORMANETAL.: GRACE-DAINASNOW-DOMINATEDBASIN X - 5 the Mississippi River basin, there was limited opportunity to study the impact of GRACE 68 data assimilation on snow pack characterization. 69 More recently, Su et al. [2010] studied the impact of GRACE data assimilation on TWS 70 estimates in North America for the express purpose of improved snow pack estimation. 71 They found that GRACE assimilation improved SWE estimation in many of the North 72 American basins where snowfall is a significant contributor to the hydrologic cycle. How- 73 ever, Su et al. [2010] also found that many issues remain to be addressed, including: 1) 74 the cause of model degradation in some high-latitude basins as a result of GRACE assim- 75 ilation, 2) the impact of GRACE observational error on DA results, and 3) the impact of 76 GRACE assimilation on components of TWS other than snow. 77 This study expands on the work by Zaitchik et al. [2008] and Su et al. [2010] via 78 extended examination of GRACE DA performance within a snow-dominated hydrologic 79 basin. Namely, additional verification activities using independent, ground-based data 80 sets are explored, a number of different GRACE products are tested during assimilation, 81 the impact of GRACE measurement error on DA results is investigated, an analysis of DA 82 innovation sequences is included, and a longer period of record is utilized, which allows 83 for a better assessment of inter-annual variability. 84 The following sections introduce the methods used in the assimilation framework (sec- 85 tion 3), highlight the study domain (section 4), discuss the GRACE measurements and 86 forward model used during the assimilation (section 5), highlight the independent data 87 88 sets used for modelverificationvalidation (section 6), present modelassimilation results (section 7), and conclude with summarized findings and implications (section 8). 89 D R A F T October 21, 2011, 1:30pm D R A F T X - 6 FORMANETAL.: GRACE-DAINASNOW-DOMINATEDBASIN 3. Data Assimilation Framework A DA framework is an effective means of merging model estimates with measurements 90 that often yields an improved estimate beyond that of the model or measurements alone 91 [McLaughlin, 2002]. Not only does DA provide a conditioned estimate that accounts 92 for both model and measurement uncertainty, but it offers the potential to effectively 93 downscalethemeasurementsinspaceandtimeviautilizationofthefiner-scaleinformation 94 associated with the prognostic model formulation, its parameters, and its forcing data 95 [Reichle et al., 2001; Zaitchik et al., 2008]. 96 Selection of the most appropriate DA system depends on feasibility, robustness, and 97 computational efficiency. In that regard, we choose to employ an Ensemble Kalman 98 Smoother (EnKS) in part because of its ability to handle non-linear models coupled 99 with its flexible, modular structure [Dunne and Entekhabi, 2006] as well as the ability 100 to leverage Zaitchik et al. [2008] as a precursor study. In general, an EnKS has two ba- 101 sic components: 1) a physically-based, forward model to propagate the model states as 102 an ensemble in order to provide background error covariances, and 2) an update scheme 103 that combines the model states and the satellite measurements in a way that accounts 104 for their respective uncertainties. The work conducted in this current study adapts the 105 EnKS presented in Zaitchik et al. [2008] for a snow-dominated basin thereby contributing 106 to the methodological development of GRACE DA (see section 5.3). The EnKS is first 107 introduced below whereas the assimilated measurements and forward model are discussed 108 in section 5. 109 D R A F T October 21, 2011, 1:30pm D R A F T FORMANETAL.: GRACE-DAINASNOW-DOMINATEDBASIN X - 7 3.1. Ensemble Kalman Smoother The prior (unconditioned) estimate of the model states, xi−xi−, is derived from a prog- 110 t τ nostic land surface model. This is illustrated in the left-hand side (i.e., Step 1) of Figure1. 111 112 The nonlinear model, Ft(·), propagates the posterior (conditioned) model states, xit+−1xiτ+−1, 113 forward in time, t, from t−1 to tone month to the next (i.e., from τ −1 to τ) using an ensemble of N realizations with prescribed model errors wi as 114 t xi− = F xi+ ,wi for i ∈ N. (1) τ t τ−1 t (cid:0) (cid:1) We adopt the convention where bold lowercase symbols denote vectors, bold uppercase 115 symbols denote matrices, non-bold symbols denote scalars, and calligraphic symbols rep- 116 resent operators. Uncertainties in the model are defined by the model error term, wi with 117 t covariance Q . In the ensemble framework, model errors are represented by perturbations 118 t that are applied to model states and forcings (section 5). 119 Next, the prior model states are updated using the observations available for the time 120 period of interest τ ∈ [t ,t ] (where t and t are the beginning and end of the assimilation 121 o f o f 122 periodwindow, i.e., first and last day of the month in this application). This is illustrated in the right-hand side (i.e., Step 2) of Figure 1. The following linear update equation is 123 employed as 124 xi+ = xi− +K y +vi −Hxi− , (2) τ τ τ τ τ (cid:2) (cid:3) where K is the Kalman gain matrix, y is the measurement vector, and H is the pre- 125 τ τ dicted measurement model that linearly maps the model states into measurement space. 126 Randomperturbations, vi, representingmeasurementerrorareaddedtothemeasurement 127 D R A F T October 21, 2011, 1:30pm D R A F T X - 8 FORMANETAL.: GRACE-DAINASNOW-DOMINATEDBASIN vector [Burgers et al., 1998]. The Kalman gain, K , is a weighted average between the 128 τ uncertainty of the prior states and the measurements such that 129 K = P−HT H P−HT +R −1, (3) τ τ τ τ τ τ (cid:0) (cid:1) whereP− isthebackgrounderrorcovariancecomputedfromxi− for i ∈ [1N], andRisthe 130 τ τ measurement error covariance. The analysis increments, x+ −x−, are applied evenly over 131 τ τ each day of the month as illustrated in Step 2 of Figure 1. The update procedure ignores 132 non-Gaussian characteristics and relies only on the first two moments of the distribution. 133 In practice, however, it may only be feasible to accurately compute the first and second 134 moments of the system state [Khare et al., 2008]. Additional details regarding the EnKS 135 update procedure applied in Equation (2) are found in Figure 5 of Zaitchik et al. [2008] 136 as well as in section 5.3 further below. 137 4. Study Domain The study domain used here is the Mackenzie River basin (MRB) located in north- 138 western Canada (Figure 2) and consists of 4 individual sub-basins. Sub-basin delineation 139 was based on topographic control and adhered to the topology of the river network. Each 140 sub-basin was extracted from the original GRACE product in order to produce sub-basin- 141 averaged TWS estimates. The smallest sub-basin is 280,000 km2, which is larger than the 142 minimum area of roughly 150,000 km2 that can be resolved by GRACE at mid-latitudes 143 [Rowlands et al., 2005; Swenson et al., 2006]. Additional details regarding the GRACE 144 measurements and measurement preprocessing activities are found in section 5.1 and sec- 145 tion 5.2, respectively. 146 D R A F T October 21, 2011, 1:30pm D R A F T FORMANETAL.: GRACE-DAINASNOW-DOMINATEDBASIN X - 9 As a whole, MRB is ∼1.8×106 km2 in drainage area (∼1.6×106 km2 for land areas only; 147 see Table 1) with the main branch of the Mackenzie River running from the highlands in 148 the southwestern corner of the domain northward toward the Arctic Ocean. The snow 149 classificationschemeofSturm et al.[2010]suggestsMRBsnowtypeisdominatedbytaiga- 150 type snow with smaller areas of tundra- and alpine-type snow found in the northwest and 151 southern regions, respectively (see Figure 2b). 152 5. Assimilated Measurements and Forward Model 5.1. GRACE Measurements Background Several different GRACE hydrology products were investigated in this study. TWS 153 anomalies from 1) the Space Geodesy Research Group (GRGS) product [Bruinsma et al., 154 2010; Horwath et al., 2011], 2) the Tellus product available from the NASA Jet Propul- 155 sion Laboratory (Tellus) [Wahr et al., 2004; Swenson and Wahr, 2006], and 3) the mass 156 concentration product generated at the NASA Goddard Space Flight Center (MasCon) 157 [Rowlands et al., 2005, 2010]. Each product utilizes the same Level 1 range-rate measure- 158 ments from GRACE, but is processed in a different manner in order to yield mass change 159 estimates in terms of equivalent water thickness. 160 Each product is available as gridded TWS anomalies (i.e., deviations from the temporal 161 ◦ ◦ mean at each location). The GRGS and Tellus products are provided on a ∼ 1 × 1 162 ◦ ◦ grid whereas the MasCon product is provided on a ∼ 4 × 4 grid. Each product was 163 subsequently converted into sub-basin-averaged total TWS values by adding the location- 164 165 specific, temporal meanlong-term average TWS from the land surface model. More infor- mation on GRACE measurement preprocessing is provided in section 5.2 and the land 166 surface model is provided in section 5.3. 167 D R A F T October 21, 2011, 1:30pm D R A F T X - 10 FORMANETAL.: GRACE-DAINASNOW-DOMINATEDBASIN 5.2. GRACE Measurement Preprocessing Conversion of the GRACE TWS anomalies into sub-basin-averaged TWS estimates 168 that are compatible with modeled TWS values begins with generating a single-replicate 169 of the forward model for the period 1 September 2002 to 1 September 2009. No model 170 171 errors are prescribed in this simulation unlike that shown in Equation (1). TemporallyLong- term (i.e., 2002-2009) averaged , sub-basin-averaged estimates of TWS derived from the 172 forward model are subsequently added to the sub-basin-averaged monthly GRACE TWS 173 anomalies, which yields monthly estimates of TWS for each modeled sub-basin that are 174 eventually assimilated using Equation (2). Additional details on the utilization of the 175 GRACE measurements in Equation (2) are found in Zaitchik et al. [2008]. 176 One notable aspect of GRACE preprocessing is the consideration of a secular trend 177 associated with post-glacial rebound (PGR). The Tellus product accounts for PGR using 178 the methods of Paulson et al. [2007]. However, the GRGS and MasCon products do not 179 account for PGR. Therefore, model output from Paulson et al. [2007] is applied here to 180 the GRGS and MasCon products in a similar manner as done for the Tellus product. 181 Preliminary DA results suggest PGR is overestimated by the model of Paulson et al. 182 [2007] in both the Slave and Peace+Athabasca sub-basins, but this cannot be verified 183 as the exact amount of PGR in these regions is unknown. Unfortunately, PGR models 184 are difficult to validate due to a lack of independent data, thus the errors are not well 185 quantified. Therefore, in an effort to better understand the impacts of PGR estimates 186 on GRACE DA performance within the MRB, two different versions of each GRACE 187 product were used in the DA experiments: 1) PGR correction applied using Paulson et 188 al. [2007] and 2) PGR correction not applied (i.e., PGR correction was removed from the 189 D R A F T October 21, 2011, 1:30pm D R A F T