Accounting for residual errors in atmosphere–ocean background models applied in satellite gravimetry

The Atmosphere and Ocean non-tidal De-aliasing Level-1B (AOD1B) product is widely used in precise orbit determination and satellite gravimetry to correct for transient effects of atmosphere–ocean mass variability that would otherwise alias into monthly mean global gravity ﬁelds. The most recent release is based on the global ERA5 reanalysis and ECMWF operational data together with simulations from the general ocean circulation model MPIOM consistently forced with ﬁelds from the corresponding atmospheric dataset. As background models are inevitably imperfect, residual errors will consequently propagate into the resulting geodetic products. Accounting for uncertainties of the background model data in a statistical sense, however, has been shown before to be a useful approach to mitigate the impact of residual errors leading to temporal aliasing artefacts. In light of the changes made in the new release RL07 of AOD1B, previous uncertainty assessments are deemed too pessimistic and thus need to be revisited. We here present an analysis of the residual errors in AOD1B RL07 based on ensemble statistics derived from different atmospheric reanalyses, including ERA5, MERRA2 and JRA55. For the oceans, we investigate the impact of both the forced and intrinsic variability through differences in MPIOM simulation experiments. The atmospheric and oceanic information is then combined to produce a new time-series of true errors, called AOe07, which is applicable in combination with AOD1B RL07. AOe07 is further complemented by a new spatial error variance–covariance matrix. Results from gravity ﬁeld recovery simulation experiments for the planned Mass-Change and Geosciences International Constellation (MAGIC) based on GFZ’s EPOS software demonstrate improvements that can be expected from rigorously implementing the newly available stochastic information from AOD1B RL07 into the gravity ﬁeld estimation process.


Introduction
For over two decades now, the satellite gravimetry missions GRACE (Tapley et al. 2004) and GRACE-FO (Landerer et al. 2020) have been monitoring and are continuing to monitor large-scale mass changes on Earth.The twin satellites are tracking ice mass loss in both Greenland (Velicogna and Wahr 2005;Sasgen et al. 2020) and Antarctica (Velicogna et al. 2014(Velicogna et al. , 2020)), changes in terrestrial water storage (Rodell et al. 2018) including the severity of drought (Boergens et al. 2020), and also sea-level change and ocean bottom pressure variations related to internal ocean dynamics (Hamlington et al. 2020;Dobslaw et al. 2020).All those processes are characterised by spatial divergence in mass transports in the Earth system that are well resolved by the monthly gravity field solutions obtained from satellite gravimetry.Mass changes in e.g. the atmosphere and ocean, however, also have significant variations at much shorter, i.e. sub-monthly, time scales.Without prior information, these high-frequency mass transport signals would degrade the monthly gravity field solutions through the effects of temporal aliasing.For this reason, they are usually accounted for in the gravity field estimation by application of a priori background model data.
Non-tidal variations in the atmosphere and oceans are routinely subtracted in satellite gravimetry processing through the Atmosphere and Ocean De-Aliasing Level-1B (AOD1B) data product (Shihora et al. 2022a, b) specifically prepared within the US-German Science Data System of the GRACE and GRACE-FO missions.AOD1B was recently updated to its most recent release AOD1B RL07 and is expected to be used as a background model in the next GRACE and GRACE-FO Level-2 releases.AOD1B RL07 is based on 3-hourly atmospheric data from the ERA5 reanalysis (Hersbach et al. 2023) by the European Centre for Medium-Range Weather Forecasts (ECMWF) as well as simulated ocean bottom pressure (OBP) variations from the MPIOM ocean model (Jungclaus et al. 2013) forced with ERA5 atmospheric data.Even though background models are occasionally updated and thereby improved over time, they will necessarily remain imperfect.As a result, high-frequency signals not removed from the GRACE and GRACE-FO sensor data will lead to residual temporal aliasing artefacts in the monthly solutions.In fact, the errors due to imperfect dealiasing are considered to be among the largest contributors to the overall GRACE and GRACE-FO error (Flechtner et al. 2016).
There are different approaches for mitigating the impact of residual aliasing errors in GRACE data processing.Most notably, several studies have shown that including an estimation of the uncertainty of the background model data can help improve the quality of the gravity field solutions.Zenner et al. (2010) and Kvas et al. (2019) suggested that including uncertainty estimations allows for a weighting of the measurements according to the associated model error.As a result, measurements associated with a larger model uncertainty have a reduced impact on the final gravity field solutions and therefore mitigate some of the effects from residual temporal aliasing.Similarly, employing the uncertainty estimates of ocean tide models has been shown to have a positive impact on the gravity solutions in dedicated performance simulation studies (Abrykosov et al. 2021).
While assessments of the residual errors of AOD1B have been performed in the past (Dobslaw et al. 2016;Poropat et al. 2019), they were only based on AOD1B RL05.This assessment is now believed to be no longer representative for the uncertainties of AOD1B RL07 given the numerous changes made over the last two releases (Shihora et al. 2023b).In this study, we thus derive a new estimation of the non-tidal atmosphere and ocean background model errors associated with AOD1B RL07 that can be readily used in the gravity field estimation process of satellite gravimetry as well as in simulation studies.Our update is especially timely in light of the ongoing efforts towards developing future generations of satellite gravimetry missions, which include both double-pair constellations and novel quantum gravity concepts (Schlaak et al. 2022;Zhou et al. 2023).
We start this work by assessing the signal content in AOD1B RL07 and how the represented variability has changed, especially for the oceanic component in Sect. 2. This is done by comparing the AOD1B update to the ITSG2018 daily gravity field solutions (Kvas et al. 2019;Mayer-Gürr et al. 2018).We then develop an estimation of the uncertainties in the atmospheric component of AOD1B through a comparison of the employed ERA5 reanalysis data to other state-of-the-art atmospheric reanalyses (Sect.3).In Sect.4, we focus on the uncertainties in the oceanic component using ensemble simulations where we quantify both the impact of the atmospheric forcing and the impact of the intrinsic variability.The derivation of a new time-series of true errors representative of the uncertainties within AOD1B RL07 is presented in Sect.5, and the computation of new error variance-covariance matrices for the application in simulation studies is described in Sect.6.The paper concludes with early application examples of the newly derived stochastic information for GRACE-like simulations (Sect.7) and a summary in Sect.8.

