Determining and disentangling the spatial and temporal variability of a field from sparse in-situ measurements is a long-standing problem in space physics. Addressing this problem, along with providing broader coverage in measurements, is commonly raised as a major justification for designing and launching multi-platform satellite missions (e.g., Escoubet et al., 2001; Burch et al., 2016). The Geospace Dynamics Constellation (GDC) mission, that is currently in the formulation phase, is one of the latest of such missions. GDC’s goal is to significantly increase our understanding of how the coupled ionosphere-thermosphere (IT) system reacts to external energy input from the Sun and the magnetosphere. Specifically, it will aim to address the following overarching science goals: 1) Understand how the high latitude ionosphere-thermosphere system responds to variable solar wind/magnetosphere forcing; and 2) Understand how internal processes in the global ionosphere-thermosphere system redistribute mass, momentum, and energy. Successfully addressing these objectives is challenging as the IT system is known to vary locally and globally, with scale-sizes ranging from sub-kilometer to thousands of kilometers and on timescales that range from seconds or minutes to hours or days.

Completely characterizing the dynamics of the IT system and determining the physical processes underlying its variability would ideally require a very dense (in altitude, longitude, and latitude) sampling of various plasma and neutral parameters over a wide range of spatial and temporal scales. As a practical matter, the number of spacecraft in such a constellation that can provide accurate comprehensive measurements is finite—for example, GDC will fly six spacecraft. This immediately introduces important challenges. On the one hand, relatively fine spatial and temporal resolutions, when compared to the physical scales, are required to properly capture the behavior of the system locally; and on the other hand, the global nature of the GDC mission would require a broad coverage of observing platforms in local time and longitude. Satisfying both conditions with a limited number of in-situ platforms can only be possible via a flexible constellation that evolves to observe the full range of relevant scales. A flexible satellite architecture, in turn, requires sophisticated algorithms that are themselves flexible and optimal in combining observations from individual points.

In this paper, we introduce a software package ‘LSGC-AS-LEO 1.0’, which is an implementation of the least-squares gradient calculation approach developed by De Keyser et al. (2007) and De Keyser (2008) with the focus on its application for the GDC mission. The gradients characterize variations of the field in space and time, helping to disentangle the spatial and temporal ambiguity in measurements in an inherently quantitative manner. In the following section, we briefly describe GDC’s expected constellation architecture and its evolution during different phases of the mission. In Section 3, we present a brief overview of the gradient calculation technique. In Section 4, we evaluate the performance of the technique on simulated GDC measurements of the IT parameters, and discuss its significance for the GDC mission. Section 5 includes the summary, conclusions, and a description of the future steps.

GDC’s constellation consists of six spacecraft in a circular orbit at the altitude of 350–400 km with an inclination of about 81–82°. At a given time, the spacecraft will sample plasma and neutral parameters on a relatively thin spherical shell, the altitude of which will slowly vary between 350 and 400 km due to orbital drag and subsequent reboosts via onboard propulsion. The satellites will be deployed via a single launch vehicle with small inclination differences. Over time, differential precession will lead to separation of the six orbital planes in local time, leading to an ever-widening instantaneous longitudinal coverage. Different phases of the mission can then be defined that are appropriate for the investigation of the IT processes over local (referred to as phase 1, with scale-lengths < 2 h in local time), regional (phase 2, with cross-track (primarily longitudinal at the equator) scale-lengths between 2 and 9 LT hours), and global (phase 3, towards the end of the mission with scale-lengths > 9 h in local time) scales.

Figure 1 shows the expected evolution of the satellite configuration during days 92 (top), 190 (middle), and 764 (bottom) from the start of the science phase of the mission. These correspond to the beginning of phase 1, end of phase 1, and sometime in the late regional/early global phase of the mission, respectively. Shown in each panel are instantaneous positions of all six spacecraft (shown as colored circles) for a scenario where the constellation is in the northern hemisphere and traveling northward. Thin vertical and horizontal lines show 5-degree increments in longitude and latitude. Dashed lines in colors matching each circle show the near future trajectory of that GDC spacecraft, demonstrating near-latitudinal direction around the equator, and predominantly longitudinal direction above ∼75° geographic latitude. Figure 1 clearly demonstrates the dynamic nature of the constellation in space and time and the complexities that one may face when attempting to determine spatio-temporal variability of the ionosphere and thermosphere from GDC measurements. For example, the hexagonal shape of the constellation seen in Panel A1 near the equator may be suitable to simultaneously determine local variability in latitude, longitude, and time; while at higher latitudes (see Panel A6), where the trajectories of the individual spacecraft cross, little information can be obtained in latitude due to lack of broad latitudinal sampling. In this case, the form of the constellation which briefly approximates a ‘pearls on a string’ configuration would be better suited for investigating variability in longitude and time.

