Mass change from GRACE: a simulated comparison of Level-1B analysis techniques

Stuart B. Andrews; Philip Moore; Matt. A. King

doi:10.1093/gji/ggu402

Mass change from GRACE: a simulated comparison of Level-1B analysis techniques

Andrews, Stuart B.; Moore, Philip; King, Matt. A. 2015-01-01 00:00:00 Downloaded from https://academic.oup.com/gji/article/200/1/503/752750 by DeepDyve user on 20 July 2022 Geophysical Journal International Geophys. J. Int. (2015) 200, 503–518 doi: 10.1093/gji/ggu402 GJI Gravity, geodesy and tides Mass change from GRACE: a simulated comparison of Level-1B analysis techniques 1 1 1,2 Stuart B. Andrews, Philip Moore and Matt. A. King School of Civil Engineering and Geosciences, Newcastle University, Newcastle upon Tyne, NE1 7RU, United Kingdom. E-mail: [email protected] School of Land and Food, University of Tasmania, Private Bag 76, Hobart 7001, Australia Accepted 2014 October 13. Received 2014 October 10; in original form 2014 March 31 SUMMARY Spherical harmonic and mascon parameters have both been successfully applied in the recov- ery of time-varying gravity ﬁelds from Gravity Recovery and Climate Experiment (GRACE). However, direct comparison of any mass ﬂux is difﬁcult with solutions generated by differ- ent groups using different codes and algorithms. It is therefore opportune to compare these methodologies, within a common software base, to understand potential limitations associated with each technique. Here we use simulations to recover a known monthly surface mass distri- bution from GRACE KBRR data. The ability of spherical harmonic and mascon parameters to resolve basin-level mass change is quantiﬁed with an assessment of how the noise and errors, inherent in GRACE solutions, are handled. Recovery of a noise and error free GLDAS anomaly revealed no quantiﬁable difference between spherical harmonic and mascon parameters. Ex- pansion of the GLDAS anomaly to degree and order 120 shows that both spherical harmonic and mascon parameters are affected by comparable omission errors. However, the inclusion of realistic KBRR noise and errors in the simulations reveals the advantage of the mascon parameters over spherical harmonics at reducing noise and errors in the higher degree and order harmonics with an rms (cm of EWH) to the GLDAS anomaly of 10.0 for the spherical ◦ ◦ harmonic solution and 8.8 (8.6) for the 4 (2 ) mascon solutions. The introduction of a con- straint matrix in the mascon solution based on parameters that share geophysical similarities is shown to further reduce the signal lost at all degrees. The recovery of a simulated Antarctic mass loss signal shows that the mascon methodology is superior to spherical harmonics for this region with an rms (cm of EWH) of 8.7 for the 2 mascon solution compared to 10.0 for the spherical harmonic solution. Investigating the noise and errors for a month when the satellites were in resonance revealed both the spherical harmonic and mascon methodologies are able to recover the GLDAS and Antarctic mass loss signal with either a comparable (spherical harmonic) or improved (mascon) rms compared to non-resonance periods. Key words: Inverse theory; Satellite geodesy; Gravity anomalies and Earth structure; Time variable gravity; Global change from geodesy; Antarctica. is measured using a K/Ka band instrument which has an accu- 1 INTRODUCTION racy of a few microns (Dunn et al. 2003). The ﬁrst time-derivative The Gravity Recovery and Climate Experiment (GRACE) twin provides a K/Ka band range rate (KBRR) measurement, which is satellite mission (Tapley et al. 2004b) was launched in 2002 with generally used during GRACE data processing and has a precision –1 the aim of mapping the Earth’s time-varying gravity (TVG) ﬁeld of 0.1 µms (Luthcke et al. 2006b). The GRACE satellites are every ∼30 days (Wahr et al. 2004) at spatial scales of hundreds of equipped with Global Positioning System (GPS) receivers to posi- kilometres. GRACE consists of two near-polar satellites separated tion the satellites and on-board accelerometers to measure the non- along-track by ∼220 km. This separation results in the two satel- conservative forces that affect the satellites allowing, shortly after lites experiencing small orbital perturbations, caused by local mass launch, the detection of changes in mass corresponding to 1.4 cm anomalies, at different times. It is through measuring the chang- of equivalent water height (EWH) over a radius of 750 km every 30 ing distance between the satellites, caused by these perturbations, days (Wahr et al. 2004). Analysis improvements have since driven that GRACE allows us to monitor spatial and temporal geophysical this uncertainty/radius downward. The data from these instruments mass variations in the gravity ﬁeld. The range between the satellites are available to the scientiﬁc community as Level-1B data, namely The Authors 2014. Published by Oxford University Press on behalf of The Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. Downloaded from https://academic.oup.com/gji/article/200/1/503/752750 by DeepDyve user on 20 July 2022 504 S. B. Andrews, P. Moore and M. A. King the raw satellite products after a low-pass ﬁlter and time correction matrices allows users to constrain the higher degrees of the TVG (Chen et al. 2008). ﬁelds to the mean ﬁeld (e.g. Bruinsma et al. 2010;Save et al. The ability to accurately recover the Earth’s temporal mass 2012) but these are not typically provided as part of the Level-2 anomalies every 30 days to a few hundred kilometres allows the solutions and are only available to users generating solutions from monitoring of global geophysical processes between successive 30 Level-1B data. day solutions and provides a means to assess and improve geo- Mass concentration (mascon) parameters, estimated using only physical models. GRACE data has been successfully used for es- GRACE Level-1B KBRR data, have also been applied to the mod- timating ice mass balance trends (e.g. King et al. 2012; Shepherd elling of TVG, primarily at the NASA Goddard Space Flight Centre et al. 2012; Luthcke et al. 2013), hydrological change (e.g. Han (GSFC; Rowlands et al. 2005, 2010; Luthcke et al. 2006b, 2008; et al. 2009; Awange et al. 2011; Ramillien et al. 2011), ocean mass Sabaka et al. 2010). The estimation of mascon parameters, ﬁrst variability (e.g. Chambers et al. 2010; Chambers & Bonin 2012), applied by Muller & Sjogren (1968) to model lunar gravity vari- Glacial Isostatic Adjustment (GIA; e.g. Riva et al. 2009) and study- ations, has allowed the Earth’s TVG to be modelled at a spatial ing Earthquake displacements (Han et al. 2008). Without the spatial and temporal resolution unachievable using conventional spherical and temporal resolution provided by GRACE, the separation of the harmonics due to the embedded smoothing resulting in less sig- signals from many of these geophysical phenomena would not have nal loss (Rowlands et al. 2010). The mascon methodology beneﬁts been possible (Swenson & Wahr 2002). from constraint equations applied during the least squares inversion Since the launch of GRACE several processing methodologies which removes the requirement of post-processing, making possible have been applied to the GRACE Level-1B data with the aim of es- solutions with a 10 day temporal resolution and a spatial sampling ◦ ◦ timating the temporal variations of mass anomalies while reducing of 4 or 2 (Rowlands et al. 2010) with the most recent published the impact that noise and errors have on the solutions. The pres- mascon solutions using a 1 spatial sampling (Luthcke et al. 2013). ence of north–south ‘stripes’ in the solution results from the limited The mascon constraint matrix is anisotropic in nature as it can be longitudinal sampling of GRACE, which reduced observability of applied to constrain distinct geographical regions while also taking the sectorial and near-sectorial coefﬁcients (Klees et al. 2008), and into account the full noise covariance matrix (Luthcke et al. 2013). result from the orbital conﬁguration of GRACE; the KBRR data The use of mascon parameters has allowed basin scale interannual are only collected along the line of the orbit. Stripes also suggest a variability to be extracted from the GRACE time-series (Rowlands high degree of correlation in the GRACE errors (Swenson & Wahr et al. 2010) due to the increased spatial sampling and temporal res- 2006). If not correctly dealt with, measurement and processing er- olution of the solutions. Mascon TVG estimates, unlike spherical rors and unmodelled mass variations in the GRACE solution could harmonic coefﬁcients, can be either local or global in extent as the be interpreted as real mass change (Wahr et al. 2006). Random KBRR data associated with the mascon is dominated by mass ﬂux and systematic errors, present in the instrument data and back- below the satellite (Luthcke et al. 2008). Mascon parameters can be ground models, can dominate the solution at short wavelengths and estimated from short arcs of KBRR data through the estimation of increase rapidly with the spherical harmonic degree (Klees et al. short arc baseline parameters (Rowlands et al. 2002). However this 2008) while errors in geophysical models can cause poor estimates requires GRACE accelerometer data to deal with non-conservative of certain spherical harmonic coefﬁcients and alias into the solu- forces and the satellite orbits to have been previously well deter- tions (Seo et al. 2008). In addition, the ground coverage of the mined using GNV1B (GPS Navigation Data Level-1B) positioning satellites also affects the quality of the solutions. Periods where the data. Estimating mascon parameters from Level-1B data is differ- GRACE satellites are in near 3, 4 and 7 day repeat orbits modify ent from the mascons used in Jacob et al. (2012) and Schrama the noise characteristics of the GRACE solutions (Save et al. 2012). et al. (2014), where mascon parameters are calculated from Level-2 Noise and errors in the GRACE solutions need to be appropriately spherical harmonics data. treated and mitigated before any meaningful mass anomalies can be Both spherical harmonic and mascon methodologies are applied estimated. to the recovery of TVG, but most users of GRACE data are heav- ily reliant on pre-processed Level-2 spherical harmonic products. The most common methodology used to estimate the temporal variations of mass anomalies is the recovery of spherical harmonic It is therefore important that an assessment of these methodologies (Stokes’) coefﬁcients (Tapley et al. 2004b;Wahr et al. 2004). The is undertaken to provide users with an understanding of their as- GRACE data processing centres, namely Centre for Space Re- sociated limitations and to assess the ability of each technique to search (CSR), Jet Propulsion Laboratory (JPL) and GeoForschungs resolve basin-level mass changes at a variety of spatial scales. The Zentrum (GFZ), each provide monthly estimates of the global TVG need for this assessment is further compounded as these processing ﬁeld in the form of spherical harmonic ﬁelds to degree and order 60– methodologies are applied by different groups using different codes 90. These Level-2 products are normally computed directly from and algorithms, making it hard to directly compare any subsequent Level-1B KBRR and GPS data, but Level-2 solutions have been suc- mass ﬂux analysis. Pritchard et al. (2011) highlighted that differ- cessfully computed using GRACE KBRR data only (Luthcke et al. ences within Level-1B processing applied at different groups is the 2006a). Due to the complex nature and computational requirement main source of the discrepancies in TVG solutions, with solutions of Level-1B data processing, Level-2 products are widely used by generated using the same fundamental processing software agreeing researchers who undertake studies of mass ﬂux analysis. A num- within estimated errors (see also Shepherd et al. 2012). ber of post-processing strategies have been proposed to remove While there are a limited number of studies which have compared the noise and errors in the solutions including the use of Gaus- different approaches within the same software base (e.g. Rowlands sian smoothing (e.g. Wahr et al. 1998, 2004; Tapley et al. 2004a), et al. 2010; Awange et al. 2011), no study, as far as we are aware, ﬁltering using empirical orthogonal functions (EOF; e.g. Wouters has been undertaken that simulates the recovery of a known sig- &Schrama 2007), destriping (e.g. Swenson & Wahr 2006)and nal using both mascon and spherical harmonic solutions within forward modelling (e.g. Wouters et al. 2008; Schrama & Wouters the same software. Accordingly, we undertake a comparison of 2011). Gaussian smoothing is probably the most common process- solutions generated through the estimation of mascon parameters ing strategy (Rowlands et al. 2010). Access to solution normal and spherical harmonic coefﬁcients. Simulations will provide an Downloaded from https://academic.oup.com/gji/article/200/1/503/752750 by DeepDyve user on 20 July 2022 Mass change from GRACE: simulated comparison 505 accuracy assessment and quantify the capability of each technique to resolve basin-level mass changes at a variety of spatial scales while understanding how the methodologies handle the noise and errors inherent at higher degree and orders. Actual noise and errors, obtained during GRACE Level-1B data processing, will be used to create simulated stripes for use during the simulations. While we only simulate the recovery of a single signal, the assessment of the methodologies aims to provide users with an understanding of their associated limitations, as correctly recovering individual monthly solutions is imperative when using these solutions for de- riving longer term trends. We will use simulations to compare a mascon and spherical har- monic solution using a standard Gaussian smoothing ﬁlter for the post-processing of the spherical harmonic solution and a spatial con- straint matrix applied to the mascon solution. The aim of this paper is to investigate the effect of the fundamental processing methodol- ogy applied to the GRACE Level-1B data on the recovered monthly time variable gravity ﬁeld. It is worth noting that different post- processing strategies applied to the spherical harmonic solutions will generate different mass ﬂux estimates, however only the dif- ference in the fundamental Level-1B processing strategy will be considered in this work. The outline of the paper is as follows. In Section 2 we describe the Figure 1. Distribution of 2 equal area mascons around Antarctica. Antarc- formation of mascon parameters before discussing the processing tic drainage basins (Zwally et al. 2012) are also shown (red). methodology in Section 3. The results of the simulations will be presented in Section 4. Throughout this paper special attention will Using eq. (1) each mascon in the grid can be expanded into be paid to the formation and behaviour of the mascon parameters to spherical harmonics as aid understanding of the methodology. Spherical harmonic solutions l l max (m) (m) are widely dealt with in the literature and their estimation is not m ˜ ˜ δ (φ, λ) = R ρ C cos mλ + S sin mλ e w lm lm explicitly described here. l=0 m=0 ( ) × P cos φ , (2) lm (m) (m) ˜ ˜ where C and S are the dimensionless normalized surface 2 FORMATION OF THE MASCON lm lm harmonics and P normalized associated Legendre polynomials. PARAMETERS lm The sum over degree l and order m in eq. (2) is truncated tol = 60. max The mascon parameters are based on the GSFC methodology which We utilize eq. (2) to modify eqs (11) and (12) of Wahr et al. (1998)to (m) involves calculating a scaling factor for a set of ‘delta’ spherical obtain (delta) spherical harmonic gravitational coefﬁcients, C lm harmonic coefﬁcients that represent a uniform layer of surface mass (m) and S , that describe the required mass spread uniformly over a lm over an area which can be added to the mean background ﬁeld and global set of mascons to give used to represent the mass ﬂux at a certain time. Rowlands et al. 3p 1 + k w l (m) (m) (2010) showed that the only difference between these differential C = H(t) C lm lm 4π R ρ 2l + 1 coefﬁcients is the shape, size and location of the area. Each separate e av area is a mascon as deﬁned by estimating a scale factor H(t) at time 3p 1 + k (m) w l (m) S = H(t) S . (3)) t for a uniform mass of one unit of EWH over the mascon area, lm lm 4π R ρ 2l + 1 e av namely 5517 kg ⎧ ⎫ In eq. (3), ρ is the average density of the Earth ( )and k av 3 l ⎨ 1, over mascon⎬ the load Love number of degree l (Farrell 1972) to account for the δ (φ, λ) = R ρ , (1) deformation of the load on an elastic Earth. The mascon parameter 0, elsewhere ⎩ ⎭ (H) is a scale factor of the spherical harmonic coefﬁcients used to describe that mascon (Rowlands et al. 2010). During the orbital where the mascon is centred on geodetic latitude φ, longitude λ; integration we use the spherical harmonic coefﬁcients from eq. (3) R is the radius of the Earth and ρ the density of fresh water e w –3 to calculate H by least-squares utilizing the partial derivatives of the (1000 kg m ). KBRR observation data (cf. Luthcke et al. 2013), that is In this paper, we use global equal area masons, where each mas- ◦ ◦ con is rectangular in latitude and longitude spanning 2 or 4 in l max l ∂KBRR ∂KBRR ∂C I I lm latitude with the longitude span selected to give approximately the A = = IJ ∂ H ∂C ∂ H ◦ ◦ ◦ ◦ J J lm same area as a 2 by 2 or 4 by 4 mascon located at the equator. l=1 m=0 ◦ ◦ ◦ ◦ While not strictly 2 by 2 or 4 by 4 , for simplicity we refer to ∂KBRR ∂S lm them as such hereafter. The only exception is at the poles where + , (4) ∂S ∂ H ◦ J lm polar caps of 2 in latitude are used. Due to the convergence of longitude at the poles, the number of mascons in a given latitude where A are elements of the design matrix, ∂KBRR /∂ H is IJ I J band will reduce poleward as seen in Fig. 1 for a grid of 2 mascons the partial derivative of the KBRR observation i with respect to the J J over Antarctica. mascon parameter j, ∂KBRR /∂C and ∂KBRR /∂S are I I lm lm Downloaded from https://academic.oup.com/gji/article/200/1/503/752750 by DeepDyve user on 20 July 2022 506 S. B. Andrews, P. Moore and M. A. King the partial derivatives of the KBRR observations with respect to the and S = 0.0001 were tested but were found to over/under constrain J J delta gravitational coefﬁcients while ∂C /∂ H and ∂S /∂ H the simulated solutions. However, the choice of S is dependent on J J lm lm are the partial derivatives of the delta gravitational coefﬁcients with the values chosen for other parameters in eqs (8) and (9) as the respect to the mascon parameter j.Let xˆ be the least squares solution weight matrix V is affected by the values assigned to S, D and T in to the KBRR linearized observation equations Ax − b − v = 0 eqs (8) and (9). The same scale factor is used for all constraints in with observation weight matrix W and Gaussian error v . The solu- the V matrix. tion vector xˆ includes GRACE orbital parameters and either mascon The need for temporal smoothing increases as the number of days parameters, H, or corrections H to mascons parameters H at the included in the mascon solution increases. For consistency with the th k stage of an iterative procedure (cf. Sabaka et al. 2010). The vector spherical harmonic approach, no temporal constraint is applied to b is formed from the KBRR calculated minus observed residuals. the mascon monthly solutions presented here, leading to a modiﬁed Analogously, we incorporate spatial and temporal constraints on version of the constraint matrix, namely the mascons through pseudo-observations Cx − q − v = 0 with 2 d IJ 1− V = S × e . (9) weight matrix V,where q is a constant vector incorporating both IJ non-zero constraints and the component −CH if the procedure is In our solutions, constraints are applied between all mascon pairs iterative. Then, xˆ is obtained from of the same type (land or ocean), but no constraint is applied between T T T T a pair consisting of a land and an ocean mascon. The constraint (A WA + C VC)xˆ = A Wb + C Vq. (5) matrix can be easily modiﬁed to allow different constraints between The matrix, C, incorporates the mascon spatial and/or temporal pairs of mascon, as we will show in Sections 4.6 and 4.7. constraints. For mascons, J , J , with the same geophysical property 1 2 the pseudo-observations can be written as 3 SIMULATED DATA PROCESSING H − H = 0, (6) J J 1 2 STRATEGY where J = J ; J , J ∈ A ,the kth mascon region. Additional 1 2 1 2 k During simulations we aim to recover a known surface mass dis- pseudo-observations are applied to ensure total mass is conserved tribution from GRACE KBRR data using Newcastle University’s (equivalent to a degree 0 term) and that the degree 1 harmonics orbit determination software Faust (Moore et al. 1999). We compute for surface mass are zero; GRACE is insensitive to degree 1 har- the effect of the surface mass distribution through a spherical har- monics as the satellites orbit around the instantaneous centre of monic expansion to a speciﬁed degree and order, with the resultant mass. The inclusion of background models during the orbital in- gravitational ﬁeld added to the mean GIF48 ﬁeld (Ries et al. 2011), tegration could result in non-zero degree 1 terms and mass not from which simulated noise and error free KBRR tracking data are being conserved in the solution vector without the use of pseudo- computed. The simulations then attempt to recover the known sig- observations; individual models do not conserve mass and may nal from the tracking data up to degree and order 60. The normal give non-zero degree 1 coefﬁcients (Wahr et al. 1998). The degree equations from each arc, here taken as 60 min, were combined to 0 and 1 pseudo-observations are prescribed a high weight in V to create monthly solutions for the two methodologies. We only use ensure the solution essentially satisﬁes these conditions. The four KBRR data in our solutions as Luthcke et al. (2006a) showed that constraint equations for degrees l = 0, 1 and order m = 0, 1 are of the gravity ﬁeld can be successfully recovered from KBRR data us- the form ing spherical harmonic parameters while estimating the Rowlands (J) (atmos) (ocean) H C =−C − C , (7) J et al. (2002) short arc baseline parameters. Mascon parameters are l,m l,m l,m usually estimated using only KBRR data. The use of KBRR data (atmos) (ocean) only removes the computational requirement of processing the GPS where the summation is over all mascons and C and C l,m l,m data which can weaken and corrupt the KBRR measurement if not are the harmonics from the GRACE dealiasing product for the ocean handled correctly (Luthcke et al. 2006a). and atmosphere. By stipulating eq. (7), the total contribution due Within our simulations, the only difference between the two pro- to the mascons for the degree 0 and 1 harmonics will offset the a cessing strategies is the gravity ﬁeld parameters being estimated priori mass of the ocean and atmosphere if used in the GRACE (spherical harmonic or mascon). Thus, the underlying processing is computations, giving zero values for these harmonics at all times. the same for the two methodologies to the extent that any differences The weight matrix V utilizes constraints between the mascon in the solution will be directly due to the parametrization (Rowlands parameters, et al. 2010) and how noise and errors propagate into the solutions. d t IJ IJ 2− − Mathematically the two methodologies are equivalent (Sabaka et al. D T V = S × e , (8) IJ 2010). where d is the distance in kilometres between the mascons, D The choice of arc length was based on experimentation in deriving IJ the correlation distance, t the time difference in days between gravity ﬁeld solutions using KBRR (and accelerometer) data with IJ the solutions, T the correlation time and S a scale factor. The scale the satellite state vectors estimated simultaneously using the short factor is required to ensure that the constraints are in proportion baseline parameters. No GPS/GNV1B data or empirical parameters to the normal matrices. Setting the constraint too small will not are required to stabilize the solutions. Analysis of 4 months of substantially reduce the noise and error in the solution, while adding solutions revealed that the average ﬁt of the KBRR residuals was −1 −1 a large constraint will oversmooth the signal and remove both noise 0.32 µms using 60 min arcs and 0.31 and 0.41 µms for 30 and error and real signal (Klosko et al. 2009). The value used for min and 90 min arcs, respectively. The 60 min arcs were favoured S, to be discussed later, depends on the units used for the normal over the 30 min arcs due to the reduction in storage requirement. matrices and constraint matrix. In Section 4.5, we found S = 0.001, 24 hr arcs were used to calculate the GRACE accelerometer bias the value used by Lemoine et al. (2007) and Rowlands et al. (2010), values before recovery of the gravity ﬁeld parameters. Alongside to be the optimal scale factor for the simulations. Values of S = 0.01 the accelerometer biases, these runs estimate the full state vector Downloaded from https://academic.oup.com/gji/article/200/1/503/752750 by DeepDyve user on 20 July 2022 Mass change from GRACE: simulated comparison 507 from KBRR, GNV1B and accelerometer data while using empirical ogy signal will reduce the leakage in the GRACE signal as seen parameters to absorb unmodelled KBRR measurement variations in Sabaka et al. (2010). However, as this assumes that sufﬁciently for quality control purposes. accurate models are available, forward modelling of hydrology will The results presented in Section 4.2 for a noise and error free not be considered in this paper. The simulation was repeated with simulation show that we are able to create KBRR tracking data the source area described to spherical harmonic degree and order that describes a simulated mass. Noise and errors are included in 120 (not shown). The results reveal that the additional degree and the simulation in Section 4.5, through creation of simulated stripes, order 61–120 has minimal effect on KBRR residuals conﬁrming the allowing us to examine how the solutions differ at higher degree result seen in Rowlands et al. (2010). and orders using the different methodologies. 4.2 Simulated recovery of the GLDAS anomaly 4 SIMULATIONS The next simulation was the recovery of a known gravity signal ◦ ◦ through estimation of spherical harmonic and 2 and 4 mascon pa- 4.1 Effect of a simulated mass on observed KBRR rameters to assess the performance of the different methodologies residuals in an ideal case. The introduced signal was the June 2006 Global Land Data Assimilation System (GLDAS) Version 1 (Rodell et al. To understand the performance of the two methodologies a number 2004) anomaly to the 2003–2010 mean derived from monthly av- of experiments were undertaken. The ﬁrst experiment examined the eraged data with a 1 spatial resolution. Since hydrology is the effect of a simulated mass on the KBRR data. The simulated mass largest seasonal contribution to the GRACE signal we are sim- was created by loading a 4 block at the equator, described from ulating the recovery of the main spatial and temporal variation, spherical harmonic degree and order 2–60, with 30 cm of EWH. The with similar characteristics to that contained within real GRACE spherical harmonic expansion was added to the mean background data. We scaled the GLDAS anomaly upwards so the degree vari- ﬁeld to create error free KBRR tracking data. The effect of the mass ances contained a similar power as the June 2006 CSR release on the KBRR residuals is plotted in Fig. 2, with t = 0 at the centre of 5 Level-2 ﬁeld (Bettadpur 2012), as not all geophysical signal is the block. The simulation was repeated with the source mass loaded present in the GLDAS anomaly. The GLDAS anomaly was derived with 10 cm of EWH. The boundary of the 4 block is plotted giving to degree and order 60 and 120 (Fig. 3) with tracking data created a spatial domain ∼440 km in diameter. from the different harmonic expansions. As the degree variances Fig. 2 shows that both simulated masses modify the KBRR resid- are identical up to degree and order 60, recovery of the GLDAS uals inside and outside the source area as found by Rowlands et al. anomaly described to degree and order 120 allows us to observe how (2010). The largest KBRR values are observed at the boundary parameter estimation is affected by omission errors. Gravity ﬁeld re- of the source area with a change from positive to negative as the covery was undertaken, for both GLDAS anomaly models, through GRACE satellite pair move over the centre. Either side of the source area is an extensive tail revealing that GRACE is sensitive to mass over a large area. Spatially, Fig. 2 covers an area of ∼2500 km. The amplitude and extent of the tail is a function of the size of the mass anomaly. Seo et al. (2006) found that the distribution of leakage is concentrated in regions with high hydrological variability. As the hydrology component is the largest seasonal signal in the GRACE data, Fig. 2 shows mass signal leakage from the source area into surrounding areas. These results, and the ﬁndings of Seo et al. (2006), lead to the conclusion that forward modelling of the hydrol- Figure 3. GLDAS mass anomaly for June 2006 as differenced from the Figure 2. The residuals resulting from the 4 block loaded with 30 cm of 2003 to 2010 mean to degree and order 60 (a); degree and order 120 (b). EWH (green) and 10 cm of EWH (blue). The boundary of the simulated Note that there is no GLDAS signal in Antarctica or Greenland or over the mass (red) is plotted for reference. oceans. Units of cm of EWH. Downloaded from https://academic.oup.com/gji/article/200/1/503/752750 by DeepDyve user on 20 July 2022 508 S. B. Andrews, P. Moore and M. A. King Figure 4. Rms of the differences between the calculated gravity ﬁeld solu- Figure 5. Degree variances of the spherical harmonic and mascon solutions. tions and the noise and error free GLDAS anomaly in cm of EWH. GLDAS The degree variances of the noise and error free input GLDAS anomaly are described to degree and order 60 (a); and degree and order 120 (b). The included for comparative purposes. distance refers to the Gaussian smoothing radius (spherical harmonic coef- ﬁcient solution) or the correlation distance (mascon solutions). The distance GLDAS anomaly is 0.05 cm of EWH, which is the same as the used in the mascons solutions (triangles) is not directly comparable to that spherical harmonic recovery with 50 km Gaussian smoothing. In- used in the harmonic solution (crosses) and do not equate to the same creasing the correlation distance in the constraint matrix increases smoothing radius. the effect of the constraint matrix and over constrains the mascon solution, reducing high frequency signal. spherical harmonic estimation to degree and order 60 (3717 pa- Fig. 4(b) shows the rms of the recovered anomaly to the input ◦ ◦ rameters) while the 4 (2 ) equal-area mascon solutions involved GLDAS data, described to degree and order 120. This allows us to estimation of 2564 (10 292) mascon parameters. The monthly so- compare the effect of omission errors on the different parameters. lution was collated from the individual 60 min arcs. No noise or The pattern is similar to the ideal case (GLDAS to degree and order errors were added to the solution at this stage. 60), but with an increase in all rms values. Methodologies that are Fig. 4(a) shows the rms in cm of EWH of the recovered anomaly to essentially truncated to degree and order 60 are, of course, unable the input GLDAS data, described to degree and order 60. This is the to fully recover an input signal described to degree and order 120 ideal case. The rms is provided for various smoothing distances. For and hence unable to capture the additional short wavelength (high spherical harmonic coefﬁcients the distance refers to the Gaussian degree) information. smoothing radius, while for the mascon solutions the distance is D Fig. 5 shows the degree variances of the spherical harmonic and of eq. (9). As such, the two distances are not directly comparable mascon solutions to degree 60. The degree variance of the input as they do not equate to the same smoothing radius. The distance GLDAS anomaly is included for reference. All three solutions match 0 km indicates that no Gaussian smoothing has been applied or that the degree variance of the input GLDAS anomaly to near degree the mascon solutions are unconstrained. 40 from where the 4 mascon solution loses power. This reduction Fig. 4(a) reveals that, without Gaussian smoothing, the spherical is the result of the under sampling of the input signal by the 4 harmonic solution is able to exactly recover the reference degree mascon solution. The spherical harmonic coefﬁcient and 2 mascon and order 60 GLDAS anomaly using perfect data. Adding Gaussian solutions are able to replicate the degree variances of the GLDAS smoothing dampens the amplitude of the recovered signal. For the anomaly. The slight reduction in power at degree ∼55–60 in the 4 mascon solution, the smallest rms is obtained when a correlation 2 mascon solution is the result of the 50 km correlation distance distance of 150 km is used in the constraint matrix, but the solution used in the constraint matrix to stabilize the inversion. The ability is not able to fully recover the input GLDAS anomaly. This is to be of the spherical harmonic and 2 mascon solutions to replicate expected as the GLDAS anomaly was expanded to spherical har- the degree variances of the GLDAS anomaly up to degree 60 in monics degree and order 60 (3717 gravity ﬁeld parameters) while Fig. 5 shows that we are able to reproduce a known input signal the 4 mascon solution involves estimation of 2564 mascon param- from the simulated KBRR tracking data, validating the processing eters. The recovery favours the spherical harmonic solution over the methodology. For the remainder of the paper, all processing uses 4 mascon solution with the latter unable to recover the gravity ﬁeld. the GLDAS anomaly described to degree and order 120 and hence The use of the constraint matrix in the mascon solution does enable allows consideration of omission errors. a small improvement in the ability to recover this short wavelength signal. The 2 mascon solution involves estimation of 10 292 mascon 4.3 Mascon centre-to-centre distance parameters and consequently nearly a factor 4 more parameters than the spherical harmonic solution. If no constraint matrix is As part of the evaluation of the force model the contribution of applied, the solution does not invert as it is rank deﬁcient by a every mascon to the accelerations of the satellites is calculated. large margin (Rowlands et al. 2010). Adding a constraint matrix Thus, as for spherical harmonics, mascon solutions presented so with ≥50 km correlation distance stabilizes the solution and allows far are global in extent with the drawback of added computational recovery of the GLDAS anomaly. The resulting rms to the input expense, at least in the case of the 2 mascon solutions. Fig. 2 Downloaded from https://academic.oup.com/gji/article/200/1/503/752750 by DeepDyve user on 20 July 2022 Mass change from GRACE: simulated comparison 509 Figure 7. Degree variances of the 2 mascon solutions for a range of Figure 6. Degree variances of the 4 mascon solutions for a range of centre-to-centre cut-off distances. The degree variances of the input GLDAS centre-to-centre cut-off distances. The degree variances of the input GLDAS anomaly (solid green) and original solution (solid black) are included for anomaly (solid green) and original solution (solid black) are included for reference. All solutions here are noise and error free. reference. All solutions here are noise and error free. revealed that a mass anomaly leaks outside of the source area, and by extension neighbouring mass anomalies will leak into the source area. Fig. 2 also showed that this leakage is not global, but spatially limited. The limited spatial extent of a mass anomaly on gravitation experienced at the satellite altitude was noted by Velicogna & Wahr (2006). This suggests that we can achieve an adequate approximation of the total contribution of all mascons to the satellite accelerations by considering a reduced subset of mascons. This will allow estimation of a local mascon solution and reduce the computational requirement in a global solution. To investigate this we repeated the simulated recovery of the GLDAS anomaly in Section 4.2, but this time we estimate only the contributions to the satellite accelerations of a subset of mascons. The number of mascons included is based on the distance between the centre of the mascon under consideration to the mascon directly below the leading GRACE satellite. We will call this the centre-to- centre distance and the mascon below the satellite the nadir mascon. Any mascon pair where the centre-to-centre distance is greater than Figure 8. The contribution of mascons to the gravitational acceleration of a prescribed value will be excluded from eq. (4) and its effect a satellite over the Amazon as a function of the centre-to-centre distance. set to zero. Thus, for each arc only a subset of mascons will be Acceleration is given as a percentage of the acceleration of the nadir mascon. included in the normal equations, reducing the processing time and storage requirement. All the individual arcs will be summed to increasing the centre-to-centre distance to 6000 km improves the ﬁt calculate a monthly solution, as previously. The solutions are global of the degree variances with the original mascon solution, at which ◦ ◦ in extent as each month all mascon parameters are estimated in the point the degree variances for both the 2 and 4 mascon solutions ◦ ◦ monthly solution. Mascon 4 and 2 solutions were generated using closely resemble the global solution. A centre-to-centre distance various centre-to-centre distances. The resulting degree variances of 4000 km or 5000 km could be used, but the solutions have less are plotted in Figs 6 and 7, respectively and compared to those from power at degree 2 and 3 harmonics (longer wavelength). the global mascon solutions and input GLDAS anomaly. The 0 km The contribution to the acceleration experienced by the GRACE centre-to-centre distance solution refers to the contribution of the satellites at distances from the nadir mascon are plotted in Fig. 8,for nadir mascon only to the acceleration experienced by the GRACE a nadir mascon located in the Amazon (similar results are observed satellites. when the nadir mascon is located elsewhere). The accelerations are Fig. 6 reveals that for a centre-to-centre cut-off distance of given as a percentage of the acceleration of the nadir mascon. Fig. 8 1000 km or less the degree variances of the original mascon so- reveals that, as expected, mascons close to the nadir mascon have a lution are not well recovered as the leakage seen in Fig. 2 has not larger contribution to the accelerations experienced by the satellites been fully accounted for. All such solutions have reduced power at than distant mascons. As the centre-to-centre distance increases, higher degrees, which dampens short wavelength information in the the contribution to the satellite acceleration reduces and approaches gravity ﬁelds. The 0 km centre-to-centre distance solution in Fig. 7 zero. Mascons within a centre-to-centre distance of up to 2000 km (2 mascon) is anomalous with more power than expected, especially have the greatest contribution but mascons 2000 km from the nadir at lower degrees due to the leakage of the spherical harmonics used mascon still contribute ∼10 per cent of the acceleration of the nadir to deﬁne the mascon parameters in eq. (3). In both Figs 6 and 7, mascon. Mascons with centre-to-centre distances up to 4000 km are Downloaded from https://academic.oup.com/gji/article/200/1/503/752750 by DeepDyve user on 20 July 2022 510 S. B. Andrews, P. Moore and M. A. King still seen to inﬂuence the accelerations experienced by the satellites. through to 6000 km at which the difference in the rms is only 0.06 cm Above 4000 km the contributions to the accelerations reduce, ex- of EWH to the original (all mascon) solution. A similar result is plaining the relatively small improvements at higher degrees seen observed for more limited testing of the 2 solution (not shown). in Figs 6 and 7. A centre-to-centre cut-off distance of 6000 km will Iterations improve the solution but are not a substitute for using be used in future processing as the difference between that and a an appropriate centre-to-centre cut-off distance and our simulations full global solution is minimal. show iterations were unable to overcome modelling deﬁciencies associated with an insufﬁcient distribution of mascons when using perfect data. Iterations will therefore not be used in the remainder of this paper. Note that this iterative procedure is not intended as a study 4.4 Iteration of the mascons of the use of iterating a correctly parametrized solution. Applying In Section 4.3 a good approximation to a true mass anomaly (sim- iterations to real GRACE observations, as discussed in Luthcke ulated with GLDAS) was obtained by including the contribution of et al. (2013), would be expected to yield improvements due to the all mascons within a centre-to-centre distance of 6000 km. How- separation of signal and noise through enhanced quality control of ever, it is also worth investigating if a similar approximation can the residuals. Sabaka et al. (2010) found that forward modelling be achieved by using a smaller centre-to-centre cut-off distance but improved signal recovery while resulting in KBRR residuals that also iterating the solution. Iteration was shown by Luthcke et al. conformed more closely to zero-mean Gaussian distributions. (2013) to substantially increase signal recovery. Here, we computed an iterated 4 mascon solution for the centre-to-centre distances of 0, 1000 and 2000 km and computed the rms of the differences to 4.5 Simulated recovery of the GLDAS anomaly with the known GLDAS anomaly (Table 1). The rms of the differences simulated stripes for the 3000, 4000, 5000 and 6000 km centre-to-centre solutions are also provided (not iterated) for comparison. Computations to this point utilized noise and error free data. We now Table 1 shows that, while iterating will improve the solution, test the effect of adding realistic noise and errors to the simulated the extent of the improvement is limited. For iterated solutions KBRR observations generated from observed noise and errors in convergence is normally reached by ∼5 steps. The rms improves processing actual GRACE data. Although not directly pertinent to with an increase in the centre-to-centre distance of the solution the simulations we outline our current two stage processing strategy for GRACE data as the ﬁrst stage provides the accelerometer bias Table 1. Rms of the differences between the iterated 4 values and state vectors for the short-arc gravity ﬁeld recovery in mascon solutions and the input GLDAS anomaly in cm of the second stage. The accelerometer bias runs are undertaken in EWH. Rms values are provided over land (and globally). 24 hr arcs. In addition to constant accelerometer biases we also The inclusion of a mascon into the solution depends on the estimate the full state vector for GRACE A and B from KBRR, distance between that mascon and the nadir mascon being GNV1B and accelerometer data in addition to KBRR bias, trend lower than the centre-to-centre distance. and once and twice per rev per orbital revolution (16 sets per arc) Centre-to-centre Iteration Rms difference to GLDAS measurement parameters. This process essentially computes the distance (km) number anomaly (cm EWH) acceleration biases from the GNV1B Cartesian coordinates while 0 1 10.5 (6.4) the KBRR parameters account for unmodelled force effects enabling 2 10.3 (6.1) bad data to be identiﬁed. 3 9.9 (6.0) To generate the KBRR noise and errors we re-processed the ac- 4 9.8 (5.9) celerometer runs estimating now 96 sets of empirical once per rev 5 9.7 (6.0) acceleration parameters along-track and cross-track instead of the 6 9.8 (6.0) KBRR parameters. The empirical accelerations absorbed the longer 7 9.8 (6.2) wavelength gravity ﬁeld mismodelling with the KBRR residuals 8 9.9 (6.3) representing a realistic sample of the noise and errors present in the 1000 1 10.4 (6.2) GRACE data. The empirical acceleration once per rev parameters 2 10.1 (6.2) also lead to slightly different values for the accelerometer biases 3 9.9 (6.2) values and estimated state vectors. The small differences in the ac- 4 9.8 (6.2) 5 9.8 (6.3) celerometer bias values between the runs are assumed to be similar 6 9.8 (6.4) to the differences between the original values and the unknown 7 9.8 (6.5) ‘true’ values. Similarly the orbit used to calculate the KBRR resid- 8 9.9 (6.7) uals will be different to the orbits used in the gravity recovery due 2000 1 9.3 (5.9) to the use of different empirical parameters. The KBRR residuals 2 9.1 (5.8) time-series will also contain errors in the orbit introduced from the 3 8.8 (5.6) inclusion of the GNV1B data as well as the KBRR data (GNV1B 4 8.7 (5.5) data is not used in the 60 min arcs for gravity recovery). Further- 5 8.6 (5.5) more, the KBRR residuals will also contain effects of aliasing of 6 8.6 (5.5) 7 8.6 (5.5) short period signals and errors in the background models which will 8 8.6 (5.5) affect the real GRACE solution. AOD1B (Atmosphere and Ocean 3000 1 8.3 (5.2) de-aliasing Level-1B) product and the GLDAS anomaly ﬁeld were 4000 1 7.7 (4.7) utilized in all computations. 5000 1 7.6 (4.6) Despite best efforts it is likely that residual geographical cor- 6000 1 7.5 (4.5) related terms are present in the KBRR residuals due to the All mascons 1 7.4 (4.4) parametrization failing to absorb all gravity ﬁeld mismodelling. Downloaded from https://academic.oup.com/gji/article/200/1/503/752750 by DeepDyve user on 20 July 2022 Mass change from GRACE: simulated comparison 511 Figure 9. Periodogram of the original KBRR residual time-series (a); periodogram created from a FFT of the phase-shifted reconstructed KBRR time-series (b); sample of the original KBRR residual time-series (c); sample of the phase-shifted reconstructed KBRR residual time-series (d). In mitigation we preserved the characterization of the error signa- tures but redistributed the signatures geographically by applying a fast Fourier transform (FFT) to the residual time-series. Modiﬁed KBRR residuals for each 24 hr were calculated by reconstructing the FFT through keeping the same amplitude of spectral frequency but assigning a random phase. Fig. 9 shows the original and modiﬁed FFT and KBRR residuals. Fig. 9(a) shows the original periodogram, while Fig. 9(b) shows the modiﬁed spectrum used to generate the required KBRR residuals. Fig. 9(c) shows a sample of the orig- inal KBRR residuals, while Fig. 9(d) shows the modiﬁed KBRR residuals. The original and modiﬁed KBRR residuals have simi- lar range and noise characteristics. Fig. 9(a) shows features that are replicated in Fig. 9(b). The low frequencies (i.e. long spatial wavelengths) were investigated and characterized as ﬂicker noise. Higher frequencies were seen to be approximately Gaussian with an increase in power between 1800 and 1900 Hz. The frequency of this signal with increased power corresponds to about twice the separa- tion distance between the GRACE satellites. Without this artefact it would have been possible to generate the spectrum using a com- bination of ﬂicker and Gaussian noise. Tests with other days and months showed that the periodogram was representative. Tracking data, in 60 min arcs, was created with the modiﬁed KBRR residuals and the gravity ﬁeld solution recovered to degree and order 60. The resulting gravity ﬁeld is shown in Fig. 10 which reveals that the recovered solution is dominated by short wavelength, high degree, north–south orientated stripes; Fig. 10(b) shows the same but with 400 km Gaussian smoothing applied. The stripes Figure 10. Simulated stripes: unsmoothed (a); with 400 km Gaussian are typical of those in the monthly GRACE solutions. A plot of smoothing (b). Units of cm of EWH. the degree variance of the simulated stripes in Fig. 11 shows that the simulated stripes begin to dominate the signal between degree 20 simulated stripes above degree 30 was found to be comparable to and 30. Degree 30 is the degree at which regularization is normally those of the unsmoothed CSR solution for 2006 June. Most of the applied to monthly solutions (Bruinsma et al. 2010;Save et al. true geophysical signal is restricted to lower degrees, below degree 2012) as required to reduce the effect of noise and errors in the ∼20 (Wouters & Schrama 2007) but short wavelength signals, as high degree coefﬁcients. The power of the degree variances of the required for hydrological and glaciological basin-level applications, Downloaded from https://academic.oup.com/gji/article/200/1/503/752750 by DeepDyve user on 20 July 2022 512 S. B. Andrews, P. Moore and M. A. King rors with consequences that mirror the errors in actual monthly GRACE ﬁelds. The degree variances of the unconstrained solutions in Fig. 11, show that the spherical harmonic and 2 mascon solution are able to reproduce the GLDAS anomaly with simulated stripes. The global values of the 0 km (50 km) solution in Table 2 (spatial constraint mascon solutions) conﬁrm that the simulated stripes are present in the spherical harmonic and 2 mascon solutions, with rms values of 38.1 and 35.2 cm, respectively. The 4 mascon solution is able to reproduce the degree variance of the GLDAS plus simulated stripes to degree ∼35, at which point it loses power. This loss of power, a consequence of undersampling of the signal, results in a rms of 14.7 cm in Table 2 for the 0 km solution. Again, the dis- tance for the spherical harmonic coefﬁcients refers to the Gaussian smoothing radius, while for the mascon solutions the distance is D of eq. (9). As previous, the distances are not directly comparable. Regardless of the strategy used to recover the simulated ﬁeld, some form of post-processing is required to reduce the effects of noise and errors in the solution. Increasing the Gaussian smoothing Figure 11. Degree variances of the simulated stripes, GLDAS with simu- distance for the spherical harmonic solution reduces the rms towards lated stripes added and the unﬁltered spherical harmonic and unconstrained the original GLDAS anomaly, with 400 km Gaussian smoothing mascon solutions. The original GLDAS model is included for comparison. producing the lowest rms to the GLDAS anomaly (Table 2). The global rms still reduces for 500 km Gaussian smoothing due to smoothing over the oceans, but by 500 km the land signal is over is contained in the higher degrees (Kusche 2007). This simulation smoothed. Applying a lower smoothing radius allows the simulated thus allows us to assess the ability of the different methodologies to stripes to dominate the solution, while a larger smoothing radius maintain the signal at lower degrees while reducing noise and error reduces the amplitude of geophysical signals and removes short in the higher degrees and maximizing recovery at these shorter ◦ ◦ wavelength information. For both the 4 and 2 mascon solutions, wavelengths. increasing the correlation distance D in the constraint matrix reduces The modiﬁed KBRR residuals were added to the original track- the rms until ∼700 km. Over all solutions in Table 2 the 2 mascon ing data that described the GLDAS anomaly in Section 4.2. These solution provides the closest match to the GLDAS anomaly when residuals emulate the spectrum of real GRACE noise and er- Table 2. Rms (cm of EWH) of the differences between the calculated gravity ﬁeld solutions and input GLDAS anomalies in the presence of simulated stripes. Rms values are provided over land (globally). The distance refers to the Gaussian smoothing radius (spherical harmonic solution) or the correlation distance (mascon solutions). Mascon solution constrained using spatial (columns 3 and 4) or basin constraints (columns 5 and 6). Rms difference to GLDAS anomaly (cm EWH) Distance (km) Spatial constraint Basin constraint applied ◦ ◦ ◦ ◦ SH coefﬁcients 4 mascon 2 mascon 4 mascon 2 mascon 0 38.1 (37.1) 14.9 (14.7) 36.4 (35.2) – – 50 37.1 (35.9) 14.9 (14.7) 31.4 (30.2) 14.9 (14.7) 31.8 (30.4) 100 31.4 (30.3) 14.1 (13.6) 18.2 (17.3) 14.2 (13.7) 19.3 (17.7) 150 24.3 (23.1) 12.2 (11.1) 13.1 (11.9) 12.6 (11.2) 14.1 (12.3) 200 17.7 (16.4) 10.9 (9.2) 11.1 (9.5) 11.3 (9.3) 11.8 (9.7) 250 13.2 (11.4) 10.1 (8.0) 10.1 (8.1) 10.5 (8.1) 10.7 (8.3) 300 10.9 (8.4) 9.6 (7.3) 9.5 (7.2) 9.9 (7.4) 10.0 (7.4) 350 10.1 (6.9) 9.3 (6.7) 9.1 (6.7) 9.6 (6.8) 9.5 (6.8) 400 10.0 (6.4) 9.1 (6.4) 8.9 (6.3) 9.3 (6.4) 9.2 (6.4) 450 10.1 (6.2) 8.9 (6.1) 8.8 (6.0) 9.2 (6.2) 9.0 (6.1) 500 10.3 (6.2) 8.9 (5.9) 8.7 (5.8) 9.0 (5.9) 8.9 (5.9) 550 10.5 (6.3) 8.8 (5.8) 8.6 (5.6) 8.9 (5.8) 8.8 (5.7) 600 10.7 (6.3) 8.8 (5.7) 8.6 (5.5) 8.9 (5.6) 8.8 (5.5) 650 10.9 (6.4) 8.8 (5.6) 8.6 (5.4) 8.8 (5.5) 8.7 (5.4) 700 11.1 (6.5) 8.8 (5.5) 8.6 (5.4) 8.8 (5.5) 8.7 (5.4) 750 11.3 (6.5) 8.8 (5.5) 8.7 (5.3) 8.8 (5.4) 8.7 (5.3) 800 11.5 (6.6) 8.8 (5.4) 8.7 (5.3) 8.8 (5.3) 8.7 (5.3) 850 11.6 (6.7) 8.8 (5.4) 8.7 (5.3) 8.8 (5.3) 8.7 (5.2) 900 11.8 (6.7) 8.9 (5.4) 8.7 (5.2) 8.8 (5.3) 8.7 (5.2) 950 11.9 (6.8) 8.9 (5.4) 8.8 (5.2) 8.8 (5.3) 8.7 (5.2) 1000 12.0 (6.9) 8.9 (5.3) 8.8 (5.2) 8.8 (5.2) 8.7 (5.2) 1050 12.1 (6.9) 8.9 (5.3) 8.8 (5.2) 8.8 (5.2) 8.7 (5.1) 1100 12.2 (7.0) 9.0 (5.3) 8.9 (5.2) 8.8 (5.2) 8.8 (5.1) Downloaded from https://academic.oup.com/gji/article/200/1/503/752750 by DeepDyve user on 20 July 2022 Mass change from GRACE: simulated comparison 513 exploited the full advantage of the mascon processing strategy; namely that the constraint matrix in eq. (5) can be adapted to con- strain mascons that share geophysical similarities as in Luthcke et al. (2006b). Thus, in addition to separating mascons into either land or ocean, we now also separate the major hydrological drainage basins deﬁned by the total runoff integrating pathways (TRIP) data set(Oki&Sud 1998). Major Greenland and Antarctic drainage basins are adopted according to Zwally et al. (2012). The result is a constraint matrix where constraints are imposed between all land mascons within the same (hydrological or glaciological) drainage basin, with mascons in different basins uncorrelated. Not all land mascons are within a major hydrological or glaciological basin. Land mascons that are not classiﬁed as belonging to a basin were assumed to be correlated, but in most instances d of eq. (8) or (9) IJ reduces the correlation to ∼0. All ocean mascons are correlated as previous. Table 2 (Basin constraint mascon solutions) summarizes the rms of the differences to the GLDAS plus simulated stripes. For both Figure 12. The degree variances of GLDAS (with and without simulated mascon solutions the use of a basin constraint, while offering a ◦ ◦ stripes) and the spherical harmonic and mascon solutions. 2 and 4 mascon similar or small improvement in the rms over land, offers an im- are the constrained mascon solutions without the basin constraints applied provement in global rms values over the solutions without basin ◦ ◦ while 2 and 4 basin are the constrained mascon solutions with basin constraints. An increase in the correlation distance (800–900 km) constraints applied. The spherical harmonic solution has 400 km Gaussian is required to obtain these improvements due to the reduced num- smoothing applied; the mascon solutions are constrained using a correlation ber of constraints based on mascons pairs that share geophysi- distance of 700 km in the constraint matrix (without basin constraints) and 900 km (with basin constraint). cal similarities. The improved global rms is due to reduction of the simulated stripes over the ocean. As ocean mascons are con- a correlation distance of ∼650 km is used in the constraint matrix. strained together, the increased correlation distance allows the effect For the 4 mascon solution we investigated the effect on the weight of the simulated stripes to be further mitigated. The improvement matrix V of varying the value of S in the eq. (9). The weight matrix is evident in the degree variances of the mascon solutions (basin) V is dependent on S and D.Using S = 0.001, the lowest rms (8.8 cm plotted in Fig. 12. Comparing mascon solutions with and without EWH) was achieved using D = 700 km. When S = 0.01 (0.0001) basin constraints, the power of the resulting mascon solutions from were used rms values of 9.8 (10.3) were obtained. For S = 0.01 the degrees 30 to 60 more closely matches the input signal when the lowest rms was 8.8 cm when D = 250 km, while for S = 0.0001 basin constraint is applied. The constraints have successfully re- an rms of 9.0 cm was obtained when D = 1500 km. We therefore duced the noise at these degrees while reducing loss of actual signal. concluded, within the scope of this limited testing, that S = 0.001 Retention of short wavelength signal is of fundamental importance to be the optimal scale factor for the simulations and will use this for improving basin-scale hydrological and glaciological recovery is all further processing. and for estimating basin-level mass ﬂux. Fig. 12 shows the degree variances of the spherical harmonic so- The ﬁnal recovered GLDAS anomalies are shown in Fig. 13 for lution with 400 km Gaussian smoothing. It is not possible to directly the spherical harmonic solution with 400 km Gaussian smoothing relate the Gaussian smoothing distance to the mascon correlation together with the 2 mascon solution using a 800 km correlation dis- distance so we show the mascon solutions for D = 700 km. This tance and basin constraints. The input GLDAS anomaly is provided value was chosen as a comparison of the degree variances of the for reference. Comparing Figs 13(c) to (a) validates the 2 mascon solutions given in Table 2 revealed it removes a comparable amount solution. The mascon methodology offers a clear improvement in of noise to the spherical harmonic solution at the higher degrees. the recovery of TVG and suppression of simulated stripes compared It also gives one of the lowest rms values to the GLDAS anomaly. with a spherical harmonic solution with Gaussian smoothing. The Fig. 12 reveals that all solutions have removed a similar amount of reduction of the simulated stripes over the ocean is particularly evi- noise and error above degree 30, and contain less power at these de- dent. Unlike spherical harmonics, mascons allow constraints within grees compared to the GLDAS signal. The mascon solutions more areas that share real world geophysical properties. Any number of closely resemble the power of the GLDAS signal between ∼degree constraints can be applied with the aim of improving the overall 15 and 30 than the spherical harmonic solution. The reduced power recovery of the TVG ﬁeld. of the spherical harmonic signal between ∼degree 15 and 30 could be partly related to the degree dependency of the Gaussian ﬁlter which Han et al. (2005) suggested could suppress information over 4.7 Simulated recovery of Antarctic mass change signal this range. The mascon solution has therefore been able to recover more of the original signal at degrees ∼15–30 while removing a One of the scientiﬁc advances provided by the GRACE mission is comparable amount of noise and error at higher degrees. recovery of continent-wide ice sheet mass change estimates, such as over Antarctica, with temporal resolutions of 30 days or less. So far our simulations have focused on recovery of a GLDAS hydro- 4.6 Basin constraints logical anomaly. Given the scientiﬁc interest in contribution of ice Mascon solutions with spatial constraints are able to remove more sheets to sea-level change and the convergence of orbital tracks at noise and error and retain more signal than spherical harmonics the pole we now simulate recovery of a realistic Antarctica mass solutions with Gaussian smoothing. However, thus far we have not signal along with the GLDAS anomaly and error/noise product Downloaded from https://academic.oup.com/gji/article/200/1/503/752750 by DeepDyve user on 20 July 2022 514 S. B. Andrews, P. Moore and M. A. King Figure 13. GLDAS anomaly for 2006 June to degree and order 120 (a); Re- covered GLDAS anomalies for the spherical harmonic solution with 400 km Gaussian smoothing to degree and order 60 (b); and 2 mascon solution with a basin constraint applied, with 800 km correlation distance (c). Units of cm EWH. described in Section 4.5. The simulated Antarctica mass signal was derived from the basin mass change rates provided in Table S1 of the supplementary material of King et al. (2012). While King et al. (2012) calculated GRACE solutions to degree and order 60, forward modelling across deﬁned basins resulted in a higher effective reso- lution. The mass change rates used have been corrected for leakage but still contain the GIA signal that is observable from GRACE. Using the same drainage basins adopted by King et al. (2012), as deﬁned using ICESat data (Zwally et al. 2012), the mass change rates for each basin were converted to a total mass change in cm of EWH using the density of water and the area of each drainage Figure 14. Simulated Antarctic mass change signal (a); spherical harmonic basin. As in King et al. (2012) basins 25 and 26, which cover the solution ﬁltered with a 350 km Gaussian smoother (b); and 2 mascon solu- northern tip of the Antarctica Peninsula, were merged into a single tion with basin constraint with correlation distance D = 800 km (c). Units basin. The resulting mass changes for each basin were then ex- of cm of EWH. panded into spherical harmonics to degree and order 60 and added to the mean background ﬁeld to create the KBRR tracking data. Simulated stripes are included as previous. The simulated Antarctic mass change signal can be seen in Fig. 14. Plotting the unﬁltered Downloaded from https://academic.oup.com/gji/article/200/1/503/752750 by DeepDyve user on 20 July 2022 Mass change from GRACE: simulated comparison 515 Table 3. Rms (cm of EWH) of the differences between the calculated gravity ﬁeld solutions and the input GLDAS anomaly with Antarctic mass change signal in the presence of different simulated noise products. Rms values are provided over land (and globally). The distance refers to the Gaussian smooth- ing radius (spherical harmonic coefﬁcient solutions) or the correlation distance (mascon solutions). All mascon solutions utilize basin constraints described in Section 4.6. For the thermal mascon solutions the value S in eq. (9) was increased from 0.001 to 0.01. Rms difference to GLDAS anomaly and Antarctic mass change (cm EWH) Distance (km) Main period Resonance period Reduced thermal control ◦ ◦ ◦ SH 2 mascon SH 2 mascon SH 2 mascon 200 17.7 (16.4) 11.8 (9.8) 20.3 (19.8) 8.1 (5.3) 45.2 (45.5) 13.9 (11.0) 250 13.2 (11.4) 10.7 (8.3) 14.8 (13.8) 8.0 (5.0) 29.9 (29.7) 12.4 (9.3) 300 10.9 (8.4) 10.0 (7.4) 11.6 (9.8) 8.0 (4.9) 20.6 (19.5) 11.7 (8.2) 350 10.1 (7.0) 9.5 (6.8) 10.3 (7.6) 8.0 (4.9) 15.6 (13.7) 11.3 (7.6) 400 10.0 (6.4) 9.2 (6.4) 10.0 (6.6) 8.1 (4.9) 13.2 (10.6) 11.0 (7.2) 450 10.2 (6.3) 9.0 (6.1) 10.1 (6.2) 8.2 (4.9) 12.2 (9.0) 10.9 (7.0) 500 10.4 (6.3) 8.9 (5.9) 10.3 (6.2) 8.3 (4.9) 11.7 (8.1) 10.8 (6.8) 550 10.6 (6.3) 8.8 (5.7) 10.6 (6.2) 8.3 (4.9) 11.6 (7.6) 10.8 (6.6) 600 10.8 (6.4) 8.8 (5.5) 10.8 (6.3) 8.4 (4.9) 11.6 (7.4) 10.8 (6.6) 650 11.0 (6.4) 8.7 (5.4) 11.1 (6.4) 8.5 (5.0) 11.6 (7.2) 10.8 (6.5) 700 11.2 (6.5) 8.7 (5.4) 11.3 (6.5) 8.6 (5.0) 11.7 (7.2) 10.8 (6.5) 750 11.4 (6.6) 8.7 (5.3) 11.5 (6.6) 8.6 (5.0) 11.8 (7.1) 10.9 (6.4) 800 11.6 (6.7) 8.7 (5.3) 11.7 (6.7) 8.7(5.0) 11.9 (7.1) 10.9 (6.4) 850 11.7 (6.7) 8.7 (5.2) 11.8 (6.8) 8.8 (5.1) 12.0 (7.1) 11.0 (6.4) 900 11.9 (6.8) 8.7 (5.2) 12.0 (6.9) 8.8 (5.1) 12.1 (7.1) 11.0 (6.4) 950 12.0 (6.9) 8.7 (5.2) 12.1 (6.9) 8.9 (5.1) 12.2 (7.1) 11.1 (6.4) 1000 12.1 (6.9) 8.7 (5.2) 12.2 (7.0) 9.0 (5.2) 12.3 (7.2) 11.1 (6.4) 1050 12.2 (7.0) 8.7 (5.1) 12.3 (7.0) 9.0 (5.2) 12.3 (7.2) 11.2 (6.5) 1100 12.3 (7.0) 8.8 (5.1) 12.4 (7.1) 9.1 (5.2) 12.4 (7.2) 11.2 (6.5) degree variances (not shown) reveals that the spherical harmonic 4.8 GRACE resonance and thermal control issues and 2 mascon solution were able to reproduce the Antarctica mass Periods where the GRACE satellites are in 3, 4 and 7 day repeat signal added to the GLDAS anomaly and simulated stripes. Appli- orbits modify the noise characteristics of the GRACE solutions cation of a single value of mass change per basin is a simpliﬁcation (Save et al. 2012), while recent GRACE solutions are affected by on the real spatial pattern of change, but this assumption is sufﬁcient poor satellite thermal control. For replicating the GRACE errors for our purposes here. and studying the correlations between the spherical harmonic coef- Solutions were estimated for the spherical harmonics and 2 ﬁcients in CSR RL05 solutions, covariance matrices are provided mascon methodologies with the rms values shown in Table 3 by CSR for the three characteristic periods in the GRACE mission: (main period). The drainage basins used in the mascon solution in early mission (2002 February–2005 May), main part of the mission Section 4.6 result in mascons within the same Antarctic drainage (2007 July–2010 December) and recent months (2011 February– ) basin being constrained together, but with mascons in different which are affected by poor thermal control. Only three covariance Antarctic basins uncorrelated. The mascon results are shown in matrices are required as the monthly variations within these periods Fig. 14(c) revealing this solution has been able to recover most of are minor as to be unimportant (John C Ries, personal communica- the input signal (Fig. 14a) while the spherical harmonic solution tion, 2014). An additional covariance matrix is available for months ﬁltered with 400 km Gaussian smoothing (Fig. 14b) has recovered where the satellite is in a repeat orbit (resonance). The effect these the spatial extent of the signal but with reduced amplitudes. This different periods have on the noise and errors is investigated through is important to consider within future GRACE studies of Antarctic recovery of the GLDAS anomaly and Antarctica mass change sig- ice mass change. The lowest rms of the differences to the input nal for two additional sets of simulated noise. The GLDAS anomaly signal (in EWH) over Antarctica is 8.7 cm for the mascon solution and Antarctica mass change signal are the same as in previous sim- (800 km correlation distance) and 10.4 cm with spherical harmonics ulations. The ﬁrst set of simulated noise relates to September 2004, (400 km Gaussian smoothing). As with the global hydrological sim- when the GRACE satellites were in 61:4 resonance (Klokocn ˇ ´ık et al. ulations, we are able to recover more of the signal using the mascon 2013). The second set relates to 2011 March, which corresponds to methodology. a period when the GRACE satellites are affected by poor thermal As a ﬁnal comparison we computed the mascon solution with control. KBRR residuals for 2004 September and 2011 March with only land/ocean constraints as applied in Section 4.5. The rms of the simulated noise were derived using the method of Section 4.5. We difference to the input Antarctica signal was 8.8 cm. This highlights computed the modiﬁed FFT and KBRR residuals for 2004 Septem- the advantage offered by the mascon methodology when real-world ber and 2011 March (as in Fig. 9). For 2004 September (not shown) geophysical properties are used to deﬁne the constraints between these revealed that the noise and errors in GRACE during resonance the mascons. We note, however, that a simple basin-wide constraint are indistinguishable from that of the main part of the mission. The may not be an accurate reﬂection of changes within a basin and more 2011 March residuals are shown in Fig. 15, revealing an increase in sophisticated treatment, such as increasing the numbers of basins the amplitude of the KBRR residuals; compare Figs 15(c) and (d) (e.g. Luthcke et al. 2013) based on ice dynamic considerations, may with Figs 9(c) and (d). This increase in amplitude is the result of further improve analysis of real GRACE data. Downloaded from https://academic.oup.com/gji/article/200/1/503/752750 by DeepDyve user on 20 July 2022 516 S. B. Andrews, P. Moore and M. A. King Figure 15. Periodogram of the poor thermal control KBRR residual time-series (a) and FFT phase shifted reconstructed KBRR time-series (b); sample of the poor thermal control KBRR residual time-series (c); sample of the reconstructed KBRR residual time-series (d). the poor thermal control, which we would expect to increase from S in eq. (9) was increased from 0.001 to 0.01 which increases the 2011 onwards. effect the constraint matrix has on the solution. For both the mascon Table 3 summarizes the rms of the difference to the GLDAS and and spherical harmonic solutions a larger constraint matrix is re- Antarctica mass loss signal in the presence of simulated stripes from quired to deal with the noise and errors resulting from the poor three periods (main, resonance and thermal control). The simulated thermal control. While the recovery in the presence of the thermal stripes for the main period are described in Section 4.5. The rms stripes is not as favourable for both methodologies, when compared values are provided for the spherical harmonic and 2 mascon solu- to the original or resonance periods, the mascon methodology is still tions. In the resonance period, increasing the Gaussian smoothing able to recover more of the input signal than the spherical harmonic until 400 km reduces the rms at which point the signal recovered methodology while also reducing the noise over the ocean. is comparable to the main period; the larger global rms suggests more of the stripes remain over the oceans. For the mascon solution the lowest rms obtained is obtained when a correlation distance of 5CONCLUSION 250 km is used, which is a much lower than the distance required for the 2006 June simulation of the main period. The rms is also much Through the use of simulations we have demonstrated the sensitivity lower, suggesting the mascon methodology is able to deal well with of GRACE to mass over a large area, with the effect on the KBRR periods of resonance. In 61:4 resonance the satellite ground-track residuals depending on the amplitude of the signal. Building on repeats every 61 orbital revolutions in the time the Earth revolves this we showed that recovery of global mass distributions require 4 times relative to the ascending node. This gives 122 repeating as- contributions to the satellite accelerations of all mascons within cending and descending passes every 4 days spaced on average every approximately 6000 km of the nadir mascon. This has implications 3 in longitude. This spatial resolution is adequate for a spherical when estimating a non-global mascon solution, and is useful in harmonic ﬁeld to degree and order 60. Similarly, a mascon correla- reducing computational effort for high resolution global mascon tion distance of 250 km, equivalent to the distance between centres solutions. We found that the use of iterations cannot compensate for of 2 equatorial mascons, provides the ties between neighbouring an insufﬁcient distribution of mascons. mascons. It would appear that the resonance facilitates separation Simulated recovery of a noise and error free GLDAS mass of spatial and temporal errors in our simulation leading to more anomaly showed that there is no quantiﬁable difference between accurate surface mass recovery than in non-resonant periods. a mass anomaly solution recovered from KBRR data using spheri- In the period of reduced thermal control both the spherical har- cal harmonics and a 2 equal area mascon solution; both recovered monic and mascon solutions are degraded. The lowest rms (11.6 cm the input GLDAS anomaly. A 4 mascon solution does not have the EWH) for the spherical harmonic solution is obtained with 600 km required resolution to recover the gravity ﬁeld solution to degree Gaussian smoothing. Even with this smoothing radius, the rms and order 60 due to the number of mascon parameters estimated obtained is larger than for the main and resonance periods. For being less than the gravity ﬁeld parameters required; this causes a the mascon solution, the lowest rms (10.8) is obtained with a corre- loss of higher degree information in the gravity ﬁeld. Recovery of lation distance D = 600 km, compared to 8.7 (8.0) for the solutions the GLDAS mass anomaly described to degree and order 120 re- with the original (resonance) noise. To achieve this rms the value vealed that all methodologies are affected by omission errors, with Downloaded from https://academic.oup.com/gji/article/200/1/503/752750 by DeepDyve user on 20 July 2022 Mass change from GRACE: simulated comparison 517 the additional short wavelength (high degree) information not being Chen, J.L., Wilson, C.R., Tapley, B.D., Blankenship, D. & Young, D., 2008. Antarctic regional ice loss rates from GRACE, Earth planet. Sci. Lett., captured. 266, 140–148. The recovery of the GLDAS anomaly with simulated stripes Dunn, C. et al., 2003. Instrument of GRACE: GPS augments gravity mea- reveals the advantage of the mascon solution over a spherical har- surements, GPS World, 14, 16–28. monic solution smoothed with a Gaussian ﬁlter when attempting to Farrell, W.E., 1972. Deformation of earth by surface loads, Rev. Geophys. reduce high degree and order noise. The use of the mascon spatial Space Phys., 10, 761–797. constraint matrix allowed more signal to be preserved up to degree Han, S.-C., Shum, C.K., Jekeli, C., Kuo, C.-Y., Wilson, C. & Seo, K.-W., ∼35, while the addition of constraints between mascon parameters 2005. Non-isotropic ﬁltering of GRACE temporal gravity for geophysical that share geophysical similarities resulted in a further reduction signal enhancement, Geophys. J. Int., 163(1), 18–25. of the lost signal at all degrees. This improvement is further evi- Han, S.-C., Sauber, J., Luthcke, S.B., Ji, C. & Pollitz, F.F., 2008. Implica- dent in the simulated recovery of the Antarctic mass signal which tions of postseismic gravity change following the great 2004 Sumatra- Andaman earthquake from the regional harmonic analysis of GRACE also conﬁrms the advantage of the mascon solution over the spher- intersatellite tracking data, J. geophys. Res.: Solid Earth, 113, B11413, ical harmonic recovery validating the use of this methodology in doi:10.1029/2008JB005705. polar regions. Inclusion of temporal constraints could result in fur- Han, S.-C., Kim, H., Yeo, I.-Y., Yeh, P., Oki, T., Seo, K.-W., Alsdorf, D. & ther improvements, but a level of assumption would be required on Luthcke, S.B., 2009. Dynamics of surface water storage in the Amazon plausible variations of mass over time. inferred from measurements of inter-satellite distance change, Geophys. The recovery of the GLDAS anomaly and Antarctica mass change Res. Lett., 36, L09403, doi:10.1029/2009GL037910. signal with simulated stripes, created during periods when GRACE Jacob, T., Wahr, J., Pfeffer, W.T. & Swenson, S., 2012. Recent contributions was in resonance and suffering from poor thermal control, revealed of glaciers and ice caps to sea level rise, Nature, 482, 514–518. that in both cases the mascon methodology was able to recover more King, M.A., Bingham, R.J., Moore, P., Whitehouse, P.L., Bentley, M.J. & of the input signal than using spherical harmonic methodology. This Milne, G.A., 2012. Lower satellite-gravimetry estimates of Antarctic sea- level contribution, Nature, 491, 586–589. advantage was especially evident during periods of resonance. The Klees, R., Revtova, E.A., Gunter, B.C., Ditmar, P., Oudman, E., Winsemius, use of the mascon methodology also reduced the effect of stripes H.C. & Savenije, H.H.G., 2008. The design of an optimal ﬁlter for monthly over the ocean. However, both the spherical harmonic and mascon GRACE gravity models, Geophys. J. Int., 175, 417–432. solutions were degraded during periods of reduced thermal control. Klokocn ˇ ´ık, J., Gooding, R.H., Wagner, C.A., Kostelecky, ´ J. & Bezdek, ˇ A., In particular, the reduced thermal control necessitated tighter con- 2013. The use of resonant orbits in satellite geodesy: a review, Surv. straints between the various mascons than during the main part of Geophys., 34, 43–72. the mission. Klosko, S., Rowlands, D., Luthcke, S., Lemoine, F., Chinn, D. & Rodell, Finally, while the values used for S and D were found to be optimal M., 2009. Evaluation and validation of mascon recovery using GRACE in this paper, they may differ for users depending on parametrization, KBRR data with independent mass ﬂux estimates in the Mississippi Basin, and are therefore provided as a guide. Furthermore, adaptations to J. Geod., 83, 817–827. Kusche, J., 2007. Approximate decorrelation and non-isotropic smoothing of the S and D values appear to be necessary as the GRACE mission time-variable GRACE-type gravity ﬁeld models, J. Geod., 81, 733–749. ages. Lemoine, F.G., Luthcke, S.B., Rowlands, D.D., Chinn, D.S., Klosko, S.M. & Cox, C.M., 2007. The use of mascons to resolve time-variable gravity from GRACE, in Dynamic Planet, pp. 231–236, eds Tregoning, P. & ACKNOWLEDGEMENTS Rizos, C., Springer. This work was funded by a National Centre for Earth Observation Luthcke, S.B., Rowlands, D.D., Lemoine, F.G., Klosko, S.M., Chinn, D. & (NCEO) award in gravimetric geodesy and partially supported by McCarthy, J.J., 2006a. Monthly spherical harmonic gravity ﬁeld solutions COST Action ES0701 ‘Improved constraints on models of Glacial determined from GRACE inter-satellite range-rate data alone, Geophys. Res. Lett., 33, L02402, doi:10.1029/2005GL024846. Isostatic Adjustment’. We are gratefully for the GRACE Level- Luthcke, S.B. et al., 2006b. Recent Greenland ice mass loss by 1B products produced by the Jet Propulsion Laboratory. MAK is drainage system from satellite gravity observations, Science, 314, a recipient of an Australian Research Council Future Fellowship 1286–1289. (project number FT110100207). The authors would like to thank Luthcke, S.B., Arendt, A.A., Rowlands, D.D., McCarthy, J.J. & Larsen, three anonymous reviewers for comments which helped improve C.F., 2008. Recent glacier mass changes in the Gulf of Alaska region this manuscript. from GRACE mascon solutions, J. Glaciol., 54, 767–777. Luthcke, S.B., Sabaka, T.J., Loomis, B.D., Arendt, A.A., McCarthy, J.J. & Camp, J., 2013. Antarctica, Greenland and Gulf of Alaska land-ice REFERENCES evolution from an iterated GRACE global mascon solution, J. Glaciol., 59, 613–631. Awange, J.L., Fleming, K.M., Kuhn, M., Featherstone, W.E., Heck, B. & Moore, P., Boomkamp, H.J., Carnochan, S. & Walmsley, R.J., 1999. Anjasmara, I., 2011. On the suitability of the 4 degree x 4 degree GRACE FAUST: multi-satellite orbital dynamics software, Adv. Space Res., 23, mascon solutions for remote sensing Australian hydrology, Rem. Sens. 785–795. Environ., 115, 864–875. Muller, P.M. & Sjogren, W.L., 1968. Mascons: lunar mass concentrations, Bettadpur, S., 2012. UTCSR Level-2 Processing Standards Document for Science, 161, 680–684. Level-2 Product Release 0005, pp. 16, GRACE 327–742. University of Oki, T. & Sud, Y.C., 1998. Design of total runoff integrating pathways Texas at Austin, CSR Publ: GR-12-xx, Rev. 4.0. (TRIP)—a global river channel network, Earth Interact., 2, 1–37. Bruinsma, S., Lemoine, J.-M., Biancale, R. & Vales, ` N., 2010. CNES/GRGS Pritchard, H.D., Luthcke, S.B. & Fleming, A.H., 2011. Understanding 10-day gravity ﬁeld models (release 2) and their evaluation, Adv. Space ice-sheet mass balance: progress in satellite altimetry and gravimetry, Res., 45, 587–601. J. Glaciol., 56, 1151–1161. Chambers, D.P. & Bonin, J.A., 2012. Evaluation of release-05 GRACE time- Ramillien, G., Biancale, R., Gratton, S., Vasseur, X. & Bourgogne, S., variable gravity coefﬁcients over the ocean, Ocean Sci., 8, 859–868. 2011. GRACE-derived surface water mass anomalies by energy in- Chambers, D.P., Wahr, J., Tamisiea, M.E. & Nerem, R.S., 2010. Ocean mass tegral approach: application to continental hydrology, J. Geod., 85, from GRACE and glacial isostatic adjustment, J. geophys. Res., 115, 313–328. B11415, doi:10.1029/2010JB007530. Downloaded from https://academic.oup.com/gji/article/200/1/503/752750 by DeepDyve user on 20 July 2022 518 S. B. Andrews, P. Moore and M. A. King Ries, J.C., Bettadpur, S., Poole, S. & Richter, T., 2011. Mean background Seo, K.-W., Wilson, C.R., Chen, J. & Waliser, D.E., 2008. GRACE’s spatial gravity ﬁelds for GRACE processing, in Proceedings of the GRACE aliasing error, Geophys. J. Int., 172, 41–48. Science Team Meeting, Austin, TX. Shepherd, A. et al., 2012. A reconciled estimate of ice-sheet mass balance, Riva, R.E.M. et al., 2009. Glacial isostatic adjustment over Antarctica from Science, 338, 1183–1189. combined ICESat and GRACE satellite data, Earth planet. Sci. Lett., 288, Swenson, S. & Wahr, J., 2002. Methods for inferring regional surface-mass 516–523. anomalies from Gravity Recovery and Climate Experiment (GRACE) Rodell, M. et al., 2004. The global land data assimilation system, Bull. Am. measurements of time-variable gravity, J. geophys. Res., 107, 2193, Meteor. Soc., 85(3), 381–394. doi:10.1029/2001JB000576. Rowlands, D.D., Ray, R.D., Chinn, D.S. & Lemoine, F.G., 2002. Short- Swenson, S. & Wahr, J., 2006. Post-processing removal of corre- arc analysis of intersatellite tracking data in a gravity mapping mission, lated errors in GRACE data, Geophys. Res. Lett., 33, L08402, J. Geod., 76, 307–316. doi:10.1029/2005GL025285. Rowlands, D.D., Luthcke, S.B., Klosko, S.M., Lemoine, F.G.R., Chinn, D.S., Tapley, B.D., Bettadpur, S., Ries, J.C., Thompson, P.F. & Watkins, M.M., McCarthy, J.J., Cox, C.M. & Anderson, O.B., 2005. Resolving mass ﬂux at 2004a. GRACE measurements of mass variability in the Earth system, high spatial and temporal resolution using GRACE intersatellite measure- Science, 305, 503–505. ments, Geophys. Res. Lett., 32, L04310, doi:10.1029/2004GL021908. Tapley, B.D., Bettadpur, S., Watkins, M. & Reigber, C., 2004b. The gravity Rowlands, D.D., Luthcke, S.B., McCarthy, J.J., Klosko, S.M., Chinn, D.S., recovery and climate experiment: mission overview and early results, Lemoine, F.G., Boy, J.P. & Sabaka, T.J., 2010. Global mass ﬂux solutions Geophys. Res. Lett., 31, 503–505. from GRACE: a comparison of parameter estimation strategies—mass Velicogna, I. & Wahr, J., 2006. Measurements of time-variable gravity show mass loss in Antarctica, Science, 311, 1754–1756. concentrations versus Stokes coefﬁcients, J. geophys. Res., 115, B01403, doi:10.1029/2009JB006546. Wahr, J., Molenaar, M. & Bryan, F., 1998. Time variability of the Earth’s Sabaka, T.J., Rowlands, D.D., Luthcke, S.B. & Boy, J.P., 2010. Improving gravity ﬁeld: hydrological and oceanic effects and their possible detection global mass ﬂux solutions from Gravity Recovery and Climate Experi- using GRACE, J. geophys. Res., 103, 30 205–30 229. ment (GRACE) through forward modeling and continuous time correla- Wahr, J., Swenson, S. & Velicogna, I., 2006. Accuracy of GRACE mass tion, J. geophys. Res., 115, B11403, doi:10.1029/2010JB007533. estimates, Geophys. Res. Lett., 33, L06401, doi:10.1029/2005GL025305. Save, H., Bettadpur, S. & Tapley, B., 2012. Reducing errors in the GRACE Wahr, J., Swenson, S., Zlotnicki, V. & Velicogna, I., 2004. Time-variable gravity solutions using regularization, J. Geod., 86, 695–711. gravity from GRACE: ﬁrst results, Geophys. Res. Lett., 31, L11501, Schrama, E.J.O. & Wouters, B., 2011. Revisiting Greenland ice sheet doi:10.1029/2004GL019779. mass loss observed by GRACE, J. geophys. Res., 116, B02407, Wouters, B. & Schrama, E.J.O., 2007. Improved accuracy of GRACE gravity doi:10.1029/2009JB006847. solutions through empirical orthogonal function ﬁltering of spherical har- Schrama, E.J.O., Wouters, B. & Rietbroek, R., 2014. A mascon ap- monics, Geophys. Res. Lett., 34, L23711, doi:10.1029/2007GL032098. proach to assess ice sheet and glacier mass balances and their uncer- Wouters, B., Chambers, D. & Schrama, E.J.O., 2008. GRACE observes tainties from GRACE data, J. geophys. Res.: Solid Earth, 119, 6048– small-scale mass loss in Greenland, Geophys. Res. Lett., 35, L20501, 6066. doi:10.1029/2008GL034816. Seo, K.W., Wilson, C.R., Famiglietti, J.S., Chen, J.L. & Rodell, M., Zwally, H.J., Giovinetto, M.B., Beckley, M.A. & Sab, J.L., 2012. 2006. Terrestrial water mass load changes from Gravity Recovery Antarctic and Greenland drainage systems, GSFC Cryospheric Sci- and Climate Experiment (GRACE), Water Resour. Res., 42, W05417, ences Laboratory, Available at: http://icesat4.gsfc.nasa.gov/cryo_data/ doi:10.1029/2005WR004255. ant_grn_drainage_systems.php, last accessed 9 December 2013. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press http://www.deepdyve.com/lp/oxford-university-press/mass-change-from-grace-a-simulated-comparison-of-level-1b-analysis-8ArkJ3QpjH