Comparing AOD1B to ITSG daily solutions
We start by examining the signal content of AOD1B RL07 in relation to residual signals remaining in previously published GRACE/GRACE-FO gravity field time-series.A comparison of the new release RL07 to its predecessor RL06 shows that the largest differences in variability are found in the oceanic domain.In contrast, the atmospheric differences over the continents are much smaller (Shihora et al. 2023b).This can be expected given the lack of assimilated observations in the ocean simulation.In turn, this also suggests that the uncertainties of AOD1B are going to be dominated by the dynamic contribution of the simulated ocean bottom pressure.To assess the degree to which the oceanic mass variations are not captured in RL07, we make use of a series of daily gravity field solutions provided by the Institute of Geodesy at Graz University of Technology (ITSG).ITSG-GRACE2018 (ITSG2018 in the following) is provided in terms of spherical harmonic coefficients up to degree and order 40 and is based on a combination of GRACE measurements and prior stochastic information in a Kalman smoother framework (Kvas et al. 2019;Mayer-Gürr et al. 2018).These daily solutions incorporated AOD1B RL06 in its processing in conjunction with the previous estimate of the associated AOD1B uncertainties.They thus represent residual mass variations not captured by AOD1B RL06.Given their global coverage and connection to actual GRACE observation, the daily gravity field solutions have already been applied successfully in several oceanic applications (Bonin and Save 2020) and were also utilised in the assessment of differences in high-frequency ocean model simulations (Schindelegger et al. 2021).
For our analyses, we synthesise an equiangular onedegree grid based on spherical harmonic coefficients from daily ITSG2018 solutions for 2004-2006.Similar to the approaches of Eicker et al. (2020) and Schindelegger et al. (2021), we transform a binary land-ocean mask from spherical harmonics onto the same grid and reject all grid points with a value below 0.8 to generate a coastal buffer.As a reference, the resulting standard deviation of the residual OBP signal of the ITSG2018 solutions is shown in Fig. 1 for three frequency bands as obtained from a fourth order Butterworth filter.The variability is shown for 3-10 days (Fig. 1a), 10-30 days (Fig. 1b) and 30-60 days (Fig. 1c).
For the highest frequencies, the residual OBP variations are mainly located in coastal regions as well as in the Southern Ocean in resonant basins and in the band of the Antarctic Circumpolar Current (ACC).Especially in the Bellingshausen Basin, the residual variability reaches values up to 2 hPa, i.e. 2 cm in equivalent water-height.The picture is similar for the moderate frequencies (10-30 days) although the strongest signals are now found south-west of Australia.For the longest frequencies we consider here (30-60 days), the residual OBP variability is generally much weaker, suggesting that the OBP variations at these frequencies are better captured by the AOD1B RL06 background model data which were subtracted during the satellite data processing.As a reference, we also show the standard deviation of AOD1B RL07 in the same frequency bands in Fig. 2. Comparing both figures indicates that for the shortest periods, the residual ITSG variability matches the overall variability of the background model, i.e. the residual circulation signal is proportional to the overall signal.This is especially visible in the Southern Ocean.In other parts, such as the northern part of the Pacific and especially for longer periods, that correspondence is significantly reduced.
Next, we consider the impact of the new release AOD1B RL07 on the residual OBP variations.As the ITSG2018 time-series already considers the AOD1B RL06 background model data, we compare the ITSG signal content only to the update of AOD1B, i.e. the difference (RL07-RL06).We then assess the impact of the model update by computing explained variances using: where AOD1B is the update to the AOD1B background model data through RL07.
The results are shown in Fig. 3 for the same three frequency bands as before.In all three bands, there are clearly regions where the update to AOD1B captures part of the residual circulation signal (red) and regions where the update does not capture the residual variability (blue).Blue areas are especially prevalent for the highest frequencies in the lower latitudes.While the same is true for the 10-30-day and 30-60-day bands, the effect is less pronounced.However, it should be noted that the negative explained variances could also indicate that the ITSG solution does not capture the variability properly, which may particularly hold for highest frequencies (Schindelegger et al. 2021).Regions where the AOD1B update captures the residual variability are in all three cases found in the band of the ACC as well as in the Arctic Ocean for the medium and long periods.Comparing the results to the amount of residual variability presented in Fig. 1, we find that the regions with negative explained variances correspond largely to areas where there is very little residual variability present in the ITSG time-series.This is especially clear for the shortest periods where the ITSG variability around the equator is close to zero.As a result, the explained variances, which are a metric relative to the ITSG variability, likely appear as highly exaggerated.When focusing only on regions with a significant amount of residual circulation signal, however, it turns out that these correspond to the areas with a positive explained variance.Examples are, for instance, the Bellingshausen Basin for the 3 -10-day band, or off the coast of south-western Australia in the 10 -30-day case.So while in many regions, especially lower latitudes, the two datasets do not correspond well, there are some regions such as parts of the Southern Ocean and the Arctic, where the variability is better captured by AOD1B RL07.Those conclusions are also consistent with previous evidence based on satellite altimetry and GRACE-FO alongtrack data reported by Shihora et al. (2022a).
Based on the results presented so far, there are local and regional improvements when considering the new release of AOD1B.However, not all of the residual oceanic mass variations present in the ITSG2018 daily solutions are captured.Hence, the uncertainty of the background model data can be expected to show significant changes compared to the earlier estimation of Dobslaw et al. (2016) which calls for a new error assessment.In the following section, we thus focus on the calculation of a new realisation of true errors for AOD1B RL07 based on model differences from atmospheric reanalyses as well as differences in ocean model simulations for subsequent use in satellite gravity data analysis.
3 Atmospheric surface pressure differences AOD1B considers the non-tidal mass variations from both the oceans and the atmosphere.For RL07, the atmospheric component is based on the ECMWF's ERA5 reanalysis data (Hersbach et al. 2023) until 2017, followed by operational ECMWF data from 2018 onward.While the reanalysis data are constrained through observations, they will still include errors induced by insufficient or conflicting observations as well as insufficient modelling of atmospheric dynamics.These uncertainties are typically distributed globally and depend for example on parametrisations, orography, etc.A common approach to address these is through the comparison of NWM fields published by different institutions.
We therefore compare the ERA5 surface pressure data to two other state-of-the-art atmospheric reanalyses: the Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA2) (Gelaro et al. 2017), from NASA's Global Modeling and Assimilation Office, and the Japanese 55-year Reanalysis (JRA55) (Kobayashi et al. 2015) from the Japan Meteorological Agency.Some characteristics regarding resolutions and data assimilation scheme are given in Table 1.
As all three reanalyses feature a different horizontal resolution, we unify all datasets by remapping to a regular 0.5 • grid following Dobslaw et al. (2016).We further resample to a six-hourly temporal resolution and subtract the mean surface pressure in each case.This is also in line with the resolution of the previous error assessment.To eliminate the impact of low frequencies and high-frequency atmospheric tidal signals which are not part of AOD1B, we subsequently apply a bandpass filter with cut-off periods of 1 and 30 days.In all three cases, the largest differences in surface pressure variations are found in the Southern Ocean and Antarctica.Differences between ERA5 & MERRA2 as well as MERRA2   5a where the surface pressure variability exceeds 10 hPa in high latitudes, differences between the reanalyses reach up to 30% for low latitudes (< 20 • ) and only 10% for higher latitudes (> 20 • ).These small differences between the reanalyses show that surface pressure variations are generally well captured by all three of them.The result is not surprising, given that the reanalyses share essentially the same physics and observational data considered for assimilation.Only in regions where the density of observations is sparse, e.g. in Antarctica, larger differences are found.Based on the differences in the reanalyses presented in this section, we thus choose to base the atmospheric component of the new uncertainty estimation on the differences between ERA5 and MERRA2.