The variability in constellation architecture along with the need to monitor its ability to estimate gradients in both latitude and longitude at a given time motivated the creation of a quantity called the ‘quality’ (or ‘Q’) value and the associated ‘Q-vector’. The Q-value can be defined as Q val = 1 + A A e , in which A is the area of the spherical polygon formed by the satellite constellation at a given time after being projected onto a 2D surface using a Mercator projection, and A _e is the area of an equilateral polygon with the same perimeter. A maximum value of Q _val = 2 indicates a ‘regularly’ distributed constellation in space, while lower values imply that the sampling constellation is anisotropically distributed, which may result in less-than-optimal spread simultaneously in both latitude and longitude—for example, a minimum value of Q _val = 1 denotes a set of co-linear spacecraft. Since Q _val does not provide information on the ‘direction’ of anisotropy associated with a constellation, it can be complimented by a second quantity Q dir = tan − 1 Δ lat Δ long which quantifies the direction of the skewness. Q _dir is defined as the inverse tangent of the ratio of the constellation’s latitudinal spread to its longitudinal spread. For a constellation consisting of at least two spacecraft Q _dir is always defined and ranges from 0° to 90°. Q _dir values of 0° and 90° indicate that the constellation only samples along the longitudinal and latitudinal directions, respectively. Using both components of the Q-vector, i,e., Q _val and Q _dir , one can interpret the evolution of the constellation over time and its ability to provide sampling in different spatial directions. In the context of this work, Q-vector is used to monitor the constellation architecture and quantify the latitudinal and longitudinal spread of the spacecraft array. This provides insight into the array’s ability to compute gradients in those specific directions. It is worth mentioning that a more mathematically rigorous approach to evaluate the effective constellation dimensions based on error considerations and through investigation of the eigenvalues and eigenvectors of the spacecraft position tensor was developed by Vogt et al. (2020). There, in the general 3-dimensional case, the eigenvalues give insight into the amount of variance in each of the principal directions of the spacecraft array, which are captured through the eigenvectors. The eigenvectors with the largest eigenvalues correspond to the directions most heavily sampled by the spacecraft array. The approach developed by Vogt et al. (2020) provides a more formal description of the constellation and will be adopted in the updated versions of the LSGC-AS-LEO. It is important to note that the assumption that an isotropic constellation is ideal for sampling is only strictly true when the gradients are isotropic. When this is not the case, ideal constellation configurations may have Q-values appropriate to the ratio of longitudinal to latitudinal gradients in the measured fields.

Figure 2 shows the trajectories of GDC spacecraft as a function of latitude and longitude (top panels) for three orbits corresponding to the beginnings of phases 1a (left), 2a (middle), and 3a (right). For each case, the magnitude and angle of the Q-vector is shown in the bottom panels as the constellation evolves in time. The drop in the magnitude of the Q-value is evident at high latitudes, consistent with the anisotropic constellation shape shown in Figure 1. As will be shown in Section 4, the ability of LSGC-AS-LEO 1.0 to fully determine both longitudinal and latitudinal variability in the IT system using GDC measurements may reduce for smaller Q-values below a threshold. Q-vector can, thus, be used as an empirical measure to optimize constellation design and provide an indication of where robust gradient calculations may be achievable.

Figure 1 and Figure 2 demonstrate the variability of GDC’s constellation geometry as a function of latitude and mission phase. As a result of this variability, any approach that would aim to utilize multi-point measurements from GDC to characterize the dynamics of the IT system would ideally need to be general and flexible, without many assumptions on the sampling scheme, and at the same time ensure the optimum utilization of data from all satellites. Further, a successful approach would need to take into account spatial and temporal scales over which ionospheric and thermospheric structures evolve relative to the size of the constellation, and provide measures of reliability for its calculations. The least-squares gradient calculation with adaptive scaling (LSGC-AS method, De Keyser et al., 2007; De Keyser, 2008) is such an approach and is briefly described in the next section.

