Secular variation signals in magnetic field gradient tensor elements derived from satellite-based geomagnetic virtual observatories

We present local time-series of the magnetic field gradient tensor elements at satellite altitude derived using a Geomagnetic Virtual Observatory (GVO) approach. Gradient element timeseries are computed in 4-monthly bins on an approximately equal-area distributed worldwide network. This enables global investigations of spatio temporal variations in the gradient tensor elements. Series are derived from data collected by the Swarm and CHAMP satellite missions, using vector field measurements and their along-track and east-west differences, when available. We find evidence for a regional secular variation impulse (jerk) event in 2017 in the first time derivative of the gradient tensor elements. This event is located at low latitudes in the Pacific region. It has a similar profile and amplitude regardless of the adopted data selection criteria and is well fit by an internal potential field. Spherical harmonic models of the internal magnetic field built from the GVO gradient series show lower scatter in near-zonal harmonics compared with models built using standard GVO vector field series. The GVO gradient element series are an effective means of compressing the spatio temporal information gathered by low-Earth orbit satellites on geomagnetic field variations, which may prove useful for core flow inversions and in geodynamo data assimilation studies.


INTRODUCTION
The main part of the geomagnetic field is generated in the Earth's fluid outer core by a process known as the geodynamo. Knowledge of how this core field varies with space and time provides information on the fluid flow dynamics in the liquid metal outer core. Although the temporal behaviour of the geomagnetic field is well characterized in time series from ground observatories, a global spatio-temporal study is hampered by the uneven distribution of these observatories. Even though low-Earth-orbit (LEO) satellites do provide good global coverage on timescales of weeks and longer, the direct study of the first time derivative of the core field, the secular variation (SV), from satellite measurements is not straightforward. This is because LEO satellites are not geostationary (i.e. do not have the same orbital period as the Earth's rotation, such that their position does not stay fixed as seen from ground), which leads to an ambiguity, since it is not possible to establish whether an observed field variation is caused by a temporal or spatial change of the field (Olsen & Stolle 2012) . Spherical harmonic (SH) field models derived from satellite measurements provide an established way of studying the SV field and its time derivative, the secular acceleration (SA), globally. However, such harmonic functions have global support, which means that an estimate of the SV at a specific position may be affected by noise from remote locations.
These issues lead  to introduce the concept of Geomagnetic Virtual Observatories (GVOs) in space, in which satellite magnetic measurements from within a selected region, collected during one month time windows, were used to derive a local monthly mean vector field at the satellite mean altitude. The resulting GVO time series resemble monthly mean series computed using ground observatory magnetic measurements, by providing the magnetic vector field elements at fixed locations. However, since they are based on satellite data, regular sampling in both space and time is possible. Thus, the GVO method provides a tool for compressing satellite magnetic field measurements into a dataset that contains time series distributed in a global grid, and which may also be supplied with error estimates (Hammer et al. 2021c). Olsen & Mandea (2007) used CHAMP measurements to derive GVO vector component time series, and carried out a global investigation of SV that identified a regional geomagnetic jerk event in 2003.
In the original GVO approach of , processing of the satellite measurements followed that of simple monthly field means at ground observatories, taking measurements from all local times and with all levels of geomagnetic activity, and relied on the assumption that short period Secular Variation Signals in Magnetic Field Gradient Tensor Elements 3 external fields would have zero mean over the course of one month. However, later studies revealed that external fields, especially those due to the magnetospheric ring current and ionospheric current systems, cause contamination of the retrieved internal GVO field signal (Olsen & Mandea 2007;Beggan et al. 2009;Domingos et al. 2019). In addition, insufficient local time sampling from within one month of polar orbiting satellites resulted in a bias due to the local time dependence of ionospheric and magnetosphere-ionosphere coupling currents (Shore 2013). Recently, the GVO processing algorithm has been further developed in an effort to reduce contamination from magnetospheric and ionospheric sources, and the local time sampling bias, with the aim of better isolating the field signal generated by the Earth's outer core (Hammer et al. 2021a). These new GVO vector field series have been used to study global patterns of field changes (Hammer et al. 2021a,c), for inferring fluid flows close to the core surface (Kloss & Finlay 2019;Rogers et al. 2019) and for data assimilation studies (Barrois et al. 2018;Huder et al. 2020).
In parallel to the development of these GVO-based techniques there has also been recent progress in the theory of space-based magnetic gradiometry, inspired by advances in satellite gravimetry. Initial studies have demonstrated that knowledge of the second-order 3 × 3 magnetic gradient tensor may be beneficial when seeking to retrieve small scale features of the field (both the lithospheric field and the time-dependent core field). This is possible because gradient elements effectively give more weight to shorter wavelengths, while at the same time some noise sources (e.g. unmodeled magnetospheric fields) are predominantly of long wavelength, which can result in a higher signal-to-noise ratio for short wavelengths compared to using the vector field components (Kotsiaros & Olsen 2012.
Information on these shorter spatial length scales is crucial for core dynamic studies, yet their robust determination is a major obstacle in core field modelling. This issue has motivated us to derive time series of field gradient elements using the GVO method, and to investigate whether such series allow, via for example spherical harmonic analysis, an improved retrieval of field changes. Gradient GVO series have the potential to also be used directly in core flow inversions and data assimilation studies of the geodynamo.
Assuming a potential field due to internal sources and no in-situ electrical currents, the field becomes a solenoidal irrotational vector field and the gradient tensor has the special property of being symmetric with a trace of zero. The assumption of a symmetric gradient tensor reduces the number of independent gradient tensor elements from nine to six, while a trace of zero further reduces this number to five. Each element of the magnetic gradient tensor may be considered as a directional filter providing specific information on the magnetic field structures. Thereby, certain gradient tensor elements better constrain specific spherical harmonics (Olsen & Kotsiaros 2011). According to the studies of Kotsiaros & Olsen (2012) and Kotsiaros & Olsen (2014), knowledge of the radial gradient of the radial field, written as [∇B] rr , is particularly suitable for resolving the higher degree parts and zonal harmonics. The East-West gradient of the azimuthal field, [∇B] φφ , and radial field, [∇B] rφ , are especially sensitive towards sectorial harmonics, while the North-South gradient of the radial, [∇B] rθ , and meridional, [∇B] θθ , fields are especially useful for determining near-zonal harmonics. The East-West gradient of the meridional field, [∇B] θφ does not provide significant additional information. In addition, knowledge of how certain external fields may influence certain gradient tensor elements is important to consider, for instance the magnetospheric ring current is expected to affect zonal terms constrained by the [∇B] rr element but not the [∇B] rφ element (Kotsiaros & Olsen 2014). Although it is not yet possible to directly measure the full magnetic gradient tensor in space (Nogueira et al. 2015), it is nonetheless possible to compute the tensor elements from a magnetic potential determined using satellite magnetic measurements.
In this paper, our aim is to estimate local time series of the magnetic field gradient tensor elements using the GVO method. Our primary motivation for deriving such gradient field series, is to improve the recovery of the small length-scales of the field's secular variation, compared with the traditional use of vector field data. In particular, the gradient elements are expected to enable a higher signal-tonoise ratio as compared to vector elements (Kotsiaros & Olsen 2014). This is because they are more sensitive to the small length scales of the field, and less affected by large scale external fields. We use satellite magnetic measurements to derive the GVO gradient series, and follow Hammer et al. (2021a) in implementing dark and quiet time data selection criteria and use 4-month time windows to minimize problems related to local time sampling.
In Section 2 we provide a detailed description of the satellite magnetic measurements and selection criteria used, and in Section 3 we describe the GVO method and the computation of both GVO vector and GVO gradient element time series. In Section 4.1 we present results of the GVO series for each of the SV gradient elements, and visually inspect these. In order to investigate the possible benefits of using GVO gradient data, we compare SH field models derived epoch by epoch from the GVO vector data and GVO gradient tensor data in Section 4.2. We study the detailed behaviour of the gradient tensor elements going from 2015 to 2018 with a focus on the Pacific region. Finally, Section 5 provides discussions and conclusions.

DATA
To derive the GVO time series we select vector magnetic field measurements from the CHAMP and Swarm satellite missions. We used CHAMP L3 magnetic data between July 2000 and September 2010 and Swarm Level 1b MAG-L, version 0505/0506, from all three Swarm satellites Alpha, Bravo and Charlie between January 2014 and April 2020, and sub-sample at 15s intervals (i.e. tak-Secular Variation Signals in Magnetic Field Gradient Tensor Elements 5 ing measurements every 15 second), in the vector field magnetometer (VFM) frame. Next, the magnetic data in the VFM frame are rotated into an Earth-Centered Earth-Fixed (ECEF) local Cartesian North-East-Centre (NEC) coordinate frame (for details see Olsen et al. (2006)) by using the Euler rotation angles from the CHAOS-7.2 model . Measurements from known problematic days (e.g. where satellite manoeuvres took place) were removed and gross data outliers for which the vector field components deviated more than 500nT from CHAOS-7.2 field model [http: //www.spacecenter.dk/files/magnetic-models/CHAOS-7/], see also Finlay et al. (2020) values were rejected. The measurements were then selected using a dark quiet time criteria defined here as: a) the sun is at least 10 • below horizon, b) geomagnetic activity index K p < 3 • , c) ring current index |dRC/dt| < 3nThr −1 , merging electric field (averaged over two-hours) at magnetopause E m ≤ 0.8mVm −1 , and placing constraints on the interplanetary magnetic field (IMF) (averaged over two-hours) requiring B z > 0nT and |B y | < 10nT (Ritter et al. 2004). Using the definition in Olsen et al. (2014), the merging electric field was derived using 1-min values of the solar wind, solar clock angle and IMF extracted from the OMNI database, http://omniweb.gsfc.nasa.gov.
Previous studies have demonstrated the benefits of using along-track differences of the satellite magnetic field measurements for retrieving higher spatial resolution of the core and also the lithospheric fields, as such differences filter out correlated noise caused by external sources Kotsiaros et al. 2015;Kotsiaros 2016;Finlay 2019). However, from data differences alone it is not possible to obtain robust information on the longer wavelength part of the field. Therefore, we also use the complementary data means in a similar manner as in previous studies of the core signal (Sabaka et al. 2013;Hammer 2018). From the satellite magnetic field measurements, B k (r), where k is the vector component of a given coordinate system, we use measurement means, Σd k , and differences, ∆d k as data. The differences, ∆d k = (∆d AT k , ∆d EW k ), and the means, Σd k = (Σd AT k , Σd EW k ), are taken along-track (AT) for each satellite and East-West (EW) between the Swarm Alpha (SWA) and Charlie (SWC) satellites. Here along-track differences are calculated from the 15 s differences ∆d AT k = [B k (r, t)−B k (r+δr, t+15s)] while the means are given by Σd AT k = [B k (r, t)+B k (r+δr, t+15s)]/2. The East-West differences were calculated as ∆d EW k = [B SWA k (r 1 , t 1 )−B SWC k (r 2 , t 2 )], and the means (r 2 , t 2 )]/2. Considering a given orbit of Swarm Alpha, the corresponding Swarm Charlie measurement were chosen to be that closest in colatitude provided that |∆t| = |t 1 − t 2 | < 50s .