Loading next page...

References (55)

(2011)
Mean background gravity fields for GRACE processing
B. Wouters, Don Chambers, E. Schrama (2008)
GRACE observes small‐scale mass loss in Greenland
Geophysical Research Letters, 35
S. Luthcke, H. Zwally, W. Abdalati, D. Rowlands, R. Ray, R. Nerem, F. LeMoine, J. Mccarthy, D. Chinn (2006)
Recent Greenland Ice Mass Loss by Drainage System from Satellite Gravity Observations
Science, 314
Ki-Weon Seo, C. Wilson, J. Famiglietti, Jianli Chen, M. Rodell (2006)
Terrestrial water mass load changes from Gravity Recovery and Climate Experiment (GRACE)
Water Resources Research, 42
E. Schrama, B. Wouters, R. Rietbroek (2014)
A mascon approach to assess ice sheet and glacier mass balances and their uncertainties from GRACE data
Journal of Geophysical Research: Solid Earth, 119
Shin‐Chan Han, J. Sauber, S. Luthcke, C. Ji, F. Pollitz (2008)
Implications of postseismic gravity change following the great 2004 Sumatra-Andaman earthquake from the regional harmonic analysis of GRACE intersatellite tracking data
Journal of Geophysical Research, 113
J. Wahr, S. Swenson, I. Velicogna (2006)
Accuracy of GRACE mass estimates
Geophysical Research Letters, 33
C. Dahle, F. Flechtner, C. Gruber, D. König, R. König, G. Michalak, K. Neumayer (2012)
GFZ GRACE Level-2 Processing Standards Document for Level-2 Product Release 0005
T. Oki, Y. Sud (1998)
Design of Total Runoff Integrating Pathways (TRIP)—A Global River Channel Network
Earth Interactions, 2
D. Rowlands, R. Ray, D. Chinn, F. Lemoine (2002)
Short-arc analysis of intersatellite tracking data in a gravity mapping mission
Journal of Geodesy, 76
E. Schrama, B. Wouters (2011)
Revisiting Greenland ice sheet mass loss observed by GRACE
Journal of Geophysical Research, 116
Matt King, R. Bingham, P. Moore, P. Whitehouse, M. Bentley, G. Milne (2012)
Lower satellite-gravimetry estimates of Antarctic sea-level contribution
Nature, 491
S. Swenson, J. Wahr (2006)
Post‐processing removal of correlated errors in GRACE data
Geophysical Research Letters, 33
D. Chambers, J. Bonin (2012)
Evaluation of Release-05 GRACE time-variable gravity coefficients over the ocean
Ocean Science, 8
Shin‐Chan Han, Hyungjun Kim, I. Yeo, P. Yeh, T. Oki, K. Seo, D. Alsdorf, S. Luthcke (2009)
Dynamics of surface water storage in the Amazon inferred from measurements of inter‐satellite distance change
Geophysical Research Letters, 36
(2012)
Antarctic and Greenland drainage systems
G. Ramillien, R. Biancale, S. Gratton, X. Vasseur, S. Bourgogne (2011)
GRACE-derived surface water mass anomalies by energy integral approach: application to continental hydrology
Journal of Geodesy, 85
Don Chambers, J. Wahr, M. Tamisiea, R. Nerem (2010)
Ocean mass from GRACE and glacial isostatic adjustment
Journal of Geophysical Research, 115
(2012)
UTCSR Level - 2 Processing Standards Document for Level - 2 Product Release 0005 , pp . 16 , GRACE 327 – 742
M. Rodell, P. Houser, U. Jambor, J. Gottschalck, K. Mitchell, C. Meng, K. Arsenault, B. Cosgrove, J. Radakovich, M. Bosilovich, J. Entin, J. Walker, D. Lohmann, D. Toll (2004)
THE GLOBAL LAND DATA ASSIMILATION SYSTEM
Bulletin of the American Meteorological Society, 85
H. Save, S. Bettadpur, B. Tapley (2012)
Reducing errors in the GRACE gravity solutions using regularization
Journal of Geodesy, 86
J. Awange, K. Fleming, M. Kuhn, W. Featherstone, B. Heck, I. Anjasmara (2011)
On the suitability of the 4° × 4° GRACE mascon solutions for remote sensing Australian hydrology
Remote Sensing of Environment, 115
K. Seo, C. Wilson, Jianli Chen, D. Waliser (2006)
GRACE's spatial aliasing error
Geophysical Journal International, 172
P. Muller, W. Sjogren (1968)
Mascons: Lunar Mass Concentrations
Science, 161
B. Tapley, S. Bettadpur, J. Ries, P. Thompson, M. Watkins (2004)
GRACE Measurements of Mass Variability in the Earth System
Science, 305
H. Pritchard, S. Luthcke, A. Fleming (2010)
Understanding ice-sheet mass balance: progress in satellite altimetry and gravimetry
Journal of Glaciology, 56
B. Tapley, S. Bettadpur, M. Watkins, C. Reigber (2004)
The gravity recovery and climate experiment: Mission overview and early results
Geophysical Research Letters, 31
D. Rowlands, S. Luthcke, S. Klosko, F. LeMoine, D. Chinn, J. Mccarthy, C. Cox, O. Anderson (2005)
Resolving mass flux at high spatial and temporal resolution using GRACE intersatellite measurements
Geophysical Research Letters, 32
P. Moore, H. Boomkamp, S. Carnochan, R. Walmsley (1999)
FAUST: Multi-satellite orbital dynamics software
Advances in Space Research, 23
(2002)
A reconciled estimate of icesheet mass balance
W. Farrell (1972)
Deformation of the Earth by surface loads
Reviews of Geophysics, 10
S. Luthcke, T. Sabaka, B. Loomis, A. Arendt, J. Mccarthy, J. Camp (2013)
Antarctica, Greenland and Gulf of Alaska land-ice evolution from an iterated GRACE global mascon solution
Journal of Glaciology, 59
R. Klees, E. Revtova, B. Gunter, P. Ditmar, E. Oudman, H. Winsemius, H. Savenije (2008)
The design of an optimal filter for monthly GRACE gravity models
Geophysical Journal International, 175
B. Wouters, E. Schrama (2007)
Improved accuracy of GRACE gravity solutions through empirical orthogonal function filtering of spherical harmonics
Geophysical Research Letters, 34
S. Luthcke, A. Arendt, D. Rowlands, J. Mccarthy, C. Larsen (2008)
Recent glacier mass changes in the Gulf of Alaska region from GRACE mascon solutions
Journal of Glaciology, 54
T. Sabaka, D. Rowlands, S. Luthcke, J. Boy (2010)
Improving global mass flux solutions from Gravity Recovery and Climate Experiment (GRACE) through forward modeling and continuous time correlation
Journal of Geophysical Research, 115
Awange (2011)
On the suitability of the 4 degree x 4 degree GRACE mascon solutions for remote sensing Australian hydrology
Rem. Sens. Environ., 115
I. Velicogna, J. Wahr (2006)
Measurements of Time-Variable Gravity Show Mass Loss in Antarctica
Science, 311
S. Luthcke, D. Rowlands, F. Lemoine, S. Klosko, D. Chinn, J. Mccarthy (2006)
Monthly spherical harmonic gravity field solutions determined from GRACE inter‐satellite range‐rate data alone
Geophysical Research Letters, 33
Jianli Chen, C. Wilson, B. Tapley, D. Blankenship, D. Young (2008)
Antarctic regional ice loss rates from GRACE
Earth and Planetary Science Letters, 266
A. Shepherd, E. Ivins, G. A, V. Barletta, M. Bentley, S. Bettadpur, K. Briggs, D. Bromwich, R. Forsberg, N. Galin, M. Horwath, S. Jacobs, I. Joughin, Matt King, J. Lenaerts, Jilu Li, S. Ligtenberg, A. Luckman, S. Luthcke, M. McMillan, Rakia Meister, G. Milne, J. Mouginot, A. Muir, J. Nicolas, J. Paden, A. Payne, H. Pritchard, E. Rignot, H. Rott, L. Sørensen, T. Scambos, B. Scheuchl, E. Schrama, Ben Smith, A. Sundal, J. Angelen, W. Berg, M. Broeke, D. Vaughan, I. Velicogna, J. Wahr, P. Whitehouse, D. Wingham, D. Yi, D. Young, H. Zwally (2012)
A Reconciled Estimate of Ice-Sheet Mass Balance
Science, 338
R. Riva, B. Gunter, T. Urban, B. Vermeersen, R. Lindenbergh, M. Helsen, J. Bamber, R. Wal, M. Broeke, B. Schutz (2009)
Glacial Isostatic Adjustment over Antarctica from combined ICESat and GRACE satellite data
Earth and Planetary Science Letters, 288
J. Klokočník, R. Gooding, C. Wagner, J. Kostelecký, A. Bezdek (2012)
The Use of Resonant Orbits in Satellite Geodesy: A Review
Surveys in Geophysics, 34
C. Dunn, W. Bertiger, Y. Bar-Sever, S. Desai, B. Haines, D. Kuang, G. Franklin, I. Harris, G. Kruizinga, T. Meehan, S. Nandi, D. Nguyen, T. Rogstad, J. Thomas, J. Tien, Larry Romans, M. Watkins, Sien-Chong Wu, S. Bettadpur, Jeongrae Kim (2003)
Instrument of Grace: GPS augments gravity measurements
S. Bruinsma, J. Lemoine, R. Biancale, N. Val’es (2010)
CNES/GRGS 10-day gravity field models (release 2) and their evaluation
Advances in Space Research, 45
S. Swenson, J. Wahr (2002)
Methods for inferring regional surface‐mass anomalies from Gravity Recovery and Climate Experiment (GRACE) measurements of time‐variable gravity
Journal of Geophysical Research, 107
D. Rowlands, S. Luthcke, J. Mccarthy, S. Klosko, D. Chinn, F. LeMoine, J. Boy, T. Sabaka (2010)
Global Mass Flux Solutions from GRACE: A Comparison of Parameter Estimation Strategies - Mass Concentrations Versus Stokes Coefficients
Journal of Geophysical Research, 115
J. Kusche (2007)
Approximate decorrelation and non-isotropic smoothing of time-variable GRACE-type gravity field models
Journal of Geodesy, 81
J. Wahr, S. Swenson, V. Zlotnicki, I. Velicogna (2004)
Time‐variable gravity from GRACE: First results
Geophysical Research Letters, 31
J. Wahr, Mery Molenaar, F. Bryan (1998)
Time variability of the Earth's gravity field: Hydrological and oceanic effects and their possible detection using GRACE
Journal of Geophysical Research, 103
T. Jacob, J. Wahr, W. Pfeffer, S. Swenson (2012)
Recent contributions of glaciers and ice caps to sea level rise
Nature, 482
Shin‐Chan Han, C. Shum, C. Jekeli, C. Kuo, C. Wilson, K. Seo (2005)
Non-isotropic filtering of GRACE temporal gravity for geophysical signal enhancement
Geophysical Journal International, 163
(2007)
Evaluation and validation of mascon recovery using GRACE KBRRdatawith independentmass flux estimates in theMississippi Basin
F. LeMoine, S. Luthcke, D. Rowlands, D. Chinn, S. Klosko, C. Cox (2007)
The use of mascons to resolve time-variable gravity from GRACE
S. Klosko, D. Rowlands, S. Luthcke, F. LeMoine, D. Chinn, M. Rodell (2009)
Evaluation and validation of mascon recovery using GRACE KBRR data with independent mass flux estimates in the Mississippi Basin
Journal of Geodesy, 83