Optimally combining information from multi-point measurements to quantitatively determine the local dynamics of a measured field often proves to be challenging. This topic has been the subject of many previous works that have let to the development of different algorithms and tools (e.g., Chanteur, 1998; Harvey, 1998; Paschmann and Daly, 1998; Robert et al., 1998; Amm and Viljanen, 1999; Darrouzet et al., 2006; Denton et al., 2020; Dunlop and Lühr, 2020; Fiori, 2020; Torbert et al., 2020; Bard and Dorelli, 2021; Denton et al., 2022; Zhu et al., 2022). Mathematically, the problem can be translated to one involving calculation of gradients—once the gradients of a field are known with respect to space and time its local dynamics can be better understood. Precise calculations of gradients, however, are generally difficult due to a number of reasons, the most important of which are the presence of noise in data and under-sampled variations of the measured field in between the spacecraft. On the one hand, for estimated gradients to be accurate the spacecraft should form a close constellation to fully resolve variations of the sampled field; and on the other hand, for the estimated values to be meaningful the difference between sampled values should be large compared to the measurement noise—a condition which may require a larger separation between the spacecraft. It is balancing between the two constraints that necessitates the use of sophisticated algorithms for computing gradients. In addition, the spatial scales of variation in the sampled field may vary depending on the process being studied, or the location of the measurements, which introduces further complexity and a need for the algorithm to assess its own performance. Given the complexity of the task of calculating gradients, a significant amount of effort has gone towards developing appropriate algorithms and tools. Much of these works have been performed in support of the European Space Agency’s (ESA) three-satellite Swarm mission (Friis-Christensen et al., 2008). For example, Ritter and Lühr (2006) and Ritter et al. (2013) were among the first to address current estimation in the context of a multi-satellite LEO mission. Vogt et al. (2009) and Vogt et al. (2013) provided tools and an error analysis framework for planar spacecraft configurations with special consideration of physical constraints to estimate out-of-plane contributions. Shen et al. (2012a) and Shen et al. (2012b) developed magnetic gradient estimation techniques with particular relevance to measurements in an environment dominated by field-aligned currents. Blagau and Vogt (2019) and Blagau and Vogt (2023) implemented relevant techniques in a publicly available software package based on Python. Of the mentioned works, those described by Shen et al. (2012b) and Vogt et al. (2020) stand out as general approaches to gradient calculations, both of which utilize some adaptation of the position tensor of the satellite configuration to compute spatial gradients in a least-squares sense. Very recently, in the past 2 years, a series of publications have emerged that introduce further methods to combine information from multiple in-situ platforms. These recent works include: Shen and Dunlop (2023), who introduced a geometrical method based on integral theorems to compute the linear gradients; Shen et al. (2021b), who developed a novel method to estimating both the linear and quadratic gradients using multiple spacecraft observation; Shen et al. (2021a) and Dunlop et al. (2021), who developed algorithms specifically for determining nonlinear magnetic gradients; and the book by Dunlop et al. (2021), that covers a number of multispacecraft techniques that were initially developed to analyze data from the Cluster mission (Escoubet et al., 2001). As mentioned before, here we implement the approach described by De Keyser et al. (2007) and De Keyser (2008), which is also a general gradient calculation tool, developed based on the concept of homogeneity scales in orthogonal directions, including in time. The method aims to ensure the utilization of the maximum amount of information available and allows to assess the meaningfulness of the result by providing error estimates. This approach is briefly described below. The description only includes high-level intuition behind the approach and is not meant as a substitute for the detailed mathematical description in De Keyser et al. (2007) and De Keyser (2008).

Imagine a scalar field f sampled at x _i = [x _i ; y _i ; z _i ; t _i ], for i = 1⋯N, by an arbitrary number of satellites that form a flexible constellation. At a given time, f may be approximated in the vicinity of point x ₀ = [x ₀; y ₀; z ₀; t ₀] via its Taylor series expansion: f ( x ) = f 0 + Δ x ⊤ g 0 + 1 2 Δ x ⊤ H 0 Δ x + ⋯ , where Δx = x − x ₀, f ₀ = f( x ₀), and g ₀ and H ₀ are the gradient and the Hessian at x ₀. If we assume that f changes linearly in the vicinity of x ₀, thus ignoring the terms of the Taylor expansion beyond the second term, we will reach the approximate equation f _a ( x ) = f ₀ + Δx ^⊤ g ₀. Using f _a ( x ), the field can be estimated at each sampled point x _i . Defining residues r _i = f _a ( x _i ) − f _i at each point, where f _i is the measured value at x _i , and requiring the residues to be zero leads to a system of N equations for the five unknowns—with the unknowns being f ₀ and four partial derivatives with respect to x, y, z, and t. Often times, gradients would be determined using data obtained in a period of time over which the field is continuously sampled via multiple in-situ platforms. For example, calculating gradients in a 30-s window where six satellites obtain measurements every three seconds would lead to 60 samples. In this case, N = 60 ≫ 5 and the system is highly over determined. One may proceed to select five samples to uniquely determine the unknowns, but the outcome would likely depend on the choices as well as the amount of measurement noise in the system. Furthermore, this approach would not utilize the maximum amount of information available, nor would it provide any insight on the reliability of the estimates. A more suitable approach would be to utilize all the available data points and solve the over-determined system of equations using a weighted least squares minimization procedure. In this case, equation i, corresponding to sample x _i , would be weighted inversely proportional to the total amount of error that is expected for the residue r _i . Once the weights are known, the system of equations may be written in the matrix form and solved via singular value decomposition. Error bars on the estimated gradients as well as the ill-conditioning degree of the problem can be obtained from an appropriate covariance matrix and singular values—for more details see De Keyser et al. (2007) and De Keyser (2008).