THEORY AND METHOD
The purpose of this paper is to derive time series of the magnetic field gradient tensor elements using the GVO method. In section 3.1, we begin by recalling the GVO method and how this is used to derive time series of the magnetic field vector elements. Following this, in section 3.2, we describe our new extension of the GVO method and how this can be used to derive time series of the magnetic field gradient tensor elements.

Geomagnetic Virtual Observatory Method
The Geomagnetic Virtual Observatory method allows for epoch estimates of the magnetic vector field components at a given target point (referred to as a GVO target location) to be derived using available satellite measurements enclosed within a cylinder of radius 700 km during the course of four months.
A radius of 700 km enables enough data for computing reliable and independent GVO estimates every four months (Hammer 2018). From these measurements, provided in an ECEF coordinate frame given by the spherical polar components, B obs = (B r , B θ , B φ ), magnetic field residuals are calculated as (Hammer et al. 2021a) where model fields subtracted are: the main field (MF), B M F , for SH degrees n ∈ [1, 13] determined using the CHAOS-7.2 model [http://www.spacecenter.dk/files/magnetic-models/ CHAOS-7/], see also Finlay et al. (2020); the static lithospheric field, B lit , for SH degrees n ∈ [14, 185] determined using the LCS-1 model (Olsen et al. 2017); the large-scale magnetospheric and associated Earth induced fields, B mag , as given by the CHAOS-7.2 model parameterized in time by the RC index ; and the ionospheric and associated Earth induced fields, B iono , as determined using the CIY4 model parameterized by 90-day averages of solar flux F10.7 (Sabaka et al. 2018).
Note here that we remove values of the time-dependent main field and then at a later stage, add back main field values at the GVO target position and target epoch Hammer et al. 2021a). However, the precise choice of main field used in both steps is not crucial (Hammer 2018;Hammer et al. 2021c). Note that in removing estimates of the time-dependent main field in eq.(1) estimates of the SV field (up to degree 13) are also effectively removed for each data point during the 4 month time windows, which aids a robust estimation. Adding back a main field estimate at the GVO epoch time then synchronizes the data from the considered four month window to this common epoch.
This way of correcting the individual data includes spatial gradients of the main field model used, so the resulting local potential does no longer contain the tensor information from the main field model which has been removed. It is thus necessary to reinstate the tensor information from the main field model when computing the magnetic field gradient tensor in section 3.2. Note that despite implementing a dark quiet time data selection scheme and also having removed estimates of magnetospheric and ionospheric fields together with their associated Earth-induced fields, contamination from non-core electrical currents may persist in the residual GVO field eq.(1). Such contamination could be related to ionospheric currents such as the polar electrojets and F-region currents.
Next, for each data point, both the position of the data point and the corresponding residual magnetic vector, eq.(1), are transformed from the spherical system to a right-handed local topocentric Cartesian system (x, y, z) having its origin at the GVO target location, as detailed in (Hammer 2018, p. 64). At this specific GVO location (and only at this location), x points towards geographic south, y points towards east and z points radially upwards (Hammer et al. 2021a).
At the GVO target point, the unit vectors of the local Cartesian frame coincide with the spherical polar unit vectors, i.e. (ê z ,ê x ,ê y ) = (ê r ,ê θ ,ê φ ). Assuming that the magnetic field measurements are made in a source free region, the residual field, δB, is a Laplacian potential field that is approximately quasi-stationary (Sabaka et al. 2010). This means that a magnetic scalar potential, V , is associated with the residual field, which in the local Cartesian coordinate system can be expanded as a sum of polynomials having the form C abc x a y b z c (Backus et al. 1996). In this application we expand to cubic terms following Hammer et al. (2021a) V (x, y, z) = C 100 x + C 010 y + C 001 z + C 200 x 2 + C 020 y 2 − (C 200 + C 020 )z 2 + C 110 xy + C 101 xz + C 011 yz The means and differences of the residual magnetic field are linked to this potential via appropriate design matrices constructed as described in Hammer et al. (2021a). The coefficients of the potential are determined from a robust least-squares solution, which includes a) an a prior data covariance matrix derived from standard deviations of the residuals between the data (means and differences) and estimates of an un-weighted least-squares solution, b) a diagonal weight matrix consisting of robust (Huber) weights, using a scale constant of 1.5 (e.g., Constable 1988), and c) an additional down-weighting factor of 1/2 when data comes from Swarm satellites Alpha and Charlie, taking into account that these fly side-by-side and thus provide similar measurements. From these potential coefficients, a mean residual magnetic field for the given GVO target point position and epoch is computed  We recall here, that the CHAOS-7.2 main field up to degree 13 has to be reinstated, as time-dependent main field estimates were removed by eq.(1).
Following the same procedure as in Hammer et al. (2021a), we compute a global grid of 300 GVO's located on an approximately equal area grid based on the sphere partitioning algorithm of Leopardi (2006). The global grid of 300 GVOs are listed starting from a position of latitude 89.9 • N at longitude 0 • going to a position of latitude 89.9 • S at longitude 0 • . The distance between the GVOs in this grid is ≈ 1400 km, and with a target cylinder radius of 700 km close to 80% of the measurements are used. This choice of radius avoids overlap between the GVO data cylinders, such that each satellite measurement is only used once in the global GVO grid computation, which helps to minimize possible correlated errors and biases due to the GVO binning method (Beggan et al. 2009; Shore 2013). The GVO height above ground is taken to be 370 km and 490 km (approximate mean orbital altitude) during CHAMP and Swarm times, respectively. Figure 1 presents the available number of GVOs for each epoch (the maximum possible number at each epoch is 300). Especially when considering the early CHAMP period, there are several epochs with very few available GVOs. Such strongly depleted epochs are primarily caused by our restrictions to geomagnetically quiet and nighttime data, and to selection based on CHAMP data quality flags. Table 1 presents the Huber-weighted mean and root-mean-square (rms) of residuals between the satellite measurements (after applying the corrections of eq.(1)) and the field estimates as computed from the determined potential coefficients

