## 4 Simulation Setup  
To build a constellation of GRACE-type satellites, we make use of the two main simulation tools developed by NASA's Goddard Space Flight Center (GSFC)--GEODYN [47] and SOLVE [48]. GEODYN is GSFC's state-of-the-art geodetic parameter estimation and precision orbit determination system that implements an iterative weighted least squares estimation algorithm to solve for spherical harmonic coefficients. Once several satellite pairs have been simulated with GEODYN, they can be combined using SOLVE, a large linear systems solver that combines data of various days and configurations into one gravity field solution. To create realistic simulations, the models described in Table 3 are used for all forthcoming simulations. These models, truncated at degree and order 60, have subdaily temporal resolutions. The ECMWF and GLDAS models are provided at a 3 h temporal resolution where NCEP, OMCT, MOG2D, and the ESA ice model are defined with a 6 h temporal resolution. They represent the dominant mass variations for Earth. \begin{table}
\begin{tabular}{c c c} \hline
**Parameter** & **Value** & **Unit** \\ \hline Propagation time (T) & 29 & (day) \\ Satellite pairs (\(n_{pairs}\)) & 6 & (pairs) \\ Orbit altitude (\(h_{alt}\)) & \(500\pm 5\) & (km) \\ Repeat period (RP) & 29 & (day) \\ Eccentricity (e) & 0 & (-) \\ Arg. of perigee (\(\omega\)) & 90 & (\({}^{\circ}\)) \\ Intersatellite distance (\(r_{int}\)) & 100 & (km) \\ Spatial cell (SC) & 4551 & (-) \\ Temporal cell (TC) & 16 & (per-day) \\ \hline Inclination (_i_) & [0–180] & (\({}^{\circ}\)) \\ RAAN (\(\Omega\)) & [0–360] & (\({}^{\circ}\)) \\ Mean anomaly (\(M\)) & [0–360] & (\({}^{\circ}\)) \\ \hline \end{tabular}
\end{table} 
Table 2: Fixed and variable parameters used in the optimization process. For the static gravity field, both the truth and nominal case have been selected to be equal since the static gravity field does not have any high-frequency mass variations and does not contribute to the temporal aliasing error. Different ocean tides, and atmospheric and ocean models have been chosen for both the truth and nominal cases. The difference in the ocean tide models is proportional to the tide models' error magnitude. In contrast, the difference between the atmospheric and ocean model represents the errors due to the atmosphere and ocean models. The instrumentation on-board is assumed to be a microwave-ranging instrument, as is the case with the GRACE mission. This analysis considers satellite error sources, such as intersatellite ranging measurement, accelerometer, and GPS position errors. The satellite-to-satellite ranging measurement noise is described by a K-band microwave ranging system characterized by a white noise spectral density of \(1.8~{}\upmu\mathrm{m}/\sqrt{\mathrm{Hz}}\)[12]. The spectral density of accelerometer errors in the radial and transversal directions are \((1+0.005/f)^{1/2}\times 10^{-10}~{}\mathrm{m}/\mathrm{s}^{2}/\mathrm{Hz}^{1/2}\), and for the normal axis, is \((1+0.1/f)^{1/2}\times 10^{-9}~{}\mathrm{m}/\mathrm{s}^{2}/\mathrm{Hz}^{1/2}\)[12]. The position errors due to GPS measurements are characterized by white noise with a standard deviation of 1 cm added to each of the three positions' axis [58]. A drag-compensated system will be assumed. The assumptions on instrument performance are based on identified intersatellite ranging and accelerometer technologies for SmallSat technologies from the Mass Change Designated Observable Study Team [59]. These simulations aim to recover the continental hydrology and ice mass variations. The simulations are performed in a series of steps. Given the initial conditions of a satellite pair, these are propagated using the defined set of truth models. The resulting output is the truth satellite-to-satellite range-rate data and satellite position. The previously defined range-rate measurement noise and accelerometer noise are superpositioned to the truth range-rate data. Errors to each of the three positions' axis of the satellites is added.