Much of the complexity of the weighted least-squares minimization approach mentioned above is condensed in obtaining proper weights, w _i , for each equation based on the expected error associated with the corresponding measured sample. The choice of weights as inverse of the total error is to ensure that measurements with large errors do not contribute much to the solution of the over-determined system of equations. In the absence of systematic errors, the error δf _i in estimating sample f _i via f _a ( x _i ) may include three components. One is the ‘measurement error’ which is the inherent uncertainty associated with the instrument that provides f _i ; the second is the ‘approximation error’ (or the ‘curvature error’) which exists due to ignoring the nonlinear terms of the Taylor expansion in f _a ( x ); and the third is the error due to the potential presence of unresolved small-scale fluctuations in the field. As described by De Keyser et al. (2007), the small-scale fluctuation error can be incorporated in a statistical sense, the measurement error can be modeled via a zero-mean random noise with a standard deviation specific to each instrument, and the approximation error can be modeled as δ f a ∝ f c | x i − x 0 | β 2 , where f _c is a constant related to the Hessian of the field at point x ₀ and | x _i − x ₀| _β is the normalized distance between the measurement point x _i and x ₀—such that the modeled approximation error increases with distance according to the third term of the Taylor expansion. Here, | x _i − x ₀| _β is not merely the physical distance between two points but one that takes into account the linearity of the field along each orthogonal direction (including time). In other words, | x i − x 0 | β = Σ k ( x i k − x 0 k ) 2 l k 2 , where the subscript k indicates components along an orthogonal direction k and l _k are homogeneity lengths defined below. To clarify this point imagine a scenario where the field f varies linearly over short scales along axis x but over large scales along axis y. In this case, approximating f along y via the first two terms of the Taylor expansion remains accurate over a longer distance when compared to x. As such, the approximation error would need to be scaled differently along x and y. To incorporate such conditions, De Keyser et al. (2007) utilized the concepts of ‘homogeneity lengths’. The homogeneity length, l _k , is the distance along a direction k over which a linear variation is a good approximation to a function f. The geometry defined by all homogeneity lengths thus defines the homogeneity space in which the field can be reasonably modeled by the first two terms of the Taylor expansion with relatively small errors in the form of δ f a ∝ f c | x i − x 0 | β 2 . Beyond the homogeneity space the modeled approximation error would need to increased rapidly with | x _i − x ₀| _β such that the information obtained outside of the homogeneity space contribute minimally to the computation of gradients. This can be enforced via incorporating an additional term in the approximation error.

The homogeneity scales are properties of the physical field and change in space and time. In general, their determination is not trivial. De Keyser (2008) has described an approach to automatically determine the homogeneity scales from the measurements. The approach could be best understood by noting that solving the weighted over-determined system is equivalent to minimizing χ ², where: χ 2 = 1 N ∑ i = 1 N r i 2 δ f i 2

However, one should bear in mind that in the presence of measurement and curvature errors, minimizing χ ² may not be a suitable target. In reality, the computations should be performed such that the residue r _i reflect the correct amount of expected error associated with the calculations at point x _i . In other words, one would reasonably expect that r i 2 ≈ δ f i 2 , and thus χ ² ≈ 1 (for N ≫ 1). In fact, it can be shown mathematically that under certain assumptions χ ² follows the χ ²-statistics, under which the most likely value for χ ² is ∼1.

The process of identifying the homogeneity scales, and consequently solving the over-determined system of equations, then starts from a set of initial guesses for the homogeneity scales, followed by the calculation of the total expected error for each measurement point, assigning weights to each equation, computing the gradients, determining χ ², then updating the homogeneity scales based on the value of χ ². This cycle may be incorporated in an optimization scheme to allow χ ² → 1 after a number of iterations. Several optimization approaches are introduced in De Keyser (2008). We have currently implemented the ‘common rescaling’ approach, which will be further discussed in the context of Figure 5. Future work will involve implementing the other approaches developed by De Keyser (2008) and exploring additional ones.

Motivated by the upcoming GDC mission, and with the goal of evaluating the ability of the satellite constellation to resolve local dynamics of the IT system, we have implemented the approach in an open-source software LSGC-AS-LEO 1.0, presently written in MATLAB. Despite the immediate focus of the software on GDC, LSGC-AS-LEO 1.0 is implemented in a general form, utilizing many settings that can be adjusted by the user, and as such it maintains the ability to support other missions with minimal updates to few functions. Here, we present and discuss typical outputs of the software on synthetic and simulated fields relevant to the IT system. Evaluating the results using simple synthetic test cases allows to validate the implementation, examine the robustness of the underlying algorithms, and further provides the opportunity to clarify certain important concepts of the algorithm described in Section 3 and in De Keyser et al. (2007) and De Keyser (2008). Additional test cases using realistic global simulations of the IT system along with GDC’s proposed ephemeris data further allows showcasing the significance of the software for GDC’s future measurements. Because GDC’s measurements will occur in a narrow altitude region between 350 and 400 km, our implementation of the LSGC-AS algorithm is tailored to assume spherical coordinates, and to assume zero radial variation, i.e., all measurement points are sampled at the same altitude and all the gradient calculations are performed in the spherical coordinates. Extending this procedure to account for non-zero radial variations is straightforward.