Computing the Magnetic Field Gradient Tensor within the GVO framework
We now proceed to formulate the magnetic field gradient tensor and describe how it transforms from a spherical polar coordinate system to the local topocentric Cartesian right-handed coordinate system used in the GVO method. This transformation will allow us to compute GVO time series for the magnetic field gradient tensor elements in analogy to the concept of GVO vector field time series.
We begin by expressing the magnetic field gradient tensor in the local Cartesian system of the GVO method described in Section 3.1. This is given by (see Appendix A for full details) ( This is a second-order tensor where the minus sign comes from defining the magnetic field as the negative gradient of the potential. The gradient tensor elements are denoted here by [∇B] jk , where the first subscript, j denotes the vector component under consideration and the second subscript, k, denotes direction of the field derivative. Using the local cubic potential eq.(2) estimated from the residual magnetic field as described in Section 3.1, a second-order residual field gradient tensor at the GVO target point can be derived using eq.(3) as Because the magnetic field is a solenoidal vector field, the divergence is zero, such that the trace of the gradient tensor vanishes, i.e. tr(∇δB GV O ) = 2(C 200 + C 020 ) − 2C 200 − 2C 020 = 0, reducing the number of independent elements from 9 to 8. In addition to this, because the field is a Laplacian potential field, the curl of the field vanishes and the number of independent tensor elements reduces to 5; in other words, the magnetic gradient tensor is diagonally symmetric and its trace is zero (Kotsiaros & Olsen 2012).
Following the GVO algorithm, gradient tensor estimates from a field model then have to be added back at the GVO target point for harmonic degrees n ≤ 13. The time-dependent CHAOS-7.2 main field up to degree 13 has here to be reinstated, as it was removed by eq.(1), in order to synchronize the data to the particular GVO epoch. Removal of the main field estimates results in the local potential in eq.
(2) also lacking the main field tensor information; this must therefore be added back, in order to obtain the full magnetic field gradient tensor. To do this, we will need to consider how the gradient tensor elements in spherical polar and Cartesian coordinate systems are related. The magnetic gradient tensor elements as expressed in the spherical coordinate system are given by (Olsen & Kotsiaros 2011;

Secular Variation Signals in Magnetic Field Gradient Tensor Elements 11
Kotsiaros & Olsen 2012), see also Appendix A Here the first column of the tensor contains the derivatives of the magnetic field components along the radial direction, the second column contains the derivatives along the co-latitudinal direction and the third column contains the derivatives along the longitudinal direction. The gradient element in the first column and row contains one term only, the field derivative term e.g. ∂ 2 /∂r 2 , while the other gradient tensor elements in addition to this also have an additional field term i.e. ∂/∂r, ∂/∂θ or ∂/∂φ.
The transformations relating the gradient tensor elements in the local Cartesian system to the tensor elements of the spherical system, only at the GVO target location, are in the end given by the following simple relations (see Appendix A for a full derivation).
[ Having determined the potential from the residual magnetic field eq.(1), we can compute a residual field gradient tensor by eq.(4) and add back main field gradient tensor estimates from the CHAOS-7.2 field model using eq.(5) for SH degrees n ≤ 13, using the above relations, in order to obtain the required GVO field gradient estimates ∇B GV O . Note that this procedure is analogous to the procedure applied in deriving vector field GVOs where the main vector field is added back. The above procedure is then repeated at each GVO location and for each epoch to compute all the desired GVO field gradient time series.
Field contamination from non-core sources dominates the instrumental errors of the satellites.
Thus, the primary limitation in obtaining core field GVOs are contributions due to unmodelled fields, such as the polar electrojet, and not measurement errors (e.g. Finlay et al. 2017 Bendat & Piersol 2010), where e jk,i is the residual of the ith data element, M is the number of data in a given series and µ jk is the residual mean for a given component. Hammer et al. (2021a) computed similar uncertainty estimates for the vector components using the vector field residuals with respect to the CHAOS-7.2 model. We note here, that using the CHAOS model to quantify the variability of the GVO series, provides only a crude indication of the data errors, which are assumed time independent. Moreover, correlation in the data errors are presently not accounted for, although this is expected especially at high latitudes. It is however challenging to empirically estimate non-diagonal data error covariance matrices with the short GVO time series presently available. Further work is needed on this important topic.
As with the ordinary GVO vector field time series, we estimate GVO gradient tensor time series on a global grid of 300 GVOs. We compute the SV as annual differences at each GVO for each tensor element. In order to quantify the scatter level in each series, we then fit cubic smoothing splines to the time series, with a knot spacing of 4 months and a smoothing parameter determined using a generalized cross-validation (GCV) approach (Green & Silverman 1993). Table 2 presents the mean rms differences between the GVO SV gradient tensor elements and GCV spline fits, separated into polar and non-polar regions. These rms numbers provide an indication of the scatter level in the GVO SV gradient data derived from the CHAMP and Swarm measurements first using CHAOS-7.2 ) as a main field model. Comparing the numbers between CHAMP and Swarm, we see that overall the values are lower for Swarm, i.e. Swarm gradient tensor element SV time series have a lower scatter than similar series for CHAMP. We note that the misfits are considerably lower at Secular Variation Signals in Magnetic Field Gradient Tensor Elements 13 Table 3. Mean of the rms differences (in pT/km yr −1 ) between GVO SV series and GCV cubic spline fits for six of the gradient tensor elements. Results are shown for GVO SV gradient data derived from Swarm and CHAMP data using COV-OBS.x2 model (Huder et al. 2020) as MF model in the GVO processing.

Component
rms [rr] rms [θθ] rms [φφ] rms [rθ] rms [rφ] rms [θφ] [pT/km yr We also tested how the choice of main field model (used for subtracting and adding back main field estimates) would impact the results. We produced test GVO tensor element series from both CHAMP and Swarm measurements using the main field values for SH degrees n ∈ [1, 13] of the COV-OBS.x2 model (Huder et al. 2020). Table 3 presents the mean rms differences using the COV-OBS.x2 model instead of CHAOS-7.2. This results in almost identical misfit levels to the GCV splines (i.e. scatter), between the GVO gradient series during CHAMP and Swarm times, regardless of whether CHAOS-7.2 or COV-OBS.x2 is chosen.

Field Gradient Element SV time series
We begin by investigating the temporal behaviour of the annual differences of each gradient tensor element at an example GVO location above Honolulu ground observatory in Hawaii, from which there are well known vector field records. To do this, we compute dedicated GVO gradient element series above the Honolulu observatory using the method described in Section 3.2. Here we are motivated by studies which have pointed out a change in secular acceleration of the radial component in the Pacific and Swarm (red) times at altitude 500km for a case study above Honolulu, Hawaii. Units are pT/km yr −1 .
occurring around 2017 (Sabaka et al. 2018;Finlay et al. 2020). In particular, we are interested to see if it is possible to identify this event in the GVO gradient tensor time series, and how this will display in the various tensor elements. Figure 2 present plots of the SV for each gradient tensor element above Honolulu, showing the GVOs derived from CHAMP (in blue) and Swarm (in red) measurements.
For comparison purposes we have mapped the two GVO series to a common altitude of 500 km by subtracting the SV gradient field differences between the GVO altitudes and 500 km altitude using the  Figure 3. By visual inspection, we find local regions with similar temporal changes as those observed at the Honolulu SV gradient series. In particular, a distinct "V " shaped behaviour is found in the eastern Pacific region in a band stretching from latitudes 20 • S to 20 • N and longitudes 180 • to 220 • with a possibly related opposite "Λ" shaped behaviour in the western Pacific region from latitudes 20 • S to 20 • N and longitudes 120 • E to 180 • E. These regional changes occur over a time window of 6 years reaching amplitudes of about 15 pT/km yr −1 . Note that a "V "-shaped SV gradient time series means a strong positive change in the SA, while a "Λ"-shaped time series means a strong negative change in the SA. Though more complex to interpret, the other SV gradient tensor elements (not shown) also exhibit distinctive behaviour in the Pacific region. These observed changes in the SV gradient elements indicate regional jerk-type event happening in the Pacific centred on 2017. In the ionosphere external fields tend to be organized according to the geometry of Earth's main field, and their signal in the GVO series may therefore be grouped accordingly to magnetic latitude in quasi-dipole coordinates (Laundal & Richmond 2017). Here magnetic latitude ±70 • (dark blue curve) may be used to approximate the border between North/South Polar and Auroral zones, while magnetic latitude ±50 • (light blue curve) divides the North/South Auroral and Low-to Mid-latitude zones (Hammer et al. 2021a). In all of the SV gradient element maps, higher scatter are found at GVOs located in the Polar and Auroral zones, which is consistent with noise (unmodeled fields) from ionospheric and magnetosphere-ionosphere coupling currents. In addition, some rapid variations can be seen, which decrease in amplitude when going towards equatorial latitudes. An example of this is found along longitude 150 • W, stretching from latitudes 30 • S to 60 • S, as highlighted in the side-panels of Figure 3 at three selected GVO locations (marked in the global maps by the green ellipse) An important question is whether the prominent change in SV observed centred on 2017 is robust and of internal origin. To address this, we produced a set of the GVO SV series above Honolulu during Swarm time, testing a range of geomagnetic selection criteria. We considered five cases: Case A using a dark quiet time data selection removing estimates of the magnetospheric and ionospheric fields as described in Section 2, this is our 'preferred' criteria for studying core field variations, Case B using a dark quiet time data selection removing estimates of the magnetospheric field but not removing estimates of the ionospheric field, Case C using a dark time data selection but with neither ionospheric nor magnetospheric corrections, Case D using a quiet time data selection from both day and night ("all" local times) and removal of the magnetospheric field, and Case E using data from "all" local times, without any quiet time data selection applied and without removal of for magnetospheric or ionospheric fields. Here "dark" means the sun is required to be at least 10 • below horizon, and "all" means that no such selection is applied, i.e. sunlit data are also included. SV gradient series for all five data selection cases are shown in Figure 4. The black dots corresponding to Case A are for our "preferred" criteria were also used to derive the maps in Figure 3. All selection criteria results in the same overall temporal "Λ/V " shape behaviour with an amplitude of ≈ 15 pT/km yr −1 . There is no increase in amplitude when including more disturbed data. For example, comparing Case C (purple star) with Case A/Case B (black/blue dots) should expose a signal from a magnetospheric source; however, the same "Λ/V " shape behaviour appears in all three cases. Comparing instead Case E (yellow star) with Case C and Case D should expose an ionospheric signal, which is expected to be larger during sunlit conditions; even though more scatter is seen in Case E, the same overall "Λ/V " shape is clearly visible and of similar amplitude. These results are consistent with an internal origin of the 2017 SV impulse event.

Example Spherical Harmonic Models Derived From Gradient Data
In this section we demonstrate that spherical harmonic (SH) field models with high temporal resolution (4 months) can be built from the global network of GVO gradient tensor time series. We then use these models to investigate global changes in SA during Swarm time, and in particular, we analyse the possible benefits of the GVO field gradient tensor series over more standard GVO vector field series. The magnetic field can be described as the negative gradient of internal, V int , and external, V ext , potentials (e.g. Sabaka et al. 2010)

Secular Variation Signals in Magnetic Field Gradient Tensor Elements
Here the internal and external potentials at a given epoch can be expressed by truncated spherical harmonic expansions up to degree N int and N ext , respectively [q m n (t)cos(mφ) + s m n (t)sin(mφ)] r r a n P m n (cos θ) , where r a = 6371.2km is the Earth's mean spherical radius, n and m are the SH degree and order, respectively, and P m n are the associated Schmidt semi-normalized Legendre functions. {g m n , h m n } and {q m n , s m n } are the internal and external expansion coefficients, respectively. Assuming an internal field source, the linear forward problem of determining the SH expansion coefficients can be written where d is a data vector containing the GVO epoch data (i.e. field vector components or field gradient tensor elements), G is the design matrix for an internal potential relating each model coefficient to either the GVO epoch vector or gradient data (for explicit expressions of the gradient tensor elements, we refer to Appendix A of Kotsiaros (2012)), and vector m contains the parameters of the potential, i.e. the internal SH coefficients here denoted as g m n and h m n for order m and degree n. For each GVO epoch we estimate a single epoch SH model using an Iteratively Reweighted Least Squares method which is applied in order to mitigate the impact of non-Gaussian distributed data residuals (e.g. Olsen et al. 2006) where we assign weights according to the diagonal weight matrix W = w/σ 2 consisting of Huber weights, w, having a Huber tuning constant of 1.5 (e.g., Constable 1988) and error estimates, σ 2 , of either the field vector components or field gradient tensor elements, depending on which type of data is solved for. These error estimates are derived as the total mean squared error (Bendat & Piersol 2010) of the residuals between the GVO vector/gradient tensor series and vector/gradient estimates of the CHAOS-7.2 model up to degree 16 (see Section 3.2). That is, the weight matrix consists of the inverse data covariance matrix elements multiplied by Huber weights applied in order to account for non-Gaussian data residuals. The weights are then iteratively updated until convergence. We derive models up to SH degree N int = 14 for each 4 month interval.
To investigate the 2017 region jerk event we next compute the secular acceleration change for each gradient tensor element between 2015.5 and 2018.5 at the Earth's surface Plotting global maps of this change in Figure 5 for each of the gradient elements for degrees n ≤ 9 at the Earth's surface, distinct patterns of SA change are seen to have occurred in the Pacific region Next, we seek to further inspect and compare the SH models obtained from GVO gradient series with similar models obtained using more traditional GVO vector component series. The left plot of Figure 6 presents the mean of the MF (dashed-dotted), SV (solid) and SA (dashed) spatial Lowes-Mauersberger power spectra , at the Core-mantle boundary (CMB) obtained by computing the time-averaged spectrum taken over all the epoch spectral lines. Spectral curves in blue and red are derived from GVO gradients and GVO vector series, respectively. The SV and SA power spectra derived from GVO gradient tensor series are seen to diverge less rapidly as compared to models derived from GVO vector series, and the SV/SA intersection happens at a slightly higher degree (11 compared to 9). This behaviour is consistent with the analyzes of Kotsiaros & Olsen (2014), who found that gradient data better constrain SV to higher SH degrees than vector data.
In order to characterize the influence of satellite configuration on the GVO gradient series, we consider case studies consisting of: 1) a single satellite, 2) two satellites at different altitudes which would provide a better local-time coverage as compared to the first case, 3) East-West capability from at the equator) . The right plot of Figure 6 presents the time-averaged MF (dasheddotted), SV (solid) and SA (dashed) spatial power spectra, derived from each of these three case. Note here, that all of the MF spectra overlap and are thus hidden by the red SAT-ABC spectrum. In addition, note that this red SAT-ABC spectrum in the right panel refers to the same model as the red MF-GVO Gradient spectrum in the left panel of Figure 6. We discuss these results and their implications in more detail in the next section.
Investigating further these SH models, Figure 7 shows the first time derivative of the internal expansion coefficients, computed based on simple first differences, derived from internal models using the GVO vector (blue, and also including an external expansion, in green) and GVO gradient (red) series derived using all three Swarm satellites. Example coefficients are shown for zonal, m = 0, terms (top row), tesseral, m = n, terms (middle row) and sectorial, m = n, terms (bottom row). To quantify the scatter level in the epoch coefficient series, standard deviations between the coefficient series and GCV smoothing spline fits (solid curves) are given in each case. Although robust estimation has been applied when deriving these models, outliers can be seen in the both series. A change in the sign of the trend in the SV signal is evident, especially in the sectorial coefficientsḣ 3 3 ,ḣ 4 4 ,ḣ 6 6 , but also in theḣ 3 5 coefficient centered on 2017. The scatter level in the near zonal coefficients of the GVO vector series are reduced if an external SH expansion (in green dots and curve) is included as well.
This points to the presence of external contamination in the vector series, which seems to be reduced in the gradient series. Figure 8 collects such standard deviations for SH degrees and order up to 12, from models derived using the GVO vector (left plot) and GVO gradient (right plot) series. We generally find less scatter in the GVO gradient series, and especially in the zonal and near-zonal (where m is close to zero) coefficients. For the sectorial terms the scatter levels are low for both series. Use of the vector series results in higher scatter levels for the near-zonal terms for degrees n > 2 and orders m ≤ 2. When also including an external SH expansion, i.e. using both internal and external expansion terms up to degree 14, for the GVO vector series (middle plot), we are able to reduce the scatter level in the near-zonal coefficients (middle plot of Figure 8), illustrating that using the GVO gradient elements in SH modelling helps in excluding external field signals. Note here that we focus on comparing SH models derived from GVO vector data with those derived solely using GVO gradient data (including an external SH expansion for the GVO gradient series would require vector information as well to obtaining a robust estimation). Global maps (not shown) show that much of the enhanced scatter is related to signals in polar regions being spuriously mapped into the internal field.
Finally, we show in Figure 9, the standard deviations from SH models based on the case study GVO gradient series, which were derived using Swarm Alpha only (case 1, left plot), Swarm Alpha and Bravo (case 2, center plot), and Swarm Alpha and Charlie (case 3, right plot).