MPIOM ensemble simulations
Next to the residual uncertainties in the atmospheric component, we also consider uncertainties in the ocean domain.In contrast with the atmospheric mass variations provided by NWMs, the OBP variability is not observationally constrained.Instead, it is based on free-running forward simulations with the Max-Planck-Institute for Meteorology Ocean Model (MPIOM) (Jungclaus et al. 2013) forced using atmospheric fields from the ERA5 reanalysis.More details on the configuration of the ocean model are given in Shihora et al. (2022a).Given the lack of observational constraints, the residual uncertainties in AOD1B RL07 are expected to be much larger compared to the atmospheric component as it was already the case for the previous estimation.To get an estimate of the residual uncertainty of the oceanic component of AOD1B, we set up an ensemble simulation using MPIOM.
In particular, we focus on two sources of uncertainty.The first source is based on the differences in the atmospheric reanalyses which will result in differences in the ocean dynamics and consequently also in differences in the OBP variability.In the following, we will refer to the variability induced by the atmosphere as forced variability.As a second contribution, chaotic intrinsic variability can arise through non-linear ocean processes.While they are typically associated with smaller scales, they can map into larger variations through non-linear interactions (Arbic et al. 2012;Zhao et al. 2021).We will refer to these variations as intrinsic variability going forward.