Publisher: Oxford University Press
Copyright: Copyright © 2022 The Royal Astronomical Society
ISSN: 0956-540X
eISSN: 1365-246X
DOI: 10.1093/gji/ggu402
Publisher site: See Article on Publisher Site

Abstract

Downloaded from https://academic.oup.com/gji/article/200/1/503/752750 by DeepDyve user on 20 July 2022 Geophysical Journal International Geophys. J. Int. (2015) 200, 503–518 doi: 10.1093/gji/ggu402 GJI Gravity, geodesy and tides Mass change from GRACE: a simulated comparison of Level-1B analysis techniques 1 1 1,2 Stuart B. Andrews, Philip Moore and Matt. A. King School of Civil Engineering and Geosciences, Newcastle University, Newcastle upon Tyne, NE1 7RU, United Kingdom. E-mail: [email protected] School of Land and Food, University of Tasmania, Private Bag 76, Hobart 7001, Australia Accepted 2014 October 13. Received 2014 October 10; in original form 2014 March 31 SUMMARY Spherical harmonic and mascon parameters have both been successfully applied in the recov- ery of time-varying gravity ﬁelds from Gravity Recovery and Climate Experiment (GRACE). However, direct comparison of any mass ﬂux is difﬁcult with solutions generated by differ- ent groups using different codes and algorithms. It is therefore opportune to compare these methodologies, within a common software base, to understand potential limitations associated with each technique. Here we use simulations to recover a known monthly surface mass distri- bution from GRACE KBRR data. The ability of spherical harmonic and mascon parameters to resolve basin-level mass change is quantiﬁed with an assessment of how the noise and errors, inherent in GRACE solutions, are handled. Recovery of a noise and error free GLDAS anomaly revealed no quantiﬁable difference between spherical harmonic and mascon parameters. Ex- pansion of the GLDAS anomaly to degree and order 120 shows that both spherical harmonic and mascon parameters are affected by comparable omission errors. However, the inclusion of realistic KBRR noise and errors in the simulations reveals the advantage of the mascon parameters over spherical harmonics at reducing noise and errors in the higher degree and order harmonics with an rms (cm of EWH) to the GLDAS anomaly of 10.0 for the spherical ◦ ◦ harmonic solution and 8.8 (8.6) for the 4 (2 ) mascon solutions. The introduction of a con- straint matrix in the mascon solution based on parameters that share geophysical similarities is shown to further reduce the signal lost at all degrees. The recovery of a simulated Antarctic mass loss signal shows that the mascon methodology is superior to spherical harmonics for this region with an rms (cm of EWH) of 8.7 for the 2 mascon solution compared to 10.0 for the spherical harmonic solution. Investigating the noise and errors for a month when the satellites were in resonance revealed both the spherical harmonic and mascon methodologies are able to recover the GLDAS and Antarctic mass loss signal with either a comparable (spherical harmonic) or improved (mascon) rms compared to non-resonance periods. Key words: Inverse theory; Satellite geodesy; Gravity anomalies and Earth structure; Time variable gravity; Global change from geodesy; Antarctica. is measured using a K/Ka band instrument which has an accu- 1 INTRODUCTION racy of a few microns (Dunn et al. 2003). The ﬁrst time-derivative The Gravity Recovery and Climate Experiment (GRACE) twin provides a K/Ka band range rate (KBRR) measurement, which is satellite mission (Tapley et al. 2004b) was launched in 2002 with generally used during GRACE data processing and has a precision –1 the aim of mapping the Earth’s time-varying gravity (TVG) ﬁeld of 0.1 µms (Luthcke et al. 2006b). The GRACE satellites are every ∼30 days (Wahr et al. 2004) at spatial scales of hundreds of equipped with Global Positioning System (GPS) receivers to posi- kilometres. GRACE consists of two near-polar satellites separated tion the satellites and on-board accelerometers to measure the non- along-track by ∼220 km. This separation results in the two satel- conservative forces that affect the satellites allowing, shortly after lites experiencing small orbital perturbations, caused by local mass launch, the detection of changes in mass corresponding to 1.4 cm anomalies, at different times. It is through measuring the chang- of equivalent water height (EWH) over a radius of 750 km every 30 ing distance between the satellites, caused by these perturbations, days (Wahr et al. 2004). Analysis improvements have since driven that GRACE allows us to monitor spatial and temporal geophysical this uncertainty/radius downward. The data from these instruments mass variations in the gravity ﬁeld. The range between the satellites are available to the scientiﬁc community as Level-1B data, namely The Authors 2014. Published by Oxford University Press on behalf of The Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. Downloaded from https://academic.oup.com/gji/article/200/1/503/752750 by DeepDyve user on 20 July 2022 504 S. B. Andrews, P. Moore and M. A. King the raw satellite products after a low-pass ﬁlter and time correction matrices allows users to constrain the higher degrees of the TVG (Chen et al. 2008). ﬁelds to the mean ﬁeld (e.g. Bruinsma et al. 2010;Save et al. The ability to accurately recover the Earth’s temporal mass 2012) but these are not typically provided as part of the Level-2 anomalies every 30 days to a few hundred kilometres allows the solutions and are only available to users generating solutions from monitoring of global geophysical processes between successive 30 Level-1B data. day solutions and provides a means to assess and improve geo- Mass concentration (mascon) parameters, estimated using only physical models. GRACE data has been successfully used for es- GRACE Level-1B KBRR data, have also been applied to the mod- timating ice mass balance trends (e.g. King et al. 2012; Shepherd elling of TVG, primarily at the NASA Goddard Space Flight Centre et al. 2012; Luthcke et al. 2013), hydrological change (e.g. Han (GSFC; Rowlands et al. 2005, 2010; Luthcke et al. 2006b, 2008; et al. 2009; Awange et al. 2011; Ramillien et al. 2011), ocean mass Sabaka et al. 2010). The estimation of mascon parameters, ﬁrst variability (e.g. Chambers et al. 2010; Chambers & Bonin 2012), applied by Muller & Sjogren (1968) to model lunar gravity vari- Glacial Isostatic Adjustment (GIA; e.g. Riva et al. 2009) and study- ations, has allowed the Earth’s TVG to be modelled at a spatial ing Earthquake displacements (Han et al. 2008). Without the spatial and temporal resolution unachievable using conventional spherical and temporal resolution provided by GRACE, the separation of the harmonics due to the embedded smoothing resulting in less sig- signals from many of these geophysical phenomena would not have nal loss (Rowlands et al. 2010). The mascon methodology beneﬁts been possible (Swenson & Wahr 2002). from constraint equations applied during the least squares inversion Since the launch of GRACE several processing methodologies which removes the requirement of post-processing, making possible have been applied to the GRACE Level-1B data with the aim of es- solutions with a 10 day temporal resolution and a spatial sampling ◦ ◦ timating the temporal variations of mass anomalies while reducing of 4 or 2 (Rowlands et al. 2010) with the most recent published the impact that noise and errors have on the solutions. The pres- mascon solutions using a 1 spatial sampling (Luthcke et al. 2013). ence of north–south ‘stripes’ in the solution results from the limited The mascon constraint matrix is anisotropic in nature as it can be longitudinal sampling of GRACE, which reduced observability of applied to constrain distinct geographical regions while also taking the sectorial and near-sectorial coefﬁcients (Klees et al. 2008), and into account the full noise covariance matrix (Luthcke et al. 2013). result from the orbital conﬁguration of GRACE; the KBRR data The use of mascon parameters has allowed basin scale interannual are only collected along the line of the orbit. Stripes also suggest a variability to be extracted from the GRACE time-series (Rowlands high degree of correlation in the GRACE errors (Swenson & Wahr et al. 2010) due to the increased spatial sampling and temporal res- 2006). If not correctly dealt with, measurement and processing er- olution of the solutions. Mascon TVG estimates, unlike spherical rors and unmodelled mass variations in the GRACE solution could harmonic coefﬁcients, can be either local or global in extent as the be interpreted as real mass change (Wahr et al. 2006). Random KBRR data associated with the mascon is dominated by mass ﬂux and systematic errors, present in the instrument data and back- below the satellite (Luthcke et al. 2008). Mascon parameters can be ground models, can dominate the solution at short wavelengths and estimated from short arcs of KBRR data through the estimation of increase rapidly with the spherical harmonic degree (Klees et al. short arc baseline parameters (Rowlands et al. 2002). However this 2008) while errors in geophysical models can cause poor estimates requires GRACE accelerometer data to deal with non-conservative of certain spherical harmonic coefﬁcients and alias into the solu- forces and the satellite orbits to have been previously well deter- tions (Seo et al. 2008). In addition, the ground coverage of the mined using GNV1B (GPS Navigation Data Level-1B) positioning satellites also affects the quality of the solutions. Periods where the data. Estimating mascon parameters from Level-1B data is differ- GRACE satellites are in near 3, 4 and 7 day repeat orbits modify ent from the mascons used in Jacob et al. (2012) and Schrama the noise characteristics of the GRACE solutions (Save et al. 2012). et al. (2014), where mascon parameters are calculated from Level-2 Noise and errors in the GRACE solutions need to be appropriately spherical harmonics data. treated and mitigated before any meaningful mass anomalies can be Both spherical harmonic and mascon methodologies are applied estimated. to the recovery of TVG, but most users of GRACE data are heav- ily reliant on pre-processed Level-2 spherical harmonic products. The most common methodology used to estimate the temporal variations of mass anomalies is the recovery of spherical harmonic It is therefore important that an assessment of these methodologies (Stokes’) coefﬁcients (Tapley et al. 2004b;Wahr et al. 2004). The is undertaken to provide users with an understanding of their as- GRACE data processing centres, namely Centre for Space Re- sociated limitations and to assess the ability of each technique to search (CSR), Jet Propulsion Laboratory (JPL) and GeoForschungs resolve basin-level mass changes at a variety of spatial scales. The Zentrum (GFZ), each provide monthly estimates of the global TVG need for this assessment is further compounded as these processing ﬁeld in the form of spherical harmonic ﬁelds to degree and order 60– methodologies are applied by different groups using different codes 90. These Level-2 products are normally computed directly from and algorithms, making it hard to directly compare any subsequent Level-1B KBRR and GPS data, but Level-2 solutions have been suc- mass ﬂux analysis. Pritchard et al. (2011) highlighted that differ- cessfully computed using GRACE KBRR data only (Luthcke et al. ences within Level-1B processing applied at different groups is the 2006a). Due to the complex nature and computational requirement main source of the discrepancies in TVG solutions, with solutions of Level-1B data processing, Level-2 products are widely used by generated using the same fundamental processing software agreeing researchers who undertake studies of mass ﬂux analysis. A num- within estimated errors (see also Shepherd et al. 2012). ber of post-processing strategies have been proposed to remove While there are a limited number of studies which have compared the noise and errors in the solutions including the use of Gaus- different approaches within the same software base (e.g. Rowlands sian smoothing (e.g. Wahr et al. 1998, 2004; Tapley et al. 2004a), et al. 2010; Awange et al. 2011), no study, as far as we are aware, ﬁltering using empirical orthogonal functions (EOF; e.g. Wouters has been undertaken that simulates the recovery of a known sig- &Schrama 2007), destriping (e.g. Swenson & Wahr 2006)and nal using both mascon and spherical harmonic solutions within forward modelling (e.g. Wouters et al. 2008; Schrama & Wouters the same software. Accordingly, we undertake a comparison of 2011). Gaussian smoothing is probably the most common process- solutions generated through the estimation of mascon parameters ing strategy (Rowlands et al. 2010). Access to solution normal and spherical harmonic coefﬁcients. Simulations will provide an Downloaded from https://academic.oup.com/gji/article/200/1/503/752750 by DeepDyve user on 20 July 2022 Mass change from GRACE: simulated comparison 505 accuracy assessment and quantify the capability of each technique to resolve basin-level mass changes at a variety of spatial scales while understanding how the methodologies handle the noise and errors inherent at higher degree and orders. Actual noise and errors, obtained during GRACE Level-1B data processing, will be used to create simulated stripes for use during the simulations. While we only simulate the recovery of a single signal, the assessment of the methodologies aims to provide users with an understanding of their associated limitations, as correctly recovering individual monthly solutions is imperative when using these solutions for de- riving longer term trends. We will use simulations to compare a mascon and spherical har- monic solution using a standard Gaussian smoothing ﬁlter for the post-processing of the spherical harmonic solution and a spatial con- straint matrix applied to the mascon solution. The aim of this paper is to investigate the effect of the fundamental processing methodol- ogy applied to the GRACE Level-1B data on the recovered monthly time variable gravity ﬁeld. It is worth noting that different post- processing strategies applied to the spherical harmonic solutions will generate different mass ﬂux estimates, however only the dif- ference in the fundamental Level-1B processing strategy will be considered in this work. The outline of the paper is as follows. In Section 2 we describe the Figure 1. Distribution of 2 equal area mascons around Antarctica. Antarc- formation of mascon parameters before discussing the processing tic drainage basins (Zwally et al. 2012) are also shown (red). methodology in Section 3. The results of the simulations will be presented in Section 4. Throughout this paper special attention will Using eq. (1) each mascon in the grid can be expanded into be paid to the formation and behaviour of the mascon parameters to spherical harmonics as aid understanding of the methodology. Spherical harmonic solutions l l max (m) (m) are widely dealt with in the literature and their estimation is not m ˜ ˜ δ (φ, λ) = R ρ C cos mλ + S sin mλ e w lm lm explicitly described here. l=0 m=0 ( ) × P cos φ , (2) lm (m) (m) ˜ ˜ where C and S are the dimensionless normalized surface 2 FORMATION OF THE MASCON lm lm harmonics and P normalized associated Legendre polynomials. PARAMETERS lm The sum over degree l and order m in eq. (2) is truncated tol = 60. max The mascon parameters are based on the GSFC methodology which We utilize eq. (2) to modify eqs (11) and (12) of Wahr et al. (1998)to (m) involves calculating a scaling factor for a set of ‘delta’ spherical obtain (delta) spherical harmonic gravitational coefﬁcients, C lm harmonic coefﬁcients that represent a uniform layer of surface mass (m) and S , that describe the required mass spread uniformly over a lm over an area which can be added to the mean background ﬁeld and global set of mascons to give used to represent the mass ﬂux at a certain time. Rowlands et al. 3p 1 + k w l (m) (m) (2010) showed that the only difference between these differential C = H(t) C lm lm 4π R ρ 2l + 1 coefﬁcients is the shape, size and location of the area. Each separate e av area is a mascon as deﬁned by estimating a scale factor H(t) at time 3p 1 + k (m) w l (m) S = H(t) S . (3)) t for a uniform mass of one unit of EWH over the mascon area, lm lm 4π R ρ 2l + 1 e av namely 5517 kg ⎧ ⎫ In eq. (3), ρ is the average density of the Earth ( )and k av 3 l ⎨ 1, over mascon⎬ the load Love number of degree l (Farrell 1972) to account for the δ (φ, λ) = R ρ , (1) deformation of the load on an elastic Earth. The mascon parameter 0, elsewhere ⎩ ⎭ (H) is a scale factor of the spherical harmonic coefﬁcients used to describe that mascon (Rowlands et al. 2010). During the orbital where the mascon is centred on geodetic latitude φ, longitude λ; integration we use the spherical harmonic coefﬁcients from eq. (3) R is the radius of the Earth and ρ the density of fresh water e w –3 to calculate H by least-squares utilizing the partial derivatives of the (1000 kg m ). KBRR observation data (cf. Luthcke et al. 2013), that is In this paper, we use global equal area masons, where each mas- ◦ ◦ con is rectangular in latitude and longitude spanning 2 or 4 in l max l ∂KBRR ∂KBRR ∂C I I lm latitude with the longitude span selected to give approximately the A = = IJ ∂ H ∂C ∂ H ◦ ◦ ◦ ◦ J J lm same area as a 2 by 2 or 4 by 4 mascon located at the equator. l=1 m=0 ◦ ◦ ◦ ◦ While not strictly 2 by 2 or 4 by 4 , for simplicity we refer to ∂KBRR ∂S lm them as such hereafter. The only exception is at the poles where + , (4) ∂S ∂ H ◦ J lm polar caps of 2 in latitude are used. Due to the convergence of longitude at the poles, the number of mascons in a given latitude where A are elements of the design matrix, ∂KBRR /∂ H is IJ I J band will reduce poleward as seen in Fig. 1 for a grid of 2 mascons the partial derivative of the KBRR observation i with respect to the J J over Antarctica. mascon parameter j, ∂KBRR /∂C and ∂KBRR /∂S are I I lm lm Downloaded from https://academic.oup.com/gji/article/200/1/503/752750 by DeepDyve user on 20 July 2022 506 S. B. Andrews, P. Moore and M. A. King the partial derivatives of the KBRR observations with respect to the and S = 0.0001 were tested but were found to over/under constrain J J delta gravitational coefﬁcients while ∂C /∂ H and ∂S /∂ H the simulated solutions. However, the choice of S is dependent on J J lm lm are the partial derivatives of the delta gravitational coefﬁcients with the values chosen for other parameters in eqs (8) and (9) as the respect to the mascon parameter j.Let xˆ be the least squares solution weight matrix V is affected by the values assigned to S, D and T in to the KBRR linearized observation equations Ax − b − v = 0 eqs (8) and (9). The same scale factor is used for all constraints in with observation weight matrix W and Gaussian error v . The solu- the V matrix. tion vector xˆ includes GRACE orbital parameters and either mascon The need for temporal smoothing increases as the number of days parameters, H, or corrections H to mascons parameters H at the included in the mascon solution increases. For consistency with the th k stage of an iterative procedure (cf. Sabaka et al. 2010). The vector spherical harmonic approach, no temporal constraint is applied to b is formed from the KBRR calculated minus observed residuals. the mascon monthly solutions presented here, leading to a modiﬁed Analogously, we incorporate spatial and temporal constraints on version of the constraint matrix, namely the mascons through pseudo-observations Cx − q − v = 0 with 2 d IJ 1− V = S × e . (9) weight matrix V,where q is a constant vector incorporating both IJ non-zero constraints and the component −CH if the procedure is In our solutions, constraints are applied between all mascon pairs iterative. Then, xˆ is obtained from of the same type (land or ocean), but no constraint is applied between T T T T a pair consisting of a land and an ocean mascon. The constraint (A WA + C VC)xˆ = A Wb + C Vq. (5) matrix can be easily modiﬁed to allow different constraints between The matrix, C, incorporates the mascon spatial and/or temporal pairs of mascon, as we will show in Sections 4.6 and 4.7. constraints. For mascons, J , J , with the same geophysical property 1 2 the pseudo-observations can be written as 3 SIMULATED DATA PROCESSING H − H = 0, (6) J J 1 2 STRATEGY where J = J ; J , J ∈ A ,the kth mascon region. Additional 1 2 1 2 k During simulations we aim to recover a known surface mass dis- pseudo-observations are applied to ensure total mass is conserved tribution from GRACE KBRR data using Newcastle University’s (equivalent to a degree 0 term) and that the degree 1 harmonics orbit determination software Faust (Moore et al. 1999). We compute for surface mass are zero; GRACE is insensitive to degree 1 har- the effect of the surface mass distribution through a spherical har- monics as the satellites orbit around the instantaneous centre of monic expansion to a speciﬁed degree and order, with the resultant mass. The inclusion of background models during the orbital in- gravitational ﬁeld added to the mean GIF48 ﬁeld (Ries et al. 2011), tegration could result in non-zero degree 1 terms and mass not from which simulated noise and error free KBRR tracking data are being conserved in the solution vector without the use of pseudo- computed. The simulations then attempt to recover the known sig- observations; individual models do not conserve mass and may nal from the tracking data up to degree and order 60. The normal give non-zero degree 1 coefﬁcients (Wahr et al. 1998). The degree equations from each arc, here taken as 60 min, were combined to 0 and 1 pseudo-observations are prescribed a high weight in V to create monthly solutions for the two methodologies. We only use ensure the solution essentially satisﬁes these conditions. The four KBRR data in our solutions as Luthcke et al. (2006a) showed that constraint equations for degrees l = 0, 1 and order m = 0, 1 are of the gravity ﬁeld can be successfully recovered from KBRR data us- the form ing spherical harmonic parameters while estimating the Rowlands (J) (atmos) (ocean) H C =−C − C , (7) J et al. (2002) short arc baseline parameters. Mascon parameters are l,m l,m l,m usually estimated using only KBRR data. The use of KBRR data (atmos) (ocean) only removes the computational requirement of processing the GPS where the summation is over all mascons and C and C l,m l,m data which can weaken and corrupt the KBRR measurement if not are the harmonics from the GRACE dealiasing product for the ocean handled correctly (Luthcke et al. 2006a). and atmosphere. By stipulating eq. (7), the total contribution due Within our simulations, the only difference between the two pro- to the mascons for the degree 0 and 1 harmonics will offset the a cessing strategies is the gravity ﬁeld parameters being estimated priori mass of the ocean and atmosphere if used in the GRACE (spherical harmonic or mascon). Thus, the underlying processing is computations, giving zero values for these harmonics at all times. the same for the two methodologies to the extent that any differences The weight matrix V utilizes constraints between the mascon in the solution will be directly due to the parametrization (Rowlands parameters, et al. 2010) and how noise and errors propagate into the solutions. d t IJ IJ 2− − Mathematically the two methodologies are equivalent (Sabaka et al. D T V = S × e , (8) IJ 2010). where d is the distance in kilometres between the mascons, D The choice of arc length was based on experimentation in deriving IJ the correlation distance, t the time difference in days between gravity ﬁeld solutions using KBRR (and accelerometer) data with IJ the solutions, T the correlation time and S a scale factor. The scale the satellite state vectors estimated simultaneously using the short factor is required to ensure that the constraints are in proportion baseline parameters. No GPS/GNV1B data or empirical parameters to the normal matrices. Setting the constraint too small will not are required to stabilize the solutions. Analysis of 4 months of substantially reduce the noise and error in the solution, while adding solutions revealed that the average ﬁt of the KBRR residuals was −1 −1 a large constraint will oversmooth the signal and remove both noise 0.32 µms using 60 min arcs and 0.31 and 0.41 µms for 30 and error and real signal (Klosko et al. 2009). The value used for min and 90 min arcs, respectively. The 60 min arcs were favoured S, to be discussed later, depends on the units used for the normal over the 30 min arcs due to the reduction in storage requirement. matrices and constraint matrix. In Section 4.5, we found S = 0.001, 24 hr arcs were used to calculate the GRACE accelerometer bias the value used by Lemoine et al. (2007) and Rowlands et al. (2010), values before recovery of the gravity ﬁeld parameters. Alongside to be the optimal scale factor for the simulations. Values of S = 0.01 the accelerometer biases, these runs estimate the full state vector Downloaded from https://academic.oup.com/gji/article/200/1/503/752750 by DeepDyve user on 20 July 2022 Mass change from GRACE: simulated comparison 507 from KBRR, GNV1B and accelerometer data while using empirical ogy signal will reduce the leakage in the GRACE signal as seen parameters to absorb unmodelled KBRR measurement variations in Sabaka et al. (2010). However, as this assumes that sufﬁciently for quality control purposes. accurate models are available, forward modelling of hydrology will The results presented in Section 4.2 for a noise and error free not be considered in this paper. The simulation was repeated with simulation show that we are able to create KBRR tracking data the source area described to spherical harmonic degree and order that describes a simulated mass. Noise and errors are included in 120 (not shown). The results reveal that the additional degree and the simulation in Section 4.5, through creation of simulated stripes, order 61–120 has minimal effect on KBRR residuals conﬁrming the allowing us to examine how the solutions differ at higher degree result seen in Rowlands et al. (2010). and orders using the different methodologies. 4.2 Simulated recovery of the GLDAS anomaly 4 SIMULATIONS The next simulation was the recovery of a known gravity signal ◦ ◦ through estimation of spherical harmonic and 2 and 4 mascon pa- 4.1 Effect of a simulated mass on observed KBRR rameters to assess the performance of the different methodologies residuals in an ideal case. The introduced signal was the June 2006 Global Land Data Assimilation System (GLDAS) Version 1 (Rodell et al. To understand the performance of the two methodologies a number 2004) anomaly to the 2003–2010 mean derived from monthly av- of experiments were undertaken. The ﬁrst experiment examined the eraged data with a 1 spatial resolution. Since hydrology is the effect of a simulated mass on the KBRR data. The simulated mass largest seasonal contribution to the GRACE signal we are sim- was created by loading a 4 block at the equator, described from ulating the recovery of the main spatial and temporal variation, spherical harmonic degree and order 2–60, with 30 cm of EWH. The with similar characteristics to that contained within real GRACE spherical harmonic expansion was added to the mean background data. We scaled the GLDAS anomaly upwards so the degree vari- ﬁeld to create error free KBRR tracking data. The effect of the mass ances contained a similar power as the June 2006 CSR release on the KBRR residuals is plotted in Fig. 2, with t = 0 at the centre of 5 Level-2 ﬁeld (Bettadpur 2012), as not all geophysical signal is the block. The simulation was repeated with the source mass loaded present in the GLDAS anomaly. The GLDAS anomaly was derived with 10 cm of EWH. The boundary of the 4 block is plotted giving to degree and order 60 and 120 (Fig. 3) with tracking data created a spatial domain ∼440 km in diameter. from the different harmonic expansions. As the degree variances Fig. 2 shows that both simulated masses modify the KBRR resid- are identical up to degree and order 60, recovery of the GLDAS uals inside and outside the source area as found by Rowlands et al. anomaly described to degree and order 120 allows us to observe how (2010). The largest KBRR values are observed at the boundary parameter estimation is affected by omission errors. Gravity ﬁeld re- of the source area with a change from positive to negative as the covery was undertaken, for both GLDAS anomaly models, through GRACE satellite pair move over the centre. Either side of the source area is an extensive tail revealing that GRACE is sensitive to mass over a large area. Spatially, Fig. 2 covers an area of ∼2500 km. The amplitude and extent of the tail is a function of the size of the mass anomaly. Seo et al. (2006) found that the distribution of leakage is concentrated in regions with high hydrological variability. As the hydrology component is the largest seasonal signal in the GRACE data, Fig. 2 shows mass signal leakage from the source area into surrounding areas. These results, and the ﬁndings of Seo et al. (2006), lead to the conclusion that forward modelling of the hydrol- Figure 3. GLDAS mass anomaly for June 2006 as differenced from the Figure 2. The residuals resulting from the 4 block loaded with 30 cm of 2003 to 2010 mean to degree and order 60 (a); degree and order 120 (b). EWH (green) and 10 cm of EWH (blue). The boundary of the simulated Note that there is no GLDAS signal in Antarctica or Greenland or over the mass (red) is plotted for reference. oceans. Units of cm of EWH. Downloaded from https://academic.oup.com/gji/article/200/1/503/752750 by DeepDyve user on 20 July 2022 508 S. B. Andrews, P. Moore and M. A. King Figure 4. Rms of the differences between the calculated gravity ﬁeld solu- Figure 5. Degree variances of the spherical harmonic and mascon solutions. tions and the noise and error free GLDAS anomaly in cm of EWH. GLDAS The degree variances of the noise and error free input GLDAS anomaly are described to degree and order 60 (a); and degree and order 120 (b). The included for comparative purposes. distance refers to the Gaussian smoothing radius (spherical harmonic coef- ﬁcient solution) or the correlation distance (mascon solutions). The distance GLDAS anomaly is 0.05 cm of EWH, which is the same as the used in the mascons solutions (triangles) is not directly comparable to that spherical harmonic recovery with 50 km Gaussian smoothing. In- used in the harmonic solution (crosses) and do not equate to the same creasing the correlation distance in the constraint matrix increases smoothing radius. the effect of the constraint matrix and over constrains the mascon solution, reducing high frequency signal. spherical harmonic estimation to degree and order 60 (3717 pa- Fig. 4(b) shows the rms of the recovered anomaly to the input ◦ ◦ rameters) while the 4 (2 ) equal-area mascon solutions involved GLDAS data, described to degree and order 120. This allows us to estimation of 2564 (10 292) mascon parameters. The monthly so- compare the effect of omission errors on the different parameters. lution was collated from the individual 60 min arcs. No noise or The pattern is similar to the ideal case (GLDAS to degree and order errors were added to the solution at this stage. 60), but with an increase in all rms values. Methodologies that are Fig. 4(a) shows the rms in cm of EWH of the recovered anomaly to essentially truncated to degree and order 60 are, of course, unable the input GLDAS data, described to degree and order 60. This is the to fully recover an input signal described to degree and order 120 ideal case. The rms is provided for various smoothing distances. For and hence unable to capture the additional short wavelength (high spherical harmonic coefﬁcients the distance refers to the Gaussian degree) information. smoothing radius, while for the mascon solutions the distance is D Fig. 5 shows the degree variances of the spherical harmonic and of eq. (9). As such, the two distances are not directly comparable mascon solutions to degree 60. The degree variance of the input as they do not equate to the same smoothing radius. The distance GLDAS anomaly is included for reference. All three solutions match 0 km indicates that no Gaussian smoothing has been applied or that the degree variance of the input GLDAS anomaly to near degree the mascon solutions are unconstrained. 40 from where the 4 mascon solution loses power. This reduction Fig. 4(a) reveals that, without Gaussian smoothing, the spherical is the result of the under sampling of the input signal by the 4 harmonic solution is able to exactly recover the reference degree mascon solution. The spherical harmonic coefﬁcient and 2 mascon and order 60 GLDAS anomaly using perfect data. Adding Gaussian solutions are able to replicate the degree variances of the GLDAS smoothing dampens the amplitude of the recovered signal. For the anomaly. The slight reduction in power at degree ∼55–60 in the 4 mascon solution, the smallest rms is obtained when a correlation 2 mascon solution is the result of the 50 km correlation distance distance of 150 km is used in the constraint matrix, but the solution used in the constraint matrix to stabilize the inversion. The ability is not able to fully recover the input GLDAS anomaly. This is to be of the spherical harmonic and 2 mascon solutions to replicate expected as the GLDAS anomaly was expanded to spherical har- the degree variances of the GLDAS anomaly up to degree 60 in monics degree and order 60 (3717 gravity ﬁeld parameters) while Fig. 5 shows that we are able to reproduce a known input signal the 4 mascon solution involves estimation of 2564 mascon param- from the simulated KBRR tracking data, validating the processing eters. The recovery favours the spherical harmonic solution over the methodology. For the remainder of the paper, all processing uses 4 mascon solution with the latter unable to recover the gravity ﬁeld. the GLDAS anomaly described to degree and order 120 and hence The use of the constraint matrix in the mascon solution does enable allows consideration of omission errors. a small improvement in the ability to recover this short wavelength signal. The 2 mascon solution involves estimation of 10 292 mascon 4.3 Mascon centre-to-centre distance parameters and consequently nearly a factor 4 more parameters than the spherical harmonic solution. If no constraint matrix is As part of the evaluation of the force model the contribution of applied, the solution does not invert as it is rank deﬁcient by a every mascon to the accelerations of the satellites is calculated. large margin (Rowlands et al. 2010). Adding a constraint matrix Thus, as for spherical harmonics, mascon solutions presented so with ≥50 km correlation distance stabilizes the solution and allows far are global in extent with the drawback of added computational recovery of the GLDAS anomaly. The resulting rms to the input expense, at least in the case of the 2 mascon solutions. Fig. 2 Downloaded from https://academic.oup.com/gji/article/200/1/503/752750 by DeepDyve user on 20 July 2022 Mass change from GRACE: simulated comparison 509 Figure 7. Degree variances of the 2 mascon solutions for a range of Figure 6. Degree variances of the 4 mascon solutions for a range of centre-to-centre cut-off distances. The degree variances of the input GLDAS centre-to-centre cut-off distances. The degree variances of the input GLDAS anomaly (solid green) and original solution (solid black) are included for anomaly (solid green) and original solution (solid black) are included for reference. All solutions here are noise and error free. reference. All solutions here are noise and error free. revealed that a mass anomaly leaks outside of the source area, and by extension neighbouring mass anomalies will leak into the source area. Fig. 2 also showed that this leakage is not global, but spatially limited. The limited spatial extent of a mass anomaly on gravitation experienced at the satellite altitude was noted by Velicogna & Wahr (2006). This suggests that we can achieve an adequate approximation of the total contribution of all mascons to the satellite accelerations by considering a reduced subset of mascons. This will allow estimation of a local mascon solution and reduce the computational requirement in a global solution. To investigate this we repeated the simulated recovery of the GLDAS anomaly in Section 4.2, but this time we estimate only the contributions to the satellite accelerations of a subset of mascons. The number of mascons included is based on the distance between the centre of the mascon under consideration to the mascon directly below the leading GRACE satellite. We will call this the centre-to- centre distance and the mascon below the satellite the nadir mascon. Any mascon pair where the centre-to-centre distance is greater than Figure 8. The contribution of mascons to the gravitational acceleration of a prescribed value will be excluded from eq. (4) and its effect a satellite over the Amazon as a function of the centre-to-centre distance. set to zero. Thus, for each arc only a subset of mascons will be Acceleration is given as a percentage of the acceleration of the nadir mascon. included in the normal equations, reducing the processing time and storage requirement. All the individual arcs will be summed to increasing the centre-to-centre distance to 6000 km improves the ﬁt calculate a monthly solution, as previously. The solutions are global of the degree variances with the original mascon solution, at which ◦ ◦ in extent as each month all mascon parameters are estimated in the point the degree variances for both the 2 and 4 mascon solutions ◦ ◦ monthly solution. Mascon 4 and 2 solutions were generated using closely resemble the global solution. A centre-to-centre distance various centre-to-centre distances. The resulting degree variances of 4000 km or 5000 km could be used, but the solutions have less are plotted in Figs 6 and 7, respectively and compared to those from power at degree 2 and 3 harmonics (longer wavelength). the global mascon solutions and input GLDAS anomaly. The 0 km The contribution to the acceleration experienced by the GRACE centre-to-centre distance solution refers to the contribution of the satellites at distances from the nadir mascon are plotted in Fig. 8,for nadir mascon only to the acceleration experienced by the GRACE a nadir mascon located in the Amazon (similar results are observed satellites. when the nadir mascon is located elsewhere). The accelerations are Fig. 6 reveals that for a centre-to-centre cut-off distance of given as a percentage of the acceleration of the nadir mascon. Fig. 8 1000 km or less the degree variances of the original mascon so- reveals that, as expected, mascons close to the nadir mascon have a lution are not well recovered as the leakage seen in Fig. 2 has not larger contribution to the accelerations experienced by the satellites been fully accounted for. All such solutions have reduced power at than distant mascons. As the centre-to-centre distance increases, higher degrees, which dampens short wavelength information in the the contribution to the satellite acceleration reduces and approaches gravity ﬁelds. The 0 km centre-to-centre distance solution in Fig. 7 zero. Mascons within a centre-to-centre distance of up to 2000 km (2 mascon) is anomalous with more power than expected, especially have the greatest contribution but mascons 2000 km from the nadir at lower degrees due to the leakage of the spherical harmonics used mascon still contribute ∼10 per cent of the acceleration of the nadir to deﬁne the mascon parameters in eq. (3). In both Figs 6 and 7, mascon. Mascons with centre-to-centre distances up to 4000 km are Downloaded from https://academic.oup.com/gji/article/200/1/503/752750 by DeepDyve user on 20 July 2022 510 S. B. Andrews, P. Moore and M. A. King still seen to inﬂuence the accelerations experienced by the satellites. through to 6000 km at which the difference in the rms is only 0.06 cm Above 4000 km the contributions to the accelerations reduce, ex- of EWH to the original (all mascon) solution. A similar result is plaining the relatively small improvements at higher degrees seen observed for more limited testing of the 2 solution (not shown). in Figs 6 and 7. A centre-to-centre cut-off distance of 6000 km will Iterations improve the solution but are not a substitute for using be used in future processing as the difference between that and a an appropriate centre-to-centre cut-off distance and our simulations full global solution is minimal. show iterations were unable to overcome modelling deﬁciencies associated with an insufﬁcient distribution of mascons when using perfect data. Iterations will therefore not be used in the remainder of this paper. Note that this iterative procedure is not intended as a study 4.4 Iteration of the mascons of the use of iterating a correctly parametrized solution. Applying In Section 4.3 a good approximation to a true mass anomaly (sim- iterations to real GRACE observations, as discussed in Luthcke ulated with GLDAS) was obtained by including the contribution of et al. (2013), would be expected to yield improvements due to the all mascons within a centre-to-centre distance of 6000 km. How- separation of signal and noise through enhanced quality control of ever, it is also worth investigating if a similar approximation can the residuals. Sabaka et al. (2010) found that forward modelling be achieved by using a smaller centre-to-centre cut-off distance but improved signal recovery while resulting in KBRR residuals that also iterating the solution. Iteration was shown by Luthcke et al. conformed more closely to zero-mean Gaussian distributions. (2013) to substantially increase signal recovery. Here, we computed an iterated 4 mascon solution for the centre-to-centre distances of 0, 1000 and 2000 km and computed the rms of the differences to 4.5 Simulated recovery of the GLDAS anomaly with the known GLDAS anomaly (Table 1). The rms of the differences simulated stripes for the 3000, 4000, 5000 and 6000 km centre-to-centre solutions are also provided (not iterated) for comparison. Computations to this point utilized noise and error free data. We now Table 1 shows that, while iterating will improve the solution, test the effect of adding realistic noise and errors to the simulated the extent of the improvement is limited. For iterated solutions KBRR observations generated from observed noise and errors in convergence is normally reached by ∼5 steps. The rms improves processing actual GRACE data. Although not directly pertinent to with an increase in the centre-to-centre distance of the solution the simulations we outline our current two stage processing strategy for GRACE data as the ﬁrst stage provides the accelerometer bias Table 1. Rms of the differences between the iterated 4 values and state vectors for the short-arc gravity ﬁeld recovery in mascon solutions and the input GLDAS anomaly in cm of the second stage. The accelerometer bias runs are undertaken in EWH. Rms values are provided over land (and globally). 24 hr arcs. In addition to constant accelerometer biases we also The inclusion of a mascon into the solution depends on the estimate the full state vector for GRACE A and B from KBRR, distance between that mascon and the nadir mascon being GNV1B and accelerometer data in addition to KBRR bias, trend lower than the centre-to-centre distance. and once and twice per rev per orbital revolution (16 sets per arc) Centre-to-centre Iteration Rms difference to GLDAS measurement parameters. This process essentially computes the distance (km) number anomaly (cm EWH) acceleration biases from the GNV1B Cartesian coordinates while 0 1 10.5 (6.4) the KBRR parameters account for unmodelled force effects enabling 2 10.3 (6.1) bad data to be identiﬁed. 3 9.9 (6.0) To generate the KBRR noise and errors we re-processed the ac- 4 9.8 (5.9) celerometer runs estimating now 96 sets of empirical once per rev 5 9.7 (6.0) acceleration parameters along-track and cross-track instead of the 6 9.8 (6.0) KBRR parameters. The empirical accelerations absorbed the longer 7 9.8 (6.2) wavelength gravity ﬁeld mismodelling with the KBRR residuals 8 9.9 (6.3) representing a realistic sample of the noise and errors present in the 1000 1 10.4 (6.2) GRACE data. The empirical acceleration once per rev parameters 2 10.1 (6.2) also lead to slightly different values for the accelerometer biases 3 9.9 (6.2) values and estimated state vectors. The small differences in the ac- 4 9.8 (6.2) 5 9.8 (6.3) celerometer bias values between the runs are assumed to be similar 6 9.8 (6.4) to the differences between the original values and the unknown 7 9.8 (6.5) ‘true’ values. Similarly the orbit used to calculate the KBRR resid- 8 9.9 (6.7) uals will be different to the orbits used in the gravity recovery due 2000 1 9.3 (5.9) to the use of different empirical parameters. The KBRR residuals 2 9.1 (5.8) time-series will also contain errors in the orbit introduced from the 3 8.8 (5.6) inclusion of the GNV1B data as well as the KBRR data (GNV1B 4 8.7 (5.5) data is not used in the 60 min arcs for gravity recovery). Further- 5 8.6 (5.5) more, the KBRR residuals will also contain effects of aliasing of 6 8.6 (5.5) 7 8.6 (5.5) short period signals and errors in the background models which will 8 8.6 (5.5) affect the real GRACE solution. AOD1B (Atmosphere and Ocean 3000 1 8.3 (5.2) de-aliasing Level-1B) product and the GLDAS anomaly ﬁeld were 4000 1 7.7 (4.7) utilized in all computations. 5000 1 7.6 (4.6) Despite best efforts it is likely that residual geographical cor- 6000 1 7.5 (4.5) related terms are present in the KBRR residuals due to the All mascons 1 7.4 (4.4) parametrization failing to absorb all gravity ﬁeld mismodelling. Downloaded from https://academic.oup.com/gji/article/200/1/503/752750 by DeepDyve user on 20 July 2022 Mass change from GRACE: simulated comparison 511 Figure 9. Periodogram of the original KBRR residual time-series (a); periodogram created from a FFT of the phase-shifted reconstructed KBRR time-series (b); sample of the original KBRR residual time-series (c); sample of the phase-shifted reconstructed KBRR residual time-series (d). In mitigation we preserved the characterization of the error signa- tures but redistributed the signatures geographically by applying a fast Fourier transform (FFT) to the residual time-series. Modiﬁed KBRR residuals for each 24 hr were calculated by reconstructing the FFT through keeping the same amplitude of spectral frequency but assigning a random phase. Fig. 9 shows the original and modiﬁed FFT and KBRR residuals. Fig. 9(a) shows the original periodogram, while Fig. 9(b) shows the modiﬁed spectrum used to generate the required KBRR residuals. Fig. 9(c) shows a sample of the orig- inal KBRR residuals, while Fig. 9(d) shows the modiﬁed KBRR residuals. The original and modiﬁed KBRR residuals have simi- lar range and noise characteristics. Fig. 9(a) shows features that are replicated in Fig. 9(b). The low frequencies (i.e. long spatial wavelengths) were investigated and characterized as ﬂicker noise. Higher frequencies were seen to be approximately Gaussian with an increase in power between 1800 and 1900 Hz. The frequency of this signal with increased power corresponds to about twice the separa- tion distance between the GRACE satellites. Without this artefact it would have been possible to generate the spectrum using a com- bination of ﬂicker and Gaussian noise. Tests with other days and months showed that the periodogram was representative. Tracking data, in 60 min arcs, was created with the modiﬁed KBRR residuals and the gravity ﬁeld solution recovered to degree and order 60. The resulting gravity ﬁeld is shown in Fig. 10 which reveals that the recovered solution is dominated by short wavelength, high degree, north–south orientated stripes; Fig. 10(b) shows the same but with 400 km Gaussian smoothing applied. The stripes Figure 10. Simulated stripes: unsmoothed (a); with 400 km Gaussian are typical of those in the monthly GRACE solutions. A plot of smoothing (b). Units of cm of EWH. the degree variance of the simulated stripes in Fig. 11 shows that the simulated stripes begin to dominate the signal between degree 20 simulated stripes above degree 30 was found to be comparable to and 30. Degree 30 is the degree at which regularization is normally those of the unsmoothed CSR solution for 2006 June. Most of the applied to monthly solutions (Bruinsma et al. 2010;Save et al. true geophysical signal is restricted to lower degrees, below degree 2012) as required to reduce the effect of noise and errors in the ∼20 (Wouters & Schrama 2007) but short wavelength signals, as high degree coefﬁcients. The power of the degree variances of the required for hydrological and glaciological basin-level applications, Downloaded from https://academic.oup.com/gji/article/200/1/503/752750 by DeepDyve user on 20 July 2022 512 S. B. Andrews, P. Moore and M. A. King rors with consequences that mirror the errors in actual monthly GRACE ﬁelds. The degree variances of the unconstrained solutions in Fig. 11, show that the spherical harmonic and 2 mascon solution are able to reproduce the GLDAS anomaly with simulated stripes. The global values of the 0 km (50 km) solution in Table 2 (spatial constraint mascon solutions) conﬁrm that the simulated stripes are present in the spherical harmonic and 2 mascon solutions, with rms values of 38.1 and 35.2 cm, respectively. The 4 mascon solution is able to reproduce the degree variance of the GLDAS plus simulated stripes to degree ∼35, at which point it loses power. This loss of power, a consequence of undersampling of the signal, results in a rms of 14.7 cm in Table 2 for the 0 km solution. Again, the dis- tance for the spherical harmonic coefﬁcients refers to the Gaussian smoothing radius, while for the mascon solutions the distance is D of eq. (9). As previous, the distances are not directly comparable. Regardless of the strategy used to recover the simulated ﬁeld, some form of post-processing is required to reduce the effects of noise and errors in the solution. Increasing the Gaussian smoothing Figure 11. Degree variances of the simulated stripes, GLDAS with simu- distance for the spherical harmonic solution reduces the rms towards lated stripes added and the unﬁltered spherical harmonic and unconstrained the original GLDAS anomaly, with 400 km Gaussian smoothing mascon solutions. The original GLDAS model is included for comparison. producing the lowest rms to the GLDAS anomaly (Table 2). The global rms still reduces for 500 km Gaussian smoothing due to smoothing over the oceans, but by 500 km the land signal is over is contained in the higher degrees (Kusche 2007). This simulation smoothed. Applying a lower smoothing radius allows the simulated thus allows us to assess the ability of the different methodologies to stripes to dominate the solution, while a larger smoothing radius maintain the signal at lower degrees while reducing noise and error reduces the amplitude of geophysical signals and removes short in the higher degrees and maximizing recovery at these shorter ◦ ◦ wavelength information. For both the 4 and 2 mascon solutions, wavelengths. increasing the correlation distance D in the constraint matrix reduces The modiﬁed KBRR residuals were added to the original track- the rms until ∼700 km. Over all solutions in Table 2 the 2 mascon ing data that described the GLDAS anomaly in Section 4.2. These solution provides the closest match to the GLDAS anomaly when residuals emulate the spectrum of real GRACE noise and er- Table 2. Rms (cm of EWH) of the differences between the calculated gravity ﬁeld solutions and input GLDAS anomalies in the presence of simulated stripes. Rms values are provided over land (globally). The distance refers to the Gaussian smoothing radius (spherical harmonic solution) or the correlation distance (mascon solutions). Mascon solution constrained using spatial (columns 3 and 4) or basin constraints (columns 5 and 6). Rms difference to GLDAS anomaly (cm EWH) Distance (km) Spatial constraint Basin constraint applied ◦ ◦ ◦ ◦ SH coefﬁcients 4 mascon 2 mascon 4 mascon 2 mascon 0 38.1 (37.1) 14.9 (14.7) 36.4 (35.2) – – 50 37.1 (35.9) 14.9 (14.7) 31.4 (30.2) 14.9 (14.7) 31.8 (30.4) 100 31.4 (30.3) 14.1 (13.6) 18.2 (17.3) 14.2 (13.7) 19.3 (17.7) 150 24.3 (23.1) 12.2 (11.1) 13.1 (11.9) 12.6 (11.2) 14.1 (12.3) 200 17.7 (16.4) 10.9 (9.2) 11.1 (9.5) 11.3 (9.3) 11.8 (9.7) 250 13.2 (11.4) 10.1 (8.0) 10.1 (8.1) 10.5 (8.1) 10.7 (8.3) 300 10.9 (8.4) 9.6 (7.3) 9.5 (7.2) 9.9 (7.4) 10.0 (7.4) 350 10.1 (6.9) 9.3 (6.7) 9.1 (6.7) 9.6 (6.8) 9.5 (6.8) 400 10.0 (6.4) 9.1 (6.4) 8.9 (6.3) 9.3 (6.4) 9.2 (6.4) 450 10.1 (6.2) 8.9 (6.1) 8.8 (6.0) 9.2 (6.2) 9.0 (6.1) 500 10.3 (6.2) 8.9 (5.9) 8.7 (5.8) 9.0 (5.9) 8.9 (5.9) 550 10.5 (6.3) 8.8 (5.8) 8.6 (5.6) 8.9 (5.8) 8.8 (5.7) 600 10.7 (6.3) 8.8 (5.7) 8.6 (5.5) 8.9 (5.6) 8.8 (5.5) 650 10.9 (6.4) 8.8 (5.6) 8.6 (5.4) 8.8 (5.5) 8.7 (5.4) 700 11.1 (6.5) 8.8 (5.5) 8.6 (5.4) 8.8 (5.5) 8.7 (5.4) 750 11.3 (6.5) 8.8 (5.5) 8.7 (5.3) 8.8 (5.4) 8.7 (5.3) 800 11.5 (6.6) 8.8 (5.4) 8.7 (5.3) 8.8 (5.3) 8.7 (5.3) 850 11.6 (6.7) 8.8 (5.4) 8.7 (5.3) 8.8 (5.3) 8.7 (5.2) 900 11.8 (6.7) 8.9 (5.4) 8.7 (5.2) 8.8 (5.3) 8.7 (5.2) 950 11.9 (6.8) 8.9 (5.4) 8.8 (5.2) 8.8 (5.3) 8.7 (5.2) 1000 12.0 (6.9) 8.9 (5.3) 8.8 (5.2) 8.8 (5.2) 8.7 (5.2) 1050 12.1 (6.9) 8.9 (5.3) 8.8 (5.2) 8.8 (5.2) 8.7 (5.1) 1100 12.2 (7.0) 9.0 (5.3) 8.9 (5.2) 8.8 (5.2) 8.8 (5.1) Downloaded from https://academic.oup.com/gji/article/200/1/503/752750 by DeepDyve user on 20 July 2022 Mass change from GRACE: simulated comparison 513 exploited the full advantage of the mascon processing strategy; namely that the constraint matrix in eq. (5) can be adapted to con- strain mascons that share geophysical similarities as in Luthcke et al. (2006b). Thus, in addition to separating mascons into either land or ocean, we now also separate the major hydrological drainage basins deﬁned by the total runoff integrating pathways (TRIP) data set(Oki&Sud 1998). Major Greenland and Antarctic drainage basins are adopted according to Zwally et al. (2012). The result is a constraint matrix where constraints are imposed between all land mascons within the same (hydrological or glaciological) drainage basin, with mascons in different basins uncorrelated. Not all land mascons are within a major hydrological or glaciological basin. Land mascons that are not classiﬁed as belonging to a basin were assumed to be correlated, but in most instances d of eq. (8) or (9) IJ reduces the correlation to ∼0. All ocean mascons are correlated as previous. Table 2 (Basin constraint mascon solutions) summarizes the rms of the differences to the GLDAS plus simulated stripes. For both Figure 12. The degree variances of GLDAS (with and without simulated mascon solutions the use of a basin constraint, while offering a ◦ ◦ stripes) and the spherical harmonic and mascon solutions. 2 and 4 mascon similar or small improvement in the rms over land, offers an im- are the constrained mascon solutions without the basin constraints applied provement in global rms values over the solutions without basin ◦ ◦ while 2 and 4 basin are the constrained mascon solutions with basin constraints. An increase in the correlation distance (800–900 km) constraints applied. The spherical harmonic solution has 400 km Gaussian is required to obtain these improvements due to the reduced num- smoothing applied; the mascon solutions are constrained using a correlation ber of constraints based on mascons pairs that share geophysi- distance of 700 km in the constraint matrix (without basin constraints) and 900 km (with basin constraint). cal similarities. The improved global rms is due to reduction of the simulated stripes over the ocean. As ocean mascons are con- a correlation distance of ∼650 km is used in the constraint matrix. strained together, the increased correlation distance allows the effect For the 4 mascon solution we investigated the effect on the weight of the simulated stripes to be further mitigated. The improvement matrix V of varying the value of S in the eq. (9). The weight matrix is evident in the degree variances of the mascon solutions (basin) V is dependent on S and D.Using S = 0.001, the lowest rms (8.8 cm plotted in Fig. 12. Comparing mascon solutions with and without EWH) was achieved using D = 700 km. When S = 0.01 (0.0001) basin constraints, the power of the resulting mascon solutions from were used rms values of 9.8 (10.3) were obtained. For S = 0.01 the degrees 30 to 60 more closely matches the input signal when the lowest rms was 8.8 cm when D = 250 km, while for S = 0.0001 basin constraint is applied. The constraints have successfully re- an rms of 9.0 cm was obtained when D = 1500 km. We therefore duced the noise at these degrees while reducing loss of actual signal. concluded, within the scope of this limited testing, that S = 0.001 Retention of short wavelength signal is of fundamental importance to be the optimal scale factor for the simulations and will use this for improving basin-scale hydrological and glaciological recovery is all further processing. and for estimating basin-level mass ﬂux. Fig. 12 shows the degree variances of the spherical harmonic so- The ﬁnal recovered GLDAS anomalies are shown in Fig. 13 for lution with 400 km Gaussian smoothing. It is not possible to directly the spherical harmonic solution with 400 km Gaussian smoothing relate the Gaussian smoothing distance to the mascon correlation together with the 2 mascon solution using a 800 km correlation dis- distance so we show the mascon solutions for D = 700 km. This tance and basin constraints. The input GLDAS anomaly is provided value was chosen as a comparison of the degree variances of the for reference. Comparing Figs 13(c) to (a) validates the 2 mascon solutions given in Table 2 revealed it removes a comparable amount solution. The mascon methodology offers a clear improvement in of noise to the spherical harmonic solution at the higher degrees. the recovery of TVG and suppression of simulated stripes compared It also gives one of the lowest rms values to the GLDAS anomaly. with a spherical harmonic solution with Gaussian smoothing. The Fig. 12 reveals that all solutions have removed a similar amount of reduction of the simulated stripes over the ocean is particularly evi- noise and error above degree 30, and contain less power at these de- dent. Unlike spherical harmonics, mascons allow constraints within grees compared to the GLDAS signal. The mascon solutions more areas that share real world geophysical properties. Any number of closely resemble the power of the GLDAS signal between ∼degree constraints can be applied with the aim of improving the overall 15 and 30 than the spherical harmonic solution. The reduced power recovery of the TVG ﬁeld. of the spherical harmonic signal between ∼degree 15 and 30 could be partly related to the degree dependency of the Gaussian ﬁlter which Han et al. (2005) suggested could suppress information over 4.7 Simulated recovery of Antarctic mass change signal this range. The mascon solution has therefore been able to recover more of the original signal at degrees ∼15–30 while removing a One of the scientiﬁc advances provided by the GRACE mission is comparable amount of noise and error at higher degrees. recovery of continent-wide ice sheet mass change estimates, such as over Antarctica, with temporal resolutions of 30 days or less. So far our simulations have focused on recovery of a GLDAS hydro- 4.6 Basin constraints logical anomaly. Given the scientiﬁc interest in contribution of ice Mascon solutions with spatial constraints are able to remove more sheets to sea-level change and the convergence of orbital tracks at noise and error and retain more signal than spherical harmonics the pole we now simulate recovery of a realistic Antarctica mass solutions with Gaussian smoothing. However, thus far we have not signal along with the GLDAS anomaly and error/noise product Downloaded from https://academic.oup.com/gji/article/200/1/503/752750 by DeepDyve user on 20 July 2022 514 S. B. Andrews, P. Moore and M. A. King Figure 13. GLDAS anomaly for 2006 June to degree and order 120 (a); Re- covered GLDAS anomalies for the spherical harmonic solution with 400 km Gaussian smoothing to degree and order 60 (b); and 2 mascon solution with a basin constraint applied, with 800 km correlation distance (c). Units of cm EWH. described in Section 4.5. The simulated Antarctica mass signal was derived from the basin mass change rates provided in Table S1 of the supplementary material of King et al. (2012). While King et al. (2012) calculated GRACE solutions to degree and order 60, forward modelling across deﬁned basins resulted in a higher effective reso- lution. The mass change rates used have been corrected for leakage but still contain the GIA signal that is observable from GRACE. Using the same drainage basins adopted by King et al. (2012), as deﬁned using ICESat data (Zwally et al. 2012), the mass change rates for each basin were converted to a total mass change in cm of EWH using the density of water and the area of each drainage Figure 14. Simulated Antarctic mass change signal (a); spherical harmonic basin. As in King et al. (2012) basins 25 and 26, which cover the solution ﬁltered with a 350 km Gaussian smoother (b); and 2 mascon solu- northern tip of the Antarctica Peninsula, were merged into a single tion with basin constraint with correlation distance D = 800 km (c). Units basin. The resulting mass changes for each basin were then ex- of cm of EWH. panded into spherical harmonics to degree and order 60 and added to the mean background ﬁeld to create the KBRR tracking data. Simulated stripes are included as previous. The simulated Antarctic mass change signal can be seen in Fig. 14. Plotting the unﬁltered Downloaded from https://academic.oup.com/gji/article/200/1/503/752750 by DeepDyve user on 20 July 2022 Mass change from GRACE: simulated comparison 515 Table 3. Rms (cm of EWH) of the differences between the calculated gravity ﬁeld solutions and the input GLDAS anomaly with Antarctic mass change signal in the presence of different simulated noise products. Rms values are provided over land (and globally). The distance refers to the Gaussian smooth- ing radius (spherical harmonic coefﬁcient solutions) or the correlation distance (mascon solutions). All mascon solutions utilize basin constraints described in Section 4.6. For the thermal mascon solutions the value S in eq. (9) was increased from 0.001 to 0.01. Rms difference to GLDAS anomaly and Antarctic mass change (cm EWH) Distance (km) Main period Resonance period Reduced thermal control ◦ ◦ ◦ SH 2 mascon SH 2 mascon SH 2 mascon 200 17.7 (16.4) 11.8 (9.8) 20.3 (19.8) 8.1 (5.3) 45.2 (45.5) 13.9 (11.0) 250 13.2 (11.4) 10.7 (8.3) 14.8 (13.8) 8.0 (5.0) 29.9 (29.7) 12.4 (9.3) 300 10.9 (8.4) 10.0 (7.4) 11.6 (9.8) 8.0 (4.9) 20.6 (19.5) 11.7 (8.2) 350 10.1 (7.0) 9.5 (6.8) 10.3 (7.6) 8.0 (4.9) 15.6 (13.7) 11.3 (7.6) 400 10.0 (6.4) 9.2 (6.4) 10.0 (6.6) 8.1 (4.9) 13.2 (10.6) 11.0 (7.2) 450 10.2 (6.3) 9.0 (6.1) 10.1 (6.2) 8.2 (4.9) 12.2 (9.0) 10.9 (7.0) 500 10.4 (6.3) 8.9 (5.9) 10.3 (6.2) 8.3 (4.9) 11.7 (8.1) 10.8 (6.8) 550 10.6 (6.3) 8.8 (5.7) 10.6 (6.2) 8.3 (4.9) 11.6 (7.6) 10.8 (6.6) 600 10.8 (6.4) 8.8 (5.5) 10.8 (6.3) 8.4 (4.9) 11.6 (7.4) 10.8 (6.6) 650 11.0 (6.4) 8.7 (5.4) 11.1 (6.4) 8.5 (5.0) 11.6 (7.2) 10.8 (6.5) 700 11.2 (6.5) 8.7 (5.4) 11.3 (6.5) 8.6 (5.0) 11.7 (7.2) 10.8 (6.5) 750 11.4 (6.6) 8.7 (5.3) 11.5 (6.6) 8.6 (5.0) 11.8 (7.1) 10.9 (6.4) 800 11.6 (6.7) 8.7 (5.3) 11.7 (6.7) 8.7(5.0) 11.9 (7.1) 10.9 (6.4) 850 11.7 (6.7) 8.7 (5.2) 11.8 (6.8) 8.8 (5.1) 12.0 (7.1) 11.0 (6.4) 900 11.9 (6.8) 8.7 (5.2) 12.0 (6.9) 8.8 (5.1) 12.1 (7.1) 11.0 (6.4) 950 12.0 (6.9) 8.7 (5.2) 12.1 (6.9) 8.9 (5.1) 12.2 (7.1) 11.1 (6.4) 1000 12.1 (6.9) 8.7 (5.2) 12.2 (7.0) 9.0 (5.2) 12.3 (7.2) 11.1 (6.4) 1050 12.2 (7.0) 8.7 (5.1) 12.3 (7.0) 9.0 (5.2) 12.3 (7.2) 11.2 (6.5) 1100 12.3 (7.0) 8.8 (5.1) 12.4 (7.1) 9.1 (5.2) 12.4 (7.2) 11.2 (6.5) degree variances (not shown) reveals that the spherical harmonic 4.8 GRACE resonance and thermal control issues and 2 mascon solution were able to reproduce the Antarctica mass Periods where the GRACE satellites are in 3, 4 and 7 day repeat signal added to the GLDAS anomaly and simulated stripes. Appli- orbits modify the noise characteristics of the GRACE solutions cation of a single value of mass change per basin is a simpliﬁcation (Save et al. 2012), while recent GRACE solutions are affected by on the real spatial pattern of change, but this assumption is sufﬁcient poor satellite thermal control. For replicating the GRACE errors for our purposes here. and studying the correlations between the spherical harmonic coef- Solutions were estimated for the spherical harmonics and 2 ﬁcients in CSR RL05 solutions, covariance matrices are provided mascon methodologies with the rms values shown in Table 3 by CSR for the three characteristic periods in the GRACE mission: (main period). The drainage basins used in the mascon solution in early mission (2002 February–2005 May), main part of the mission Section 4.6 result in mascons within the same Antarctic drainage (2007 July–2010 December) and recent months (2011 February– ) basin being constrained together, but with mascons in different which are affected by poor thermal control. Only three covariance Antarctic basins uncorrelated. The mascon results are shown in matrices are required as the monthly variations within these periods Fig. 14(c) revealing this solution has been able to recover most of are minor as to be unimportant (John C Ries, personal communica- the input signal (Fig. 14a) while the spherical harmonic solution tion, 2014). An additional covariance matrix is available for months ﬁltered with 400 km Gaussian smoothing (Fig. 14b) has recovered where the satellite is in a repeat orbit (resonance). The effect these the spatial extent of the signal but with reduced amplitudes. This different periods have on the noise and errors is investigated through is important to consider within future GRACE studies of Antarctic recovery of the GLDAS anomaly and Antarctica mass change sig- ice mass change. The lowest rms of the differences to the input nal for two additional sets of simulated noise. The GLDAS anomaly signal (in EWH) over Antarctica is 8.7 cm for the mascon solution and Antarctica mass change signal are the same as in previous sim- (800 km correlation distance) and 10.4 cm with spherical harmonics ulations. The ﬁrst set of simulated noise relates to September 2004, (400 km Gaussian smoothing). As with the global hydrological sim- when the GRACE satellites were in 61:4 resonance (Klokocn ˇ ´ık et al. ulations, we are able to recover more of the signal using the mascon 2013). The second set relates to 2011 March, which corresponds to methodology. a period when the GRACE satellites are affected by poor thermal As a ﬁnal comparison we computed the mascon solution with control. KBRR residuals for 2004 September and 2011 March with only land/ocean constraints as applied in Section 4.5. The rms of the simulated noise were derived using the method of Section 4.5. We difference to the input Antarctica signal was 8.8 cm. This highlights computed the modiﬁed FFT and KBRR residuals for 2004 Septem- the advantage offered by the mascon methodology when real-world ber and 2011 March (as in Fig. 9). For 2004 September (not shown) geophysical properties are used to deﬁne the constraints between these revealed that the noise and errors in GRACE during resonance the mascons. We note, however, that a simple basin-wide constraint are indistinguishable from that of the main part of the mission. The may not be an accurate reﬂection of changes within a basin and more 2011 March residuals are shown in Fig. 15, revealing an increase in sophisticated treatment, such as increasing the numbers of basins the amplitude of the KBRR residuals; compare Figs 15(c) and (d) (e.g. Luthcke et al. 2013) based on ice dynamic considerations, may with Figs 9(c) and (d). This increase in amplitude is the result of further improve analysis of real GRACE data. Downloaded from https://academic.oup.com/gji/article/200/1/503/752750 by DeepDyve user on 20 July 2022 516 S. B. Andrews, P. Moore and M. A. King Figure 15. Periodogram of the poor thermal control KBRR residual time-series (a) and FFT phase shifted reconstructed KBRR time-series (b); sample of the poor thermal control KBRR residual time-series (c); sample of the reconstructed KBRR residual time-series (d). the poor thermal control, which we would expect to increase from S in eq. (9) was increased from 0.001 to 0.01 which increases the 2011 onwards. effect the constraint matrix has on the solution. For both the mascon Table 3 summarizes the rms of the difference to the GLDAS and and spherical harmonic solutions a larger constraint matrix is re- Antarctica mass loss signal in the presence of simulated stripes from quired to deal with the noise and errors resulting from the poor three periods (main, resonance and thermal control). The simulated thermal control. While the recovery in the presence of the thermal stripes for the main period are described in Section 4.5. The rms stripes is not as favourable for both methodologies, when compared values are provided for the spherical harmonic and 2 mascon solu- to the original or resonance periods, the mascon methodology is still tions. In the resonance period, increasing the Gaussian smoothing able to recover more of the input signal than the spherical harmonic until 400 km reduces the rms at which point the signal recovered methodology while also reducing the noise over the ocean. is comparable to the main period; the larger global rms suggests more of the stripes remain over the oceans. For the mascon solution the lowest rms obtained is obtained when a correlation distance of 5CONCLUSION 250 km is used, which is a much lower than the distance required for the 2006 June simulation of the main period. The rms is also much Through the use of simulations we have demonstrated the sensitivity lower, suggesting the mascon methodology is able to deal well with of GRACE to mass over a large area, with the effect on the KBRR periods of resonance. In 61:4 resonance the satellite ground-track residuals depending on the amplitude of the signal. Building on repeats every 61 orbital revolutions in the time the Earth revolves this we showed that recovery of global mass distributions require 4 times relative to the ascending node. This gives 122 repeating as- contributions to the satellite accelerations of all mascons within cending and descending passes every 4 days spaced on average every approximately 6000 km of the nadir mascon. This has implications 3 in longitude. This spatial resolution is adequate for a spherical when estimating a non-global mascon solution, and is useful in harmonic ﬁeld to degree and order 60. Similarly, a mascon correla- reducing computational effort for high resolution global mascon tion distance of 250 km, equivalent to the distance between centres solutions. We found that the use of iterations cannot compensate for of 2 equatorial mascons, provides the ties between neighbouring an insufﬁcient distribution of mascons. mascons. It would appear that the resonance facilitates separation Simulated recovery of a noise and error free GLDAS mass of spatial and temporal errors in our simulation leading to more anomaly showed that there is no quantiﬁable difference between accurate surface mass recovery than in non-resonant periods. a mass anomaly solution recovered from KBRR data using spheri- In the period of reduced thermal control both the spherical har- cal harmonics and a 2 equal area mascon solution; both recovered monic and mascon solutions are degraded. The lowest rms (11.6 cm the input GLDAS anomaly. A 4 mascon solution does not have the EWH) for the spherical harmonic solution is obtained with 600 km required resolution to recover the gravity ﬁeld solution to degree Gaussian smoothing. Even with this smoothing radius, the rms and order 60 due to the number of mascon parameters estimated obtained is larger than for the main and resonance periods. For being less than the gravity ﬁeld parameters required; this causes a the mascon solution, the lowest rms (10.8) is obtained with a corre- loss of higher degree information in the gravity ﬁeld. Recovery of lation distance D = 600 km, compared to 8.7 (8.0) for the solutions the GLDAS mass anomaly described to degree and order 120 re- with the original (resonance) noise. To achieve this rms the value vealed that all methodologies are affected by omission errors, with Downloaded from https://academic.oup.com/gji/article/200/1/503/752750 by DeepDyve user on 20 July 2022 Mass change from GRACE: simulated comparison 517 the additional short wavelength (high degree) information not being Chen, J.L., Wilson, C.R., Tapley, B.D., Blankenship, D. & Young, D., 2008. Antarctic regional ice loss rates from GRACE, Earth planet. Sci. Lett., captured. 266, 140–148. The recovery of the GLDAS anomaly with simulated stripes Dunn, C. et al., 2003. Instrument of GRACE: GPS augments gravity mea- reveals the advantage of the mascon solution over a spherical har- surements, GPS World, 14, 16–28. monic solution smoothed with a Gaussian ﬁlter when attempting to Farrell, W.E., 1972. Deformation of earth by surface loads, Rev. Geophys. reduce high degree and order noise. The use of the mascon spatial Space Phys., 10, 761–797. constraint matrix allowed more signal to be preserved up to degree Han, S.-C., Shum, C.K., Jekeli, C., Kuo, C.-Y., Wilson, C. & Seo, K.-W., ∼35, while the addition of constraints between mascon parameters 2005. Non-isotropic ﬁltering of GRACE temporal gravity for geophysical that share geophysical similarities resulted in a further reduction signal enhancement, Geophys. J. Int., 163(1), 18–25. of the lost signal at all degrees. This improvement is further evi- Han, S.-C., Sauber, J., Luthcke, S.B., Ji, C. & Pollitz, F.F., 2008. Implica- dent in the simulated recovery of the Antarctic mass signal which tions of postseismic gravity change following the great 2004 Sumatra- Andaman earthquake from the regional harmonic analysis of GRACE also conﬁrms the advantage of the mascon solution over the spher- intersatellite tracking data, J. geophys. Res.: Solid Earth, 113, B11413, ical harmonic recovery validating the use of this methodology in doi:10.1029/2008JB005705. polar regions. Inclusion of temporal constraints could result in fur- Han, S.-C., Kim, H., Yeo, I.-Y., Yeh, P., Oki, T., Seo, K.-W., Alsdorf, D. & ther improvements, but a level of assumption would be required on Luthcke, S.B., 2009. Dynamics of surface water storage in the Amazon plausible variations of mass over time. inferred from measurements of inter-satellite distance change, Geophys. The recovery of the GLDAS anomaly and Antarctica mass change Res. Lett., 36, L09403, doi:10.1029/2009GL037910. signal with simulated stripes, created during periods when GRACE Jacob, T., Wahr, J., Pfeffer, W.T. & Swenson, S., 2012. Recent contributions was in resonance and suffering from poor thermal control, revealed of glaciers and ice caps to sea level rise, Nature, 482, 514–518. that in both cases the mascon methodology was able to recover more King, M.A., Bingham, R.J., Moore, P., Whitehouse, P.L., Bentley, M.J. & of the input signal than using spherical harmonic methodology. This Milne, G.A., 2012. Lower satellite-gravimetry estimates of Antarctic sea- level contribution, Nature, 491, 586–589. advantage was especially evident during periods of resonance. The Klees, R., Revtova, E.A., Gunter, B.C., Ditmar, P., Oudman, E., Winsemius, use of the mascon methodology also reduced the effect of stripes H.C. & Savenije, H.H.G., 2008. The design of an optimal ﬁlter for monthly over the ocean. However, both the spherical harmonic and mascon GRACE gravity models, Geophys. J. Int., 175, 417–432. solutions were degraded during periods of reduced thermal control. Klokocn ˇ ´ık, J., Gooding, R.H., Wagner, C.A., Kostelecky, ´ J. & Bezdek, ˇ A., In particular, the reduced thermal control necessitated tighter con- 2013. The use of resonant orbits in satellite geodesy: a review, Surv. straints between the various mascons than during the main part of Geophys., 34, 43–72. the mission. Klosko, S., Rowlands, D., Luthcke, S., Lemoine, F., Chinn, D. & Rodell, Finally, while the values used for S and D were found to be optimal M., 2009. Evaluation and validation of mascon recovery using GRACE in this paper, they may differ for users depending on parametrization, KBRR data with independent mass ﬂux estimates in the Mississippi Basin, and are therefore provided as a guide. Furthermore, adaptations to J. Geod., 83, 817–827. Kusche, J., 2007. Approximate decorrelation and non-isotropic smoothing of the S and D values appear to be necessary as the GRACE mission time-variable GRACE-type gravity ﬁeld models, J. Geod., 81, 733–749. ages. Lemoine, F.G., Luthcke, S.B., Rowlands, D.D., Chinn, D.S., Klosko, S.M. & Cox, C.M., 2007. The use of mascons to resolve time-variable gravity from GRACE, in Dynamic Planet, pp. 231–236, eds Tregoning, P. & ACKNOWLEDGEMENTS Rizos, C., Springer. This work was funded by a National Centre for Earth Observation Luthcke, S.B., Rowlands, D.D., Lemoine, F.G., Klosko, S.M., Chinn, D. & (NCEO) award in gravimetric geodesy and partially supported by McCarthy, J.J., 2006a. Monthly spherical harmonic gravity ﬁeld solutions COST Action ES0701 ‘Improved constraints on models of Glacial determined from GRACE inter-satellite range-rate data alone, Geophys. Res. Lett., 33, L02402, doi:10.1029/2005GL024846. Isostatic Adjustment’. We are gratefully for the GRACE Level- Luthcke, S.B. et al., 2006b. Recent Greenland ice mass loss by 1B products produced by the Jet Propulsion Laboratory. MAK is drainage system from satellite gravity observations, Science, 314, a recipient of an Australian Research Council Future Fellowship 1286–1289. (project number FT110100207). The authors would like to thank Luthcke, S.B., Arendt, A.A., Rowlands, D.D., McCarthy, J.J. & Larsen, three anonymous reviewers for comments which helped improve C.F., 2008. Recent glacier mass changes in the Gulf of Alaska region this manuscript. from GRACE mascon solutions, J. Glaciol., 54, 767–777. Luthcke, S.B., Sabaka, T.J., Loomis, B.D., Arendt, A.A., McCarthy, J.J. & Camp, J., 2013. Antarctica, Greenland and Gulf of Alaska land-ice REFERENCES evolution from an iterated GRACE global mascon solution, J. Glaciol., 59, 613–631. Awange, J.L., Fleming, K.M., Kuhn, M., Featherstone, W.E., Heck, B. & Moore, P., Boomkamp, H.J., Carnochan, S. & Walmsley, R.J., 1999. Anjasmara, I., 2011. On the suitability of the 4 degree x 4 degree GRACE FAUST: multi-satellite orbital dynamics software, Adv. Space Res., 23, mascon solutions for remote sensing Australian hydrology, Rem. Sens. 785–795. Environ., 115, 864–875. Muller, P.M. & Sjogren, W.L., 1968. Mascons: lunar mass concentrations, Bettadpur, S., 2012. UTCSR Level-2 Processing Standards Document for Science, 161, 680–684. Level-2 Product Release 0005, pp. 16, GRACE 327–742. University of Oki, T. & Sud, Y.C., 1998. Design of total runoff integrating pathways Texas at Austin, CSR Publ: GR-12-xx, Rev. 4.0. (TRIP)—a global river channel network, Earth Interact., 2, 1–37. Bruinsma, S., Lemoine, J.-M., Biancale, R. & Vales, ` N., 2010. CNES/GRGS Pritchard, H.D., Luthcke, S.B. & Fleming, A.H., 2011. Understanding 10-day gravity ﬁeld models (release 2) and their evaluation, Adv. Space ice-sheet mass balance: progress in satellite altimetry and gravimetry, Res., 45, 587–601. J. Glaciol., 56, 1151–1161. Chambers, D.P. & Bonin, J.A., 2012. Evaluation of release-05 GRACE time- Ramillien, G., Biancale, R., Gratton, S., Vasseur, X. & Bourgogne, S., variable gravity coefﬁcients over the ocean, Ocean Sci., 8, 859–868. 2011. GRACE-derived surface water mass anomalies by energy in- Chambers, D.P., Wahr, J., Tamisiea, M.E. & Nerem, R.S., 2010. Ocean mass tegral approach: application to continental hydrology, J. Geod., 85, from GRACE and glacial isostatic adjustment, J. geophys. Res., 115, 313–328. B11415, doi:10.1029/2010JB007530. Downloaded from https://academic.oup.com/gji/article/200/1/503/752750 by DeepDyve user on 20 July 2022 518 S. B. Andrews, P. Moore and M. A. King Ries, J.C., Bettadpur, S., Poole, S. & Richter, T., 2011. Mean background Seo, K.-W., Wilson, C.R., Chen, J. & Waliser, D.E., 2008. GRACE’s spatial gravity ﬁelds for GRACE processing, in Proceedings of the GRACE aliasing error, Geophys. J. Int., 172, 41–48. Science Team Meeting, Austin, TX. Shepherd, A. et al., 2012. A reconciled estimate of ice-sheet mass balance, Riva, R.E.M. et al., 2009. Glacial isostatic adjustment over Antarctica from Science, 338, 1183–1189. combined ICESat and GRACE satellite data, Earth planet. Sci. Lett., 288, Swenson, S. & Wahr, J., 2002. Methods for inferring regional surface-mass 516–523. anomalies from Gravity Recovery and Climate Experiment (GRACE) Rodell, M. et al., 2004. The global land data assimilation system, Bull. Am. measurements of time-variable gravity, J. geophys. Res., 107, 2193, Meteor. Soc., 85(3), 381–394. doi:10.1029/2001JB000576. Rowlands, D.D., Ray, R.D., Chinn, D.S. & Lemoine, F.G., 2002. Short- Swenson, S. & Wahr, J., 2006. Post-processing removal of corre- arc analysis of intersatellite tracking data in a gravity mapping mission, lated errors in GRACE data, Geophys. Res. Lett., 33, L08402, J. Geod., 76, 307–316. doi:10.1029/2005GL025285. Rowlands, D.D., Luthcke, S.B., Klosko, S.M., Lemoine, F.G.R., Chinn, D.S., Tapley, B.D., Bettadpur, S., Ries, J.C., Thompson, P.F. & Watkins, M.M., McCarthy, J.J., Cox, C.M. & Anderson, O.B., 2005. Resolving mass ﬂux at 2004a. GRACE measurements of mass variability in the Earth system, high spatial and temporal resolution using GRACE intersatellite measure- Science, 305, 503–505. ments, Geophys. Res. Lett., 32, L04310, doi:10.1029/2004GL021908. Tapley, B.D., Bettadpur, S., Watkins, M. & Reigber, C., 2004b. The gravity Rowlands, D.D., Luthcke, S.B., McCarthy, J.J., Klosko, S.M., Chinn, D.S., recovery and climate experiment: mission overview and early results, Lemoine, F.G., Boy, J.P. & Sabaka, T.J., 2010. Global mass ﬂux solutions Geophys. Res. Lett., 31, 503–505. from GRACE: a comparison of parameter estimation strategies—mass Velicogna, I. & Wahr, J., 2006. Measurements of time-variable gravity show mass loss in Antarctica, Science, 311, 1754–1756. concentrations versus Stokes coefﬁcients, J. geophys. Res., 115, B01403, doi:10.1029/2009JB006546. Wahr, J., Molenaar, M. & Bryan, F., 1998. Time variability of the Earth’s Sabaka, T.J., Rowlands, D.D., Luthcke, S.B. & Boy, J.P., 2010. Improving gravity ﬁeld: hydrological and oceanic effects and their possible detection global mass ﬂux solutions from Gravity Recovery and Climate Experi- using GRACE, J. geophys. Res., 103, 30 205–30 229. ment (GRACE) through forward modeling and continuous time correla- Wahr, J., Swenson, S. & Velicogna, I., 2006. Accuracy of GRACE mass tion, J. geophys. Res., 115, B11403, doi:10.1029/2010JB007533. estimates, Geophys. Res. Lett., 33, L06401, doi:10.1029/2005GL025305. Save, H., Bettadpur, S. & Tapley, B., 2012. Reducing errors in the GRACE Wahr, J., Swenson, S., Zlotnicki, V. & Velicogna, I., 2004. Time-variable gravity solutions using regularization, J. Geod., 86, 695–711. gravity from GRACE: ﬁrst results, Geophys. Res. Lett., 31, L11501, Schrama, E.J.O. & Wouters, B., 2011. Revisiting Greenland ice sheet doi:10.1029/2004GL019779. mass loss observed by GRACE, J. geophys. Res., 116, B02407, Wouters, B. & Schrama, E.J.O., 2007. Improved accuracy of GRACE gravity doi:10.1029/2009JB006847. solutions through empirical orthogonal function ﬁltering of spherical har- Schrama, E.J.O., Wouters, B. & Rietbroek, R., 2014. A mascon ap- monics, Geophys. Res. Lett., 34, L23711, doi:10.1029/2007GL032098. proach to assess ice sheet and glacier mass balances and their uncer- Wouters, B., Chambers, D. & Schrama, E.J.O., 2008. GRACE observes tainties from GRACE data, J. geophys. Res.: Solid Earth, 119, 6048– small-scale mass loss in Greenland, Geophys. Res. Lett., 35, L20501, 6066. doi:10.1029/2008GL034816. Seo, K.W., Wilson, C.R., Famiglietti, J.S., Chen, J.L. & Rodell, M., Zwally, H.J., Giovinetto, M.B., Beckley, M.A. & Sab, J.L., 2012. 2006. Terrestrial water mass load changes from Gravity Recovery Antarctic and Greenland drainage systems, GSFC Cryospheric Sci- and Climate Experiment (GRACE), Water Resour. Res., 42, W05417, ences Laboratory, Available at: http://icesat4.gsfc.nasa.gov/cryo_data/ doi:10.1029/2005WR004255. ant_grn_drainage_systems.php, last accessed 9 December 2013.

Journal

Geophysical Journal International – Oxford University Press

Published: Jan 1, 2015

Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 7-Day Trial for You or Your Team.

Learn More →

Mass change from GRACE: a simulated comparison of Level-1B analysis techniques

Mass change from GRACE: a simulated comparison of Level-1B analysis techniques

Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 7-Day Trial for You or Your Team.

Learn More →

Mass change from GRACE: a simulated comparison of Level-1B analysis techniques

Mass change from GRACE: a simulated comparison of Level-1B analysis techniques

References (55)

Abstract

Journal

Recommended Articles

There are no references for this article.

Our policy towards the use of cookies