DISCUSSION AND CONCLUSIONS
In this study we have extended the existing GVO concept and derived time series of the secondorder gradient tensor elements of the geomagnetic field at a global network of 300 locations. We have computed such GVO gradient time series from the mean and differences of vector magnetic field measurements, along track and in the east-west direction, from the low Earth orbiting CHAMP and Swarm satellites.
Inspecting SV gradient tensor elements for a GVO located above the Honolulu ground observatory we found evidence in the gradient series for a regional jerk-type event centered on 2017, observed as a In addition to the signature of geomagnetic jerks, global maps of the GVO SV tensor element series, also show some example rapid smaller amplitude SV fluctuations within a few years, especially noticeable in the d [∇B] rr /dt element (see left panel of Figure 3). Some of these small amplitude variations in the SV gradient series might be indicative of external field leakage. The most obvious sources of external field leakage are 1) rapid variations in polar latitude GVO series that persist to lower latitudes along the same meridional, which could be an indication of contamination by polar electrojet or field aligned currents, 2) rapid variations at mid or low latitude GVO series seen at all longitudes that could be caused by a signal having magnetospheric origin. We find no signs of distinct temporal variations along all longitudes at mid/low latitudes, suggesting that remaining magnetospheric disturbances are likely small. We do find some examples of rapid variations that seem to diminish towards lower latitudes, which could signify contamination from polar electrojet or field aligned currents systems.
Longer time series are needed in order to study the origin of such variability in the GVO gradient series.
In order to test for possible improvements in retrieving the SV signal using GVO gradient tensor data as compared to GVO vector data, we produced simple unregularized SH field models built from the GVO gradient and vector data derived using Swarm measurements. Comparing the power spectra of these models shown in Figure 6 supports the findings of Kotsiaros & Olsen (2014), that harmonics of the SV above degree 6 can be better resolved when using gradient tensor data than using vector data. In particular, analysis of the first time derivatives of the SH coefficients as presented in Figure   8, shows that especially zonal and near-zonal harmonics of models derived from GVO gradients have less scatter compared to similar models derived from GVO vector data. We find that the scatter level in the SV gradient tensor elements is higher during CHAMP times compared with Swarm times. The orbital configuration of the three Swarm satellites is advantageous for computing the gradient tensor as more data are available and Swarm Bravo, having a slightly higher altitude than Alpha and Charlie, provides information on the radial gradient which enables better potential determination. The right panel of Figure 6 shows that having two satellites, one of which is at a higher orbit, slightly lowers the spatial power (green) as compared with having one satellite (dark blue). In addition, from Figure 9, having two satellites lowers the scatter of the tesseral coefficients for m > 2. Considering instead two satellites in a side-by-side orbit, the associated spectral line (light blue in the right panel of Figure 6) is only marginally lower, but differences are most obvious from Figure   9, where the scatter level of the tesseral coefficients is smaller for orders m > 2. The zonal terms for n < 7 are found to have higher scatter, but this is expected as East-West information alone does not allow determination of zonal terms (Kotsiaros & Olsen 2012). Comparing the three case studies with the original case of using all three Swarm satellites, the associated spatial power (red in the right panel of Figure 6) is lower, while the scatter level across all SH coefficients as shown in Figure 8 is the lowest, except for lowest zonal terms which seems to be related to including the East-West gradient information as shown in the right plot of Figure 9. In order to mitigate such behavior of the low degree zonal terms, in future studies it might be worth exploring a Selective Infinite-Variance Weighting (SIVW) approach (Kotsiaros & Olsen 2014;Sabaka et al. 2013), wherein weights are applied to data subsets that are more sensitive to certain parameter subsets. By this approach, the East-West data could in the future be given less weight with regard to parts of the GVO potential that provide information of the zonal terms. The larger misfits during CHAMP times, may also result from a closer proximity of the lower flying CHAMP satellite to ionospheric current systems. From our tests we therefore conclude that the Swarm satellite trio is advantageous for deriving GVO gradient series. Having satellites at different altitudes better fills the 3D space of the GVO data cylinder, improving the recovery of the gradient quantities, for example improving the determination of the tesseral harmonics. Using all three satellites, enhances the recovery of all harmonic coefficients, and is clearly superior to having a single satellite.