For the first test case, imagine a synthetic thermospheric neutral temperature of the form T n = 1000 + 200 ⁡ sin 1.5 θ − 40 + 200 ⁡ sin 2 ϕ + 0.1 t 100 2.1 where θ, ϕ, and t are latitude, longitude, and time in units of degree, degree, and seconds, respectively. T _n is not meant to represent the actual neutral temperature in the thermosphere; instead it serves as a simple function with simultaneous variability in latitude, longitude, and time. Running the software with an appropriate label specifying the synthetic function results in the following actions: GDC’s predicted ephemeris data over a predefined period of time are imported and used to sample the synthetic T _n with an adjustable cadence; a predefined level of measurement noise is then introduced to the sampled values; datapoints [ x _i , f _i ] from all or selected spacecraft gathered within a period of time centered at t _i are then passed onto a routine in order to determine the gradients at point x _i . The gradients are calculated at an arbitrary ‘evaluation point’, which by default is defined at the geometric center of the constellation. An optimization routine is used that adaptively modifies the homogeneity scales such that computations result in an appropriate outcome for χ ². Various settings that define the details of the computations are specified prior to running the routine. These include the initial guesses for homogeneity scales, f _c , the measurement noise level, available satellites, the phase and time of the GDC mission, the evaluation point, optimization method, measurement cadence, integration time, whether calculations are performed on f or log(f), etc. The final result is an estimate of the value of T _n at the evaluation point, the longitudinal and latitudinal gradients of T _n at the evaluation point, the temporal rate of change of T _n at the evaluation point, and error bars for these quantities, as well as estimates for the homogeneity lengths.

Figure 3 shows the results of the gradient calculations for the synthetic function T _n sampled during the first 100 min of GDC’s phase 1a. The input parameters to LSGC-AS-LEO that have been used to produce this and the following plots are summarized in Table 1. In the top panel of Figure 3, T _n , as seen by six individual satellites, are shown with thin black lines. The size of the constellation in phase 1a is small compared to the scales over which T _n varies in space and time, and as a result the six spacecraft measure relatively similar values—i.e., the black lines nearly overlay each other—with the exceptions at the highest latitudes near times 14 and 60 min. Also shown in Panel A are the true (red) and the estimated (blue) values of T _n at the evaluation point, here defined as the center of the constellation. The calculated latitudinal, longitudinal, and temporal gradients along with their expected errors (i.e., their 3 − σ error bars) are shown in Panels B–D in blue. The calculated gradients closely follow the true gradients which are shown in red, with the exceptions at times ∼14 and 60 min where ∇ _lat T _n briefly deviates from the red curve and its error bars significantly increase. Noting the latitude and longitude of the evaluation point and the components of the Q-vector in Panels E and F, respectively, it becomes clear that the difficulty in determining ∇ _lat T _n (the latitudinal vector component of the gradient) around times ∼14 and 60 min is due to the shape of the constellation as the spacecraft approach the polar regions. As was shown in Panel A6 of Figure 1, at these times (primarily poleward of 75° geographic latitude in either hemisphere) the constellation spreads along longitude without much coverage in latitude, leading to poorly-constrained estimates for ∇ _lat T _n . This longitudinal elongation of the constellation is also reflected in Panel G which shows the condition number defined in De Keyser et al. (2007) as c o n d = s min 2 / s max 2 , where s _min and s _min are the minimum and maximum singular values obtained during the singular value decomposition of the over-determined system of equations. The condition number serves as a general measure of ill-conditioning of the problem and decreases by about an order of magnitude as the constellation elongates at higher latitudes, signaling large errors in parts of the calculations. Finally, the bottom two panels of Figure 3 demonstrate the performance of the implemented common rescaling optimization routine described in De Keyser (2008). The purpose of the optimization routine is to adaptively rescale the initial homogeneity lengths such that χ ² approaches a predefined value, resulting in refinements to the error estimates as well as an estimate of the scale lengths over which the calculated gradient is expected to be accurate. In this case, in order for the error bars in Panels A–D to represent the 3 − σ errors, optimization works to enforce χ ² ∼ 1/9. Panels H and I show the calculated latitudinal homogeneity length and χ ², respectively. While the initial homogeneity scales in latitude, longitude, and time has been chosen as L _lat = 60°, L _long = 160°, and L _t = 30 min, respectively, the optimization routine scales all the homogeneity lengths by a common factor that varies between ∼7–50 as the constellation moves in space and sample different features of T _n . Most notably, as the constellation approaches the polar region between times 10–20 and 55–65 min the homogeneity scales decrease, consistent with larger variations in ∇ _long T _n .