Forced variability
Excluding tides, high-frequency mass variations in the oceans are largely caused by atmospheric surface winds leading to a redistribution of water masses.These wind-driven barotropic changes are particularly pronounced in middle to high latitudes.Gradients in atmospheric surface pressure over the oceans can also drive relevant OBP variations (Ponte 1993).While at low frequencies the ocean surface can be expected to compensate the atmospheric pressure anomalies, at higher frequencies, the response is non-equilibrium (i.e.involves currents and mass motion).While atmospheric reanalyses such as ERA5 capture atmospheric dynamics with a great deal of realism, they are of course not perfect and residual uncertainties remain as shown in the previous section.This includes presented differences in surface pressure but also differences in other relevant atmospheric fields such as surface winds.
To assess the impact of these atmospheric differences on the MPIOM simulations, we perform three model experiments from 1995 to 2020.One is forced using atmospheric ERA5 data.This simulation is thus equivalent to the simulation used in AOD1B RL07 with the only difference being that for the simulation here we have used a 3-hourly forcing frequency in order to be consistent with the other two experiments.The second and third MPIOM runs use either the MERRA2 or the JRA55 reanalysis data for the atmospheric forcing.In both of these cases, we also use 3-hourly forcing.The atmospheric forcing considers contributions from atmospheric pressure, near-surface horizontal wind speed and stresses, solar radiation, precipitation, cloud cover, temperature and dew point temperature.We extract bottom pressure fields from the ocean model, subtract the temporal mean and bandpass filter the results using 1-and 30-day cut-off periods.
As a reference, we show the OBP variability, i.e. standard deviation, of the ERA5 forced run alone in Fig. 5b.Regions with the highest variability exceeding 5 hPa are found in the Southern Ocean in the region of the ACC, especially in the Bellingshausen Basin and the South-Australian Basin.Additionally, high variability is found in shelf areas as well as the Arctic Ocean where OBP variations are largely driven by barotropic pressure changes (Bingham and Hughes 2008).
We now turn to the differences between the simulations with varied forcing.Figure 6 shows the standard deviation of OBP differences between the ERA5 and MERRA2 runs (a), the ERA5 and JRA55 runs (b), and the MERRA2 and JRA55 runs (c).In all three cases, the largest differences in variability match the regions that show a high variability in the first place as shown in Fig. 5b.The largest signals are found again in the Southern Ocean, where the wind-driven barotropic variability is generally high, but also differences in the atmospheric reanalysis data are largest (see Fig. 4).In these regions, the OBP differences reach values of up to 1 hPa which amounts to 15 -20% of the variability is those regions and even over 50% east of the Drake Passage.This highlights the sensitivity of the high-frequency mass variations in the Southern Ocean to the atmospheric forcing.Comparing the individual subfigures reveals that the differences between the ERA5 and JRA55 forced runs are smaller than the difference of either to the MERRA2 forced simulation.Again, this is consistent with the differences in atmospheric surface pressure 4 where the differences between ERA5 and JRA55 tend to be the smallest.