AVAILABILITY OF DATASETS AND MATERIAL
The GVO gradient tensor data underlying this article and their associated uncertainty estimates are available from https://data.dtu.dk/, at (Hammer et al. 2021b where the summation convention has been used. The Cartesian coordinates are related to the spherical coordinates via (Riley et al. 2004, p. 363) The magnetic scalar potential, V , can be considered as a tensor of zero-order (or rank). The gradient operator in the generalized coordinates u p , where p = 1, 2, 3, having the covariant basis e p = ∂r/∂u p and contravariant basis e p can be defined as (Riley et al. 2004;Casotto & Fantino 2009) Applying the gradient operator to the potential generates a new tensor of one order higher, which is the first-order tensor (vector) describing the magnetic field where we use the comma notation V ,p = ∂V /∂u p to denote partial differentiation. The components that arise by applying an operator such as eq.(A.3) which is represented in the covariant basis, are here referred to as the "natural" components. In contrast, describing the gradient in terms of the normalized basis vectors, the "physical" components of the first order gradient can be expressed as That is, we use the "*" notation to denote the physical components represented in the normalized basis, and distinguish it from the natural components. The metric scale factor h p is determined by the elements of the metric tensor g pq = e p · e q (which completely characterize any curvilinear coordinate system) as h p = √ g pp . Note that the metric tensor also facilitates the conversion between covariant and contravariant bases (Riley et al. 2004;Casotto & Fantino 2009). Applying the gradient operator to eq.(A.3), produces the second-order gradient operator where Γ s pq denotes the Christoffel's symbols of the second kind (an array of numbers describing the derivatives of the covariant basis vector along that same basis), which can be expressed in terms of the metric tensor as (Riley et al. 2004, p. 814) Applying the operator eq.(A.6) to the magnetic potential V generates the second-order magnetic gradient tensor elements which can be written using the Christoffel's symbols (A.8) Here the semicolon notation is used to denote covariant differentiation. In addition, indices which occurs after a semicolon, means differentiation. In physical components eq.(A.8) is written as Recall here, that the "*" denotes the physical components. An essential aspect of the first-and second-order tensors is how their physical elements V * ,p or V * ;p q in one coordinate system u p transforms to a new coordinate system u p (Casotto & Fantino 2009;Riley et al. 2004, p. 811) where the partial derivatives ∂u p /∂u p are expressed by the Jacobian matrix. The Jacobian matrix times the metric scale factor term, can be regarded as a rotation matrix such that we may re-write eqs.(A.10) and (A.11) having the transformation matrix determined as