It will be insightful to repeat the gradient calculation test described above with slight modifications. In the left panels of Figure 4 we show the results from a similar run in which the measurement error introduced to T _n is now significantly increased from 0.1% to 3% of f _i . This is an important test as GDC will include various particle and field instruments with different levels of measurement accuracy for which the performance of the gradient calculation need to be tested. As can be seen in Panels A–D, an immediate consequence of the elevated noise level is larger errors for the calculated gradients, which are also reflected by the increased error bars. This is not surprising since gradient calculations are highly sensitive to the presence of background noise. Further, it is seen that, unlike in Figure 3B, the incorrect latitudinal gradients between 13–19 and 59–62 min are no longer properly covered by the error bars—this is perhaps due to inaccurate choices of initial homogeneity scales or f _c . Nevertheless, the decreased condition number at these times can still serve as a cautionary flag regarding the accuracy of the results. In the right panels of Figure 4 we show the results from a similar run with yet another modification: the time axis now corresponds to the first 100 min of GDC’s phase 2b. In this phase the constellation size (i.e., the average distance between spacecraft) increases. The individual spacecraft now sample different values of T _n with differences that are greater than the 3% background noise level—see Panel J. Consequently, the accuracy of the estimated gradients increases and their error bars significantly drop in Panels K–M. This is due to the reduced noise amplification nature of the gradient operator over longer distances, and applies to cases where the measurement error dominates over the approximation error.

Before wrapping up the evaluation of the gradient calculation approach on the synthetic field T _n , it would be helpful to revisit the impact of the optimization and the adaptive scaling of the homogeneity lengths from a different prospective. Consider the results shown in Figure 4J–R at the single time t ₀ = 74.75 min. Here, the constellation approached the polar regions with the latitude of the evaluation point at 80.1°. At this time, about 120 samples [ x _i , f _i ], collected by all six spacecraft within 2.5 min from t ₀, have been used to calculate the gradients shown in panels K–M. According to Panel H, the optimization has decreased the homogeneity lengths by a common factor of about 7, from the initial scales of L _lat = 60°, L _long = 160°, and L _t = 30 min. In Figure 5 we show the distribution of the residues |r _i | for each datapoint x _i as a function of the normalized distance | x _i − x ₀| _β . The absolute value of the residues are shown for the following two cases: the black datapoints correspond to the case where the least-squares minimization problem is solved without the optimization procedure (thus, only according to the initial homogeneity lengths); whereas the blue datapoints correspond to the case where the optimization has been turned on (thus allowing the homogeneity lengths to be rescaled). For the case of adaptive scaling turned off, the least-squares minimization approach is successful in solving the over-determined system of equations such that the residues are small and on the order of the measurement noise added to f _i (shown by the horizontal black bar). Despite the reasonably small residues, however, we note that the gradients calculated using the initial homogeneity scales do not correctly represent the actual values—in other words, the calculated gradients are legitimate, yet inappropriate, solutions of the over-determined system. With adaptive scaling enabled, the homogeneity scales decrease, | x _i − x ₀| increases, and as a result only a fraction of the datapoint that are closest to evaluation point are effectively used to calculate the gradients. The calculated gradients in latitude, longitude, and time are then modified by −96, 97, and −740%, respectively, compared to the previous case which brings them very close to the true gradients. At the same time, we observe that the distribution of |r _i | as a function of the normalized distance | x _i − x ₀| follows the following expected form: from the definition of f _a ( x ) in Section 3, one would generally expect | r i | = | f a ( x i ) − f i | ≈ δ f m , i + f c | x i − x 0 | β 2 , where δf _m,i is the measurement error at point x _i . Accordingly, in Figure 5 for smaller values of | x _i − x ₀| the residues are small, on the order of the measurement noise, while at larger distances |r _i | increases rapidly due to larger approximation errors.

With the implementation of the least-squares gradient calculation approach validated via the synthetic example, we now turn to more realistic scenarios to test the ability of the software to determine the local dynamic of the ionosphere and thermosphere using GDC measurements. Figure 6A shows a frozen-in-time snapshot of the electron density as a function of geographic latitude and longitude from a three-dimensional TIEGCM simulation, featuring the Equatorial Ionization Anomaly (EIA). EIA is a dayside/dusk ionospheric feature that is formed from a phenomenon known as the Fountain/Appleton effect (Appleton, 1946). This phenomenon is a result of a vertical ion drift near the magnetic equator that is related to the E-region dynamo (Hanson and Moffett, 1966). This results in a concentration of plasma at ± 20° in magnetic latitude with a depletion of plasma at the magnetic equator primarily in the F-region (200–450 km) (Rishbeth, 2000). The EIA is in a highly collisional domain and is coupled with the thermosphere, but due to the lack of spatial and temporal resolutions from data collected in previous missions, there are unanswered questions to the exact coupling mechanisms in this region. Here, we employ simulations of the EIA using NCAR’s High-Resolution Thermosphere-Ionosphere-Electrodynamics General Circulation Model (HR-TIEGCM). HR-TIEGCM is a global 3D numerical model that can simulate the IT region from ∼97 km to ∼600 km with a longitude/latitude grid of 0.625 ° × 0.625 °, a vertical resolution of 0.25 of a scale height, and a timestep of 60 s (Qian et al., 2014; Dang et al., 2021). This temporal and spatial resolution is on par with what the GDC mission will be capable of capturing.