Intrinsic variability
Next we turn to the estimation of the contribution due to intrinsic variability in MPIOM.Oceanic intrinsic variations emerge not through variations in the atmospheric forcing but instead emerge from mesoscale turbulence on scales of O(10−100) km and O(10−100) days (Sérazin et al. (2015) and references therein).Impact of this initial intrinsic variability is also found at much larger and longer (i.e.interannual) scales (Penduff et al. 2011) suggesting a spontaneous inverse cascade to these scales (Sérazin et al. 2018).Studies show that this intrinsic variability has a significant impact Fig. 5 Standard deviation of ERA5 surface pressure (SP) and simulated ocean bottom pressure (OBP) from MPIOM using atmospheric ERA5 forcing data.Results are bandpass filtered using a 4th-order Butterworth filter to only contain periods in the 1-30-day band Fig. 6 Standard deviation of ocean bottom pressure (OBP) differences from three different MPIOM simulations where the atmospheric forcing is varied.Subfigures show the results for differences between ERA5 vs. MERRA2 atmospheric forcing (a), ERA5 vs. JRA55 forcing (b) as well as MERRA2 vs. JRA55 forcing (c).All OBP fields are bandpass filtered using a 4th-order Butterworth filter to only contain periods in the 1-30-day band on a number of oceanic variables such as sea-level (Sérazin et al. 2015), ocean heat content (Penduff et al. 2019), transports (Cravatte et al. 2021) or OBP (Zhao et al. 2021) on various time-scales.In order to disentangle the impact of the atmospheric forcing from the intrinsic variations, a common approach is to perform parallel ocean simulations with identical forcing that only differ in some small initial perturbations.Differences in the resulting variability can then be used to examine intrinsic variations (Sérazin et al. 2015).We employ the same approach in the following.
We perform three MPIOM simulations that are all based on ERA5 atmospheric forcing data.All simulations start in the year 1994 from a transient ERA5 run started in 1960 based on a 2000 year long spin-up simulation with daily climatological forcing.One simulation, the reference, uses the "correct" initial conditions for the year 1994.For the second simulation, we start the run in 1994 but employ the initial conditions based on the year 1993, thus shifting the initial state by 1 year.The third run uses the initial conditions after only 1000 years of the spin-up.This represents a larger difference in the initial conditions for comparison.We then compute differences in 3-hourly OBP between all three simulations and band-pass filter the results to only include frequencies in the 1-30 day range.The results in terms of standard deviations of OBP differences are given in Fig. 7. Figure 7a shows the difference between the reference simulation and the simulation with the initial conditions shifted by 1 year, whereas Fig. 7b shows the difference between the reference and the simulation using the initial conditions after half the spin-up.Figure 7c shows the difference between the two simulations using shifted initial conditions.We explicitly note here the difference in scale, which is in this case given in Pa, compared to previous figures.
The impact of the initial conditions varies strongly with latitude.Differences are largest in the Southern Ocean along the band of the ACC and reach values of up to 30 Pa.In contrast with the impact of the atmospheric forcing, the differences are thus smaller and amount only to 30-50% of the forced variability in the Southern Ocean.Other regions show a very limited sensitivity to the choice of initial conditions.Most of the open ocean and also coastal areas are largely unaffected.Comparing the three standard deviation differences with each other shows that the spatial distribution and magnitude are very similar suggesting that the distance in time by which the initial conditions are shifted is not relevant here.As the MPIOM configuration used here and in AOD1B Fig. 7 Standard deviation of ocean bottom pressure (OBP) differences from three different MPIOM simulations where the initial conditions are varied.Subfigures show the results for differences between using the correct initial conditions for 1994 vs. a 1-year shift in initial condi-tions (a), 1994 vs. initial conditions after half of the spin-up simulation (i.e.1000 years of spin-up) (b) as well as a one-year shift vs. using the ocean state after half of the spin-up (c).All OBP fields are bandpass filtered to only contain periods in the 1-30-day band is not eddy-permitting, there is likely little mesoscale activity and consequently also a reduced amount of intrinsic variability (Penduff et al. 2019).This might explain why the impact of the initial conditions is comparatively small compared to the forced variability and compared to results presented by Zhao et al. (2021).

