## 1 Introduction  
Since 2002, the Gravity Recovery and Climate Experiment (GRACE) satellite mission has provided the scientific community with invaluable gravimetric data [1, 2]. By measuring Earth's monthly time-variable gravity field from 2002 to 2017, it allowed us to infer natural fluctuations and anthropogenic fingerprints on the Earth's global water cycle. This includes areas such as determining changes in continental hydrology and ocean bottom pressure [3, 4]. It allowed us to characterize the global melting of glaciers and ice sheets through the measurement of their mass loss (\(\sim\)600 Gt/year) [5]--the rate at which ocean mass is contributing to sea-level rise [6]. The societal and scientific value of mass transport data obtained from GRACE has been recognized in such a manner that a successor to GRACE-FO was recommended as one of five designated missions within the decade following GRACE (2017-2027) in the National Academy of Sciences Decadal Survey for Earth Science and Applications from Space [7]. The successive mission, GRACE Follow-On (GRACE-FO) [8], has continued the monthly mapping of Earth's time-variable gravity field since 2018. Due to GRACE and GRACE-FO's single polar pair mission architecture, full global coverage for a single pair can take up to 30 days for spatial resolutions of a few hundred kilometers, thus leading to a spatial and temporal resolution of the detected signals of \(\sim\)300 km\({}^{2}\) and 30 days, respectively. This means that signals at higher spatial and temporal resolutions, such as rapid flooding events, cannot be quantified [9]. These mission architectures are limited in their spatiotemporal performance. They are affected by temporal aliasing errors, which are caused due to undersampling of submonthly geophysical signals, corrupting the gravity solution. A higher frequency sampling of the gravity field would be a straightforward method to reduce such errors. Recently, there has been a request by the community for data products with increased spatiotemporal resolution [10]. Past mission analyses have tackled this problem by comparing the improvement of a two-pair GRACE-like (i.e., low-low satellite-satellite tracking) constellation to a single-pair configuration [11; 12; 13; 14].  
An \(n_{pair}\) GRACE-like satellite constellation is a possible solution to reach higher spatiotemporal resolutions and resolve submonthly gravity field signals. Due to technological advances in small satellite technology, many scientific missions can be fulfilled with small satellite hardware [15; 16; 17; 18; 19]. Several groups [20; 21] have been exploring the possibility of collecting GRACE-type intersatellite ranging measurements from smallsats or cubesats, which would dramatically decrease the costs of the satellites. This is a critical path to pursue, as the implementation of new gravity-sensing technologies (atom interferometry, etc.) are still some years away [22]. A constellation composed of small satellites could allow for higher spatial and temporal resolution as well as multiple decades of observations since it lacks single-point failures that are present in single pair satellite constellations. Should one satellite or satellite pair stop functioning, it would not disrupt the rest of the constellation to carry on the mission, albeit potentially with limited capability [23; 24; 25; 26; 27]. The relatively large number of satellites in the constellation would additionally drive the cost down by enabling the production of instrumentation in larger numbers [29; 30]. Creating the appropriate tools and verifying them by performing corresponding simulations in a high-fidelity environment is necessary to optimize such an \(n_{pair}\), GRACE-like mission design. This paper aims to design a methodology that, given a set of constraints, optimizes the orbital configuration of a GRACE-like constellation for the recovery of temporal gravity field variations.