Secular Variation Signals in Magnetic Field Gradient Tensor Elements 33
scale factor ratios between the two coordinate systems using the physical components. Thus equations (A.10) and (A.11) (equivalently eqs.(A.12) and (A.13)) allow us to transform the tensors in one coordinate system, for instance the global (x,ỹ,z), to another, for instance (r, θ, φ). Note here, that for the equivalence of eqs.(A.12) and (A.13) in the natural components, the transformation matrix is defined by eq.(A.14) but without the matrix D, i.e. without the metric correction. Let us now consider the two transformations: a) Transformation from the global Cartesian (x,ỹ,z) to the spherical system (r, θ, φ) b) Transformation from the spherical system (r, θ, φ) to the local system (z, x, y) First, we specify the inner products of the basis vectors, the covariant metric tensors for the Cartesian system , which is valid for both the global and local Cartesian coordinate systems having basis vectors

Secular Variation Signals in Magnetic Field Gradient Tensor Elements 35
The first-order tensor (i.e. the magnetic field vector) in the spherical polar coordinates is given by eq.(A.5) using the metric scale factors eq.(A.18) The physical components of first-order magnetic tensor elements in the local Cartesian system, can be derived from eqs.(A.12) and (A.14), using the Jacobian matrix eq.(A.22). Here the transformation matrix R becomes the identity matrix because the Jacobian matrix, (∂u p /∂u p = ∂(r, θ, φ)/∂(z, x, y)), is the inverse of the D matrix. Therefore the physical elements in the local Cartesian system are identical to those in the spherical polar coordinates given by the relations In addition, from eqs.(A.5) and (A.18), and using the fact that eq.(A.17) also holds for the local Cartesian system having (h z = h x = h y = 1), we obtain the relations V ,z = V ,r , V ,x = 1 r V ,θ , V ,y = 1 rsin θ V ,φ . (A.25) The second-order tensor in the spherical polar coordinates is given by eq.(A.8) using the Christoffel's symbols from eq.(A.19) and metric scale factors eq.(A.18) ∇∇V = V ,rrêrêr + 1 r V ,θr − 1 r 2 V ,θ ê rêθ + 1 rsinθ V ,φr − 1 r 2 sinθ V ,φ ê rêφ + 1 r V ,rθ − 1 r 2 V ,θ ê θêr + 1 r 2 V ,θθ + 1 r V ,r ê θêθ + 1 r 2 sinθ V ,φθ − cosθ r 2 sin 2 θ V ,φ ê θêφ + 1 rsinθ V ,rφ − 1 r 2 sinθ V ,φ ê φêr + 1 r 2 sinθ V ,θφ − cosθ r 2 sin 2 θ V ,φ ê φêθ + · · · + 1 r 2 sin 2 θ V ,φφ + 1 r V ,r + cosθ r 2 sinθ V ,θ ê φêφ . , we can infer the relations between the gradient tensor elements in the local Cartesian system and the gradient tensor described in the spherical system. We can do this by recognizing, that the components in front of the elementary unit tensor (e.g. e rêr ) of eq.(A.26), are the physical components of the ∇∇V tensor, so that: V ,xy = 1 r 2 sinθ V ,θφ − cosθ r 2 sin 2 θ V ,φ V ,yy = 1 r 2 sin 2 θ V ,φφ + 1 r V ,r + cosθ r 2 sinθ V ,θ . (A.33) In order to express the gradient tensor in Cartesian coordinates, we note that the metric tensor becomes the identity matrix meaning that the metric scale factors h p becomes unity, and all of the Christoffel's symbols becomes zero such that the gradient tensor is given by eq.(A.8) where the minus sign comes from defining the field as the negative gradient of the potential.

Secular Variation Signals in Magnetic Field Gradient Tensor Elements 37
We recall that the semicolon notation denotes tensor elements following Casotto & Fantino (2009), and not the second order spatial derivatives. However, in the case of the gradient tensor in Cartesian coordinates, these two are equivalent cf. eq.(A.34).