Our goal in this example is to evaluate the performance of LSGC-AS-LEO 1.0 to calculate the spatial gradients associated with the EIA using GDC data. In order to properly display the calculated gradients in the context of the two-dimensional image on Panel A, we briefly ignore the temporal variations of the electron density and assume that temporal gradients are zero at all locations. Upon running the software, the simulated data are imported and interpolated every 20 s along the trajectories of the GDC spacecraft during the first few hours of phase 1a. 4% measurement noise is then introduced to the data, which is comparable to the expected performance (i.e., 2% accuracy, and 1% precision) of GDC’s Atmospheric Electrodynamics probe for THERmal plasma (AETHER) instrument. The gradients are calculated with the initial homogeneity scales of L _lat = 10°, L _long = 50°, and L _t = 2,500 s. Overplotted on the electron density map on Panel A is the trajectory of the evaluation point (the continuous black curve), as well as the spatial gradient vectors (black arrows). One can verify that the gradient vectors always point nearly perpendicular to the electron density contours. Shown on the right panels, in a similar format as those in Figure 3, are the electron densities seen by individual satellites (Panel B), the latitudinal (C), longitudinal (D), and temporal (E) gradients, latitude and longitude of the evaluation point (F) at the center of the constellation, L _lat (G), and χ ² (H). While in Panel A gradient vectors are calculated and shown for about 5 GDC orbits through the simulated field, the right panels only show time series for the first 120 min (∼1.3 orbits) from the beginning of the interval starting at round −30° latitude and −106° longitude. The blue circles overplotted on the trajectory of the evaluation point in Panel A mark temporal timesteps of 10 min, allowing to compare the time series plots with the two-dimensional map in Panel A. The gradients are calculated accurately with the exception of the temporal gradients which suffer from very large error bars. The reason for these error estimates will be discussed later when visiting Figure 7. The spatial gradients shown in Panels C and D allow us to visualize the underlying electron density field without referring to the two-dimensional map on Panel A. For example, from the large negative longitudinal gradients around time 40 min, one could conclude that the electron density peak seen in Panel B at that time is due to a primarily longitudinal structure. This can indeed be confirmed by visiting Panel A and noting that observed electron density enhancement around time 40 min is due the constellation traversing along the eastern edge of the EIA at the longitude of ∼60° and the latitude range of −20 to +40°.

In Figure 7 we repeat a similar test as that in Figure 6 while focusing on the constellation architecture at later phases of the mission. Here, the left panels correspond to the first 100 min of phase 2a and the panels on the right correspond to the beginning of the phase 2b. The figures show electron densities observed by individual satellites (Panel A and F), the calculated gradients (B–D and G–I), and the latitude and longitude of the evaluation point (F and J). In both examples we allow the simulated electron density field to evolve in time, such that the true temporal gradient are non-zero. At the beginning of phase 2a, it appears that the spatial gradients associated with the EIA are often properly captured. However, by the beginning of phase 2b, the size of the constellation has increased such the certain smaller scale features are no longer resolved, and at times the sharpest gradients are underestimated—e.g., around time 140–150 and 185–195 min. As a consequence of these errors, the estimated electron density at the evaluation point, at the center of the constellation, also deviate from the true values in Panel F. With respect to the temporal gradients, in both Panels D and J, the estimated values are significantly greater than true gradients and are accompanied by very large error bars, indicating that the estimates are likely not meaningful.