AOe07 time-series
In the previous sections, we have analysed and compared atmospheric mass variations from different state-of-the-art reanalyses, which gives an estimate of the residual atmospheric uncertainty.Similarly we have examined the oceanic uncertainty in the forced variability and intrinsic chaotic variability in MPIOM through ensemble simulations.Next, we combine the individual components to derive a single timeseries of true errors representative for AOD1B RL07 that can be used in the processing of either satellite gravimetry or in simulation studies for MAGIC.
The final time-series should merge the atmospheric information over the continents with the oceanic uncertainties elsewhere.For the atmosphere, we choose to use the differences between the ERA5 and the MERRA2 reanalyses.This combination is compatible with the ERA5 based AOD1B RL07 but also shows the larger differences in atmospheric mass variations and is thus deemed a better estimate of the residual uncertainties.For the ocean component, we combined the impact of the forced variability and the intrinsic variability in OBP.This is done through a new MPIOM simulation forced with atmospheric MERRA2 data that is also based on a shift in the initial conditions by one year.The uncertainty information over the oceans is then given by the difference between the new simulation to the reference simulation based on ERA5 data.Atmospheric tides in the surface pressure data and atmospherically induced tides in the simulated ocean bottom pressure are estimated and subtracted in the same way as for AOD1B RL07 and as described by Balidakis et al. (2022).The atmospheric and oceanic components are subsequently combined into a single time-series of 6-hourly fields and then highpass filtered using a cut-off frequency of 30 days in order to only represent the high frequencies relevant for satellite gravimetry.
There is, however, an issue that needs to be addressed especially for the ocean component.As the uncertainty estimation is based on differences between OGCM simulations with MPIOM only, the uncertainties are likely underestimated in their magnitude.As indicated by Quinn and Ponte (2011), model differences as they are used here tend to underestimate the high-frequency OBP variations, especially when they are based on simulations using the same ocean model.Similarly, the impact of the intrinsic variability greatly depends on the resolution of the ocean model (Sérazin et al. 2015).We here use, just as in AOD1B RL07, MPIOM's TP10L40 configuration based on a 1-degree tri-polar grid and thus expect the intrinsic variations to be under-represented in the ensemble configuration.While the spatio-temporal pattern of OBP variability in the Southern Ocean is captured, its magnitude and certain peculiarities (e.g. the Argentine Gyre) are likely underestimated.Similar to the approach of the previous uncertainty estimation (Dobslaw et al. 2016), we thus calculate a scaling factor for the oceanic component.This is done through a comparison of the oceanic uncertainty estimate to the residual OBP variations in the ITSG2018 daily solutions after RL07 has been subtracted.We then select a global scaling factor to adjust the uncertainty time-series up to the ITSG variability without exceeding it.Note that the global factor of 2.4 is applied to the oceanic component of the time-series only.For the atmospheric contribution over the continents, the strong dependence on the assimilated barometer data means that no scaling is required.We note, however, that this procedure does not introduce any actual GRACE or GRACE-FO gravity information from the ITSG solutions into the error time-series.This approach merely compares the overall magnitude in variability.
The thus derived final error time-series, labelled AOe07, can then be considered an estimation of the residual uncertainties in the AOD1B RL07 background model.In Fig. 8, we show the standard deviation of the new AOe07 time-series (a).In addition, Fig. 8b shows the previous uncertainty estimation AOerr, which was developed by Dobslaw et al. (2016) to represent the model deficiencies of AOD1B RL05.AOerr was based on operational and ERA-Interim atmospheric data from the ECWF as well as simulated OBP from the OMCT ocean model.The AOerr time-series is available for the 12year period from 1995 to 2006 as part of the ESA ESM.
Globally, the new time-series features a smaller variability, reflecting the improvements made to the AOD1B background model data over the years.The reduction is both visible over the oceans and over the continents.The discrepancy in the accuracy of the background model data between the atmospheric component over land and the simulated ocean bottom pressure is still reflected in the new time-series.
The largest uncertainties are found in the coastal areas where the amplitude of the barotropic high-frequency variations is also comparatively high.Additionally, increased uncertainties are also found in the Southern Ocean, as described in Sect.4.1.Compared to the previous release, the magnitude is smaller in the Southern Ocean although the regions with enhanced uncertainties remain similar.In contrast with the previous estimation, however, the uncertainties in the Arctic Ocean are significantly reduced which can be attributed to MPIOMs much better representation of the Arctic Ocean based on a tri-polar grid.
Technically, AOe07 is available as a 6-hourly series of fully normalised Stokes coefficients from a spherical harmonic expansion up to degree and order (d/o) 180.The time-series covers 26 years from 1995 to 2020 and thus extends the previous version by 14 years.The data can be accessed, together with the previous uncertainty assessment, via the ESA ESM repository under ftp://ig2-dmz.gfzpotsdam.de/ESAESM/.

Variance-covariance matrices
One way to include the error estimation of background models is through the use of a variance-covariance matrix (VCM) which represents the spatio-temporal uncertainties.Such a VCM can then be used either in the gravity field estimation process or in dedicated simulation studies as shown by Abrykosov et al. (2021) for the case of ocean tides.In the latter case, the approach offers an opportunity to signifi-cantly improve the gravity field retrieval performance if the non-tidal background stochastic modelling is improved as well.In this section, we derive a new VCM based on the updated AOe07 time-series that captures the spatio-temporal uncertainties of the non-tidal atmosphere and ocean highfrequency mass variations.
The calculation of the VCM is based on the computation of both variances as well as covariances between the Stokes coefficients via where X l 1 ,m 1 stands for the C and S Stokes coefficients of degree l 1 and order m 1 while X represents the temporal mean value.Variances are computed analogously using l 1 = l 2 and m 1 = m 2 .Based on Eq. 1, a fully populated stationary VCM is calculated up to d/o 40.This then matches the resolution of current GRACE daily solutions such as ITSG2018.Additionally, a second diagonal matrix, containing thus only variances, is calculated up to the full d/o 180.Like the AOe07 timeseries, both matrices are publicly available under Shihora et al. (2023a).We note that the VCMs computed in this way do not include any regularisation as it is implemented, e.g. for ocean tides (Abrykosov et al. 2021).In contrast with the tidal case, the AO VCM is based on a much larger number of epochs (i.e. 26 years of 6-hourly data) and in preliminary tests does not pose any problems in the application.Nonetheless, it should be kept in mind that the VCM is not regularised.In the following, we will turn to possible applications in GRACE-like and Next-Generation-Gravity-Mission (NGGM) simulations.

Application in simulation studies
In the context of satellite gravimetry, the background model uncertainties (as given by AOe07 and the VCM) can be applied in a variety of ways.Specifically, within the GRACE and GRACE-FO data processing, it can be applied in the gravity field estimation process to weight observations according to the uncertainty of the associated background model information.Secondly, AOe07 can be used in dedicated satellite gravimetry simulation studies.Especially the estimation of gravity field retrieval errors is regularly performed in the context of future gravity mission scenarios and their comparison.Here we present an example application of AOe07 in a Mass-Change and Geosciences International Constellation (MAGIC) simulation scenario that consists of a polar pair at 488 km and an inclined pair at 397 km.As a preparatory step, we perform a GRACE-FO baseline simulation for the year 2002 using GFZs EPOS-OC software (Zhu et al. 2004).The time-variable source model for the simulation is based on the ESA Earth-System-Model (ESA ESM) which also provides the so called DEAL coefficients (Dobslaw et al. 2016).The DEAL coefficients contain the unperturbed de-aliasing model that is applied as a background model in case that perfect model-based AO predictions are assumed.To arrive at a realistically perturbed background model for the simulations, the DEAL coefficients are either perturbed using the old AOerr product, which then results in a perturbed background model that corresponds to AOD1B RL05, or the new AOe07 series, which is designed to represent the capabilities of AOD1B RL07.
In the simulation, we perform two gravity field recoveries using either DEAL + AOerr or DEAL + AOe07 as a background model and subsequently subtract the HIS component (i.e.terrestrially stored water, continental ice-sheets and solid Earth) of the ESA ESM to derive gravity field retrieval errors.We note that the retrieval errors obtained in this way contain both the AO-aliasing error and the HIS aliasing.Alternatively, the gravity field recoveries can be performed using ether DEAL + AOerr + HIS or DEAL + AOe07 + HIS, which then excludes any HIS aliasing errors in the final retrieval errors.
Figure 9 illustrates the gravity field retrieval error for the new AOe07 error time-series (red, dashed) in comparison with the previous AOerr time-series (blue, solid).The mean HIS signal is given in black as indicator for the actual geophysical signal that needs to be captured by a satellite mission.Thin lines represent monthly results for 2002, while the bold lines indicate the yearly mean results.Subfigure (a) shows the results including only the AO-aliasing error, while (b) includes the AO-aliasing contribution as well as the HISaliasing.
The simulation results suggest an overall reduction in the retrieval error of a few mm EWH, which we consider as relevant.For the mean results shown here, the RMSE is reduced by about 30 mm.While some reduction is also seen in parts at low degrees, the reduction is best visible for degrees above about d/o 30, and as a result, the mean HIS signal is recoverable to slightly higher spatial resolution.While there are slight variations from month to month, the results are fairly consistent for all computed months of the year 2002, indicating that when AO errors are considered, the simulated gravity retrieval error is reduced using the new AOe07 error estimation.Most likely, this can be attributed to the overall reduced amplitude of the new error time-series compared to its predecessor.Comparing the two subfigures shows that the impact of the HIS-aliasing is much smaller than the AOcontribution.This also underlines why we specifically focus on the AO contribution in the error estimation.It also corroborates the decision of the GRACE project team to abstain from the application of a hydrological de-aliasing model for Level-2 gravity field processing.

Summary and conclusions
For this study, we have produced and examined a new quantification of the remaining uncertainties in the latest release 07 of the non-tidal Atmosphere and Ocean De-Aliasing Level-1B (AOD1B) product.The newly updated uncertainty estimation is called AOe07 and can be used both for the processing of monthly gravity field solutions for mitigating the impact of residual temporal aliasing, and for dedicated simulation studies of future satellite gravimetry mission concepts.
We have shown that the new release of AOD1B does indeed improve the representation of the oceanic highfrequency variability by comparing both RL07 and RL06 to the daily ITSG2018 gravity field solutions in regions where ITSG displays a significant amount of variability.Considering that the represented variability of AOD1B is less accurate over the oceans compared to the atmospheric contribution over land (Shihora et al. 2023b), an update of the residual uncertainty estimation was deemed necessary.
The new estimation is based on model inter-comparisons both for the atmosphere and the ocean.For the atmospheric part, we consider model differences between the latest ECMWF reanalysis ERA5, which is the basis of AOD1B As, however, the oceanic component is based on model differences from the same ocean model, the uncertainties are most likely underestimated.Similar to the previous release by Dobslaw et al. (2016), this issue has been addressed through a time-invariant global scaling factor that is applied to the entire oceanic domain.The scaling factor is determined based on a comparison of the model differences to the residual variability in the daily ITSG data.While this is not an ideal approach, it allows the magnitude of the variability to be adequately captured.However, we note that the ratio of intrinsic versus forced variability is smaller than in comparable studies (Zhao et al. 2021) and the contribution of the intrinsic variability is under-represented in certain regions.We attribute this to the limited spatial resolution of the MPIOM configuration, which does not resolve mesoscale processes which are the ultimate driver of the intrinsic variability.We believe that, given the ongoing progress towards the next generation of satellite gravimetry missions and the associated studies, the timely availability of applicable background model uncertainties is important, so that AOe07 has been made available nonetheless.In a possible future refinement, additional investigations that specifically address the impact of the ocean model resolution as well as exploring differences based on different ocean models including high-resolution shallow-water codes should be considered along with a larger number of model runs to better capture the effects of the intrinsic variability.In that way, an even more realistic uncertainty estimate could be provided.
Based on the AOe07 time-series, we also provide a variance-covariance matrix that represents explicitly the (time-averaged) spatial correlations for further use in simulation studies.A successful example of the application of VCMs is given in Abrykosov et al. (2021) for the case of ocean tides.Lastly, we have demonstrated the application of the new data-set in an exemplary satellite gravimetry simulation for the future MAGIC constellation.The results, which compare the gravity field retrieval error when using either the new AOe07 data or the previous error estimation, show a general improvement in the monthly retrieval by about 30 mm EWH, especially at higher degrees, and thus demonstrate a better signal recovery by a few degrees.We recommend that simulation studies are based on the ESA ESM DEAL coefficients in combination with AOe07 to arrive at a realistically perturbed de-aliasing model that is compatible with recent developments in background models.

Fig. 3
Fig.3Amount of variance explained by the update of AOD1B, i.e. the difference (RL07-RL06) using the residual circulation signal from ITSG2018 daily gravity field solutions as a base time-series.Results

Fig. 8
Fig. 8 Standard deviation of the new time-series of true errors adapted for AOD1B RL07 (a) and the previous error estimation (b).In both cases, 12 years from 1995 to 2006 are used

Fig. 9
Fig. 9 Simulated degree amplitudes of gravity retrieval errors in mm EWH considering only the AO noise contribution.Results using the previous AOerr (new AOe07) estimation are given in blue (red).Thin