The difficulty in determining ∇ _t N _e in this case arises from several factors which include: 1) the constellation architecture during the chosen phases of the mission, including the large (∼7 km/s) velocity of the satellites, that limits the duration of time over which a given spatial region is visited; 2) the ratio of the spatial gradients to ∇ _t N _e at a given point; and 3) the relatively large measurement noise introduced to the samples. For example, consider the largest true temporal gradients in Figure 7D that reach a maximum of ∼ 4 × 1 0 7 m − 3 / s at time 64 min. If persisted during a 2-min interval over which the gradients are computed δN _e due to ∇ _t N _e may reach to ∼ 5 × 1 0 9 m − 3 . This is an order of magnitude smaller than the measurement noise introduced at these times, and two orders of magnitude smaller than δN _e due to the latitudinal gradients. Nevertheless, as will be shown in Figure 8, temporal gradients of IT parameters can still be reasonably obtained using LSGC-AS-LEO 1.0 and GDC data depending on a number of factor, including the phase of the mission, measurement noise level, and the form of the variations of the measured field. It should also be emphasized that in low-measurement-noise conditions ( < 1%), applicable to high-precision electric and magnetic field data and their products such as E × B ion drifts, or plasma densities obtained from the background plasma frequency, even small temporal gradients can be reliably captured by the current approach. Finally, it should be noted that, as can be seen in the predicted GDC ephemerides (Rev C data files: https://ccmc.gsfc.nasa.gov/RoR_WWW/GDC_support/Proposer_Resources/GDC_EphemerisRevC.zip; description of data files: https://ccmc.gsfc.nasa.gov/static/files/GDC-SCI-DESC-0005RevC.pdf), there are other sampling configurations (e.g., ‘Follow the Leader’ configurations in Phases 2 and 3) in which GDC will have longer revisit times, allowing the study of slow temporal variations.

In Figure 8 we return to the constellation geometry in phase 1a and evaluate the performance of the software in disentangling the spatial and temporal variability in a vector quantity: the zonal (U, shown on the left) and meridional (V, shown on the right) neutral winds as simulated by TIEGCM. Here, in order to demonstrate the ability of the approach to capture temporal gradients only 1% noise has been added to the simulated values. For reference, the expected performance of GDC’s neutral wind measurements by the Modular Spectrometer for Atmosphere and Ionosphere Characterization (MoSAIC) instrument is the accuracy and precision of no more than 4.5 m/s. With this level of measurement noise, all the gradients are determined with reasonably small error bars, including periods of non-zero ∇ _t U _n and ∇ _t V _n between times ∼120–150, as well as the most prominent variations of spatial gradients near the polar regions. While the specified 4.5 m/s accuracy of neutral wind measurement by the MoSAIC instrument translates to measurement noise levels as low as 2.25% in our example, an effective measurement noise of 1% or below is achievable by averaging multiple datapoints from MoSAIC that will be obtained within the 20-s measurement cadence used in our example shown in Figure 8. MoSAIC will provide independent samples of ion and neutral densities, temperatures, compositions, and velocities every two seconds or faster, providing a minimum of 10 samples within 20 s, providing the possibility to effectively reduce random measurement noise level by a factor of ∼3, to below 1%.

It is important to recognize that the ability of a mission like GDC to measure gradients relies on a reasonably good match between the spatiotemporal separation of the satellites and the phenomena being studied. Early phases of the mission are best for studying small-scale, sharp variations, and later phases are best for studying larger-scale, more gradual variations. As an example consider the simulated neutral temperature field shown in Figure 9A. The image shows variation of T _n in space at a fixed time, with the minimum and maximum values of ∼850 and 1200 K, respectively. For reference, the trajectory of the GDC constellation during the first orbit of phase 2b are also overplotted on the image. In the right panels of Figure 9, the gradients computed during phase 2b of the mission are shown to successfully capture the spatial variations of the underlying field. It should also be emphasized that, while the examples shown in Figures 3–9 demonstrate that the least-squares gradients calculation approach is a powerful tool capable of estimating the gradients of the IT parameters using GDC’s future measurements, the utilization of the approach is not necessarily trivial or free from limitations and caveats. For example, the performance of the algorithm currently relies on appropriate choices for several parameters, such as the homogeneity lengths, the determination of which may prove to be difficult in a real scenario. Further, for the approach to be applicable in a routine and robust manner, the sensitivity of the chosen parameter values needs to be studied. Such investigations will be the subject of future works.

In preparation for the Geospace Dynamics Constellation mission, we have implemented the least-square gradient calculation approach described by De Keyser et al. (2007) and De Keyser (2008) in a open-source software LSGC-AS-LEO 1.0. Calculation of the gradients of the ionospheric and thermospheric parameters would allow to disentangle the spatial and temporal variability of the measured fields in an inherently quantitative manner. The approach that is implemented is flexible and, thus, highly suitable for GDC’s dynamic constellation architecture. Via a number of examples, and by utilizing synthetic and simulated IT fields, we have validated the performance of the software and demonstrated its power to explore GDC’s ability to resolve spatial and temporal structures that may exist in the ionosphere and the thermosphere. We have also established the usefulness of the software to evaluate suitability of various constellation geometries and assess the impact of measurement sensitivities on addressing GDC’s science objectives. For example, from the test cases discussed one may draw the following several conclusions: 1) computation of the temporal gradients of neutral and plasma variables, while are sensitive to the measurement noise level, are possible with GDC measurements; 2) The spatial gradients of the equatorial ionization anomaly can be reasonably resolved during phase 1 of the mission, while at the later phases the gradients are likely to be underestimated. On the other hand, in the presence of measurement noise, computing the gradients of the neutral temperature would likely be more difficult in the earliest phases of the mission due to the small gradients and large homogeneity lengths involved; 3) Gradients of the neutral wind can be well determined in the earliest phases of the mission even at the highest latitudes where the constellation skews in longitude.

We will continue to enhance the software in several directions. The future work will include further refinements to the optimization scheme and implementation of direction-dependent scaling for the homogeneity lengths, implementation of constrains (e.g., divergence free or curl free) in the algorithm for vector fields, enhancement of computation efficiency, and performing additional tests on simulated data to further probe GDC’s ability to address its science questions.