Concurrent Empirical Magnetic Reconstruction of Storm and Substorm Spatial Scales Using Data Mining and Virtual Spacecraft

Data mining (DM) has ushered in a new era of empirical magnetic reconstructions of the magnetosphere via application of the k-nearest neighbors (kNN) method. In this approach, the combined magnetosphere storm-substorm state is characterized by the Sym-H and AL indices, their time derivatives, and the solar wind electric field vBzIMF. However, using the DM reconstructions to account for the substorm contributions to the ring current as well as describing storm-time substorms remains a problem. The inner region r ≤ 12RE, where the ring current develops, has a much higher density of data than the tail region 12RE ≤ r ≤ 22RE, where substorms operate. This results in two models inconsistent in their scales dictated by the corresponding data densities. The inner model reconstructs storm time dynamics, including the formation of the westward and eastward ring current and pressure distributions. The outer model captures substorm features, including the thinning and rapid dipolarization of the tail sheet during the growth and expansion phases, respectively. However, the substorm model is insufficient to reconstruct the eastward ring current while the storm model cannot fully reproduce substorm effects because it overfits in the tail region. This issue is addressed by constructing a hybrid model which is fit using virtual magnetic field observations generated by sampling the other two models. The resulting merged resolution model concurrently captures the spatial scales associated with both storms in the inner region and substorms in the near-tail region. Hence it is particularly useful for investigation of the storm-substorm relationship, including storm-time substorms and the impact of individual substorm injections to the buildup of the storm-time ring current.


INTRODUCTION
Storms and substorms represent two major modes of magnetospheric activity and the resulting space weather (e.g., [1,2]), which is reflected in the corresponding low-and high-latitude geomagnetic indices, e.g., Sym-H [3] and AL [4]. Substorms frequently occur during the development of storm main phases and to a lesser extent during the recovery phase [5]. Although the original idea of substorms viewed them as building blocks of storms [6], which is reflected in the names of these phenomena, that view is now strongly modified [7,8], as it is clear that storms and substorms are related because of their common driver, the solar wind [9,10]. Despite this, storm-time in-situ observations during substorm expansions show charged particles injected deep into the inner magnetosphere [11,12]. Furthermore, statistical analyses of particle measurements from the RBPSICE and HOPE instruments on the Van Allen Probes identified an energization of the ring current during substorm expansions [13,14]. Additionally, empirical magnetic field models found that the substorm expansion phase was correlated with an increase in the amount of current contained within the symmetric ring current (SRC) during both storm and non-storm substorms [15]. As both types of disturbances have dramatic impacts on the global configuration of the geomagnetic field and associated current systems, a concurrent description of storms and substorms and their phases is needed.
Since the development of the earliest empirical magnetic field models [16], the storm state of the magnetosphere was considered in their construction by binning the magnetometer data by storm activity level. Over time, empirical models became more complex in their description of the storm state, by incorporating additional current systems and making them functionally dependent on storm activity indices and solar wind conditions [17][18][19][20]. Data mining (DM) ushered in an entirely different approach to the empirical modeling of storms, the first of which was termed the TS07D model [21,22]. The DM technique refits the model for each snapshot in time to a small subset of the entire magnetometer database. This subset of data is identified by mining the whole database for other time intervals when the magnetosphere was presumably in a similar storm state configuration, characterized using the storm index Sym-H, its derivative, and the solar wind electric field parameter vB IMF z . The DM algorithm employed in TS07D is the k-Nearest Neighbor (kNN) method [23], which is described in section 2.2. The first applications of the TS07D model were to investigate the global magnetic field and current system configuration of storms driven by different solar wind phenomena [22,24].
In contrast to the storm state, until recently, inclusion of the description of the substorm state of the magnetosphere within empirical magnetic field models has been limited. Eventoriented approaches sought to describe the substorm growth phase by incorporating a magnetotail thin current sheet (TCS) into storm models and then adjusting the TCS to match particular event observations [25,26]. Wire-models attempted to describe the substorm expansion by hand tailoring the substorm current wedge (SCW) [27][28][29]. Showcasing the flexibility of the DM approach, TS07D was customized to picture geomagnetic substorms [15,30] by including the substorm index AL and its time derivative as binning parameters and by incorporating a TCS into the model structure. Termed the SST19 model, it was successful in reconstructing the primary substorm features of the magnetosphere including the stretching of the magnetotail associated with the enhancement of the TCS in the near-tail during the growth phase and the rapid dipolarization of the magnetotail along with the formation of the substorm current wedge during the expansion phase. It also revealed a connection between the substorm expansion phase and the ring current enhancement, seen as ≈ 1 MA increase in the dayside ring current during a non-storm substorm.
However, the resolution of the SST19 model was insufficient to fully resolve the innermost magnetosphere, in particular, it was unable to reconstruct the eastward component of the ring current. Earthward (r 4R E ) of the dominant westward ring current, the azimuthal component of the ring current changes sign becoming eastwardly oriented [31]. Assuming force-balance and pressure isotropy, the boundary between the eastward and westward oriented ring current identifies the location of the plasma pressure peak [32].
An advantage of the DM approach, which reflects a postulate of machine learning techniques, is that the addition of more high quality data sets allows for increasing model complexity. Thus, an obvious remedy to enable the model to more fully resolve the inner magnetosphere is to increase the resolution of the equatorial field, taking advantage of the numerous spacecraft missions that fly through this region (e.g., Van Allen Probes, THEMIS, Polar, Cluster). Indeed, several applications of the TS07D family of models have demonstrated this [33][34][35][36], but unlike SST19, the focus of these studies was largely limited to the inner magnetosphere. As this study reveals, the higher equatorial resolution used for those inner magnetosphere investigations overfits the near-tail, introducing numerous artifacts in the reconstruction of the magnetic fields and electric currents. On one hand, this is a key advantage of the TS07D approach, that is, the resolution of the model can be customized to the particular region of interest, but on the other hand, it is also a shortcoming in that no single resolution is adequate for the entirety of the spatial domain of the magnetosphere. This dilemma is discussed in detail in section 2.4.
In section 3 we present a simple yet effective solution for how to concurrently reconstruct both the inner magnetosphere necessary for storm features and the near-tail needed for the description of substorms. Two separate models are constructed with varying equatorial resolutions, one customized for the inner and the other for the near-tail regions. Both models are then sampled in their respective regions to synthesize a distribution of virtual magnetic field observations dense enough to fit a third model, effectively merging the two other models while maintaining the divergenceless property of the magnetic and current density fields. The focus of this study will be a pair of nonstorm time substorms that occurred on 8 March 2009. These were previously analyzed by Stephens et al. [15] and as such represent a good test case for comparing the merged model to the SS19 model, as is described in section 3.2.

Magnetic Field Architecture
Empirical magnetic field models are designed with two primary considerations: the spatial structure of the current systems and their dynamical evolution over time. For describing the spatial structure, it is useful to model each current system individually as a sub-model termed a module [20]. The total magnetic field B tot is then the sum of magnetic field of each module along with the internal field, e.g., B tot = B FAC + B PRC + B SRC + B tail + B MP + B int corresponding to the magnetic field from the field-aligned currents (FACs), partial and symmetric ring current (PRC and SRC), the cross-tail current, and the magnetopause (or Chapman-Ferraro) currents, respectively. The internal field, B int , is readily determined with ground magnetometers and as such is not in the scope of this research and the IGRF model [37] is used to represent it. Earlier models sought to define these modules by hand, crafting a mathematical description based on the theoretical picture of the current system (e.g., [19]). Each of these modules will have non-linear parameters that determine their spatial scales and linear amplitude coefficients controlling their intensity. For example, a magnetic field module could take the form B ′ (r) = a 1 B(r; β 1 ), where a 1 is the amplitude coefficient and β 1 is a non-linear parameter defining the module's mathematical structure (e.g., the current system's spatial scale size or thickness). The total model's set of a i s and β j s are then fit to the available magnetometer data [20]. The dynamical evolution of the current systems can thus be introduced by simply making a i s and β j s functions of time. Some of the earliest models achieved this straightforwardly by binning the magnetometer data by the K p storm index and performing separate fits for each bin [16,38]. The proceeding models instead opted to make a i s and β j s as functions of solar wind conditions and geomagnetic indices [17,18]. Again, the mathematical structure of the functions was hand-tailored.
The TS07D [22] and derivative models [15,[33][34][35] utilized a wholly different approach that sought to eliminate many of the hand tailored elements, motivated by the principle that the data should dictate the model instead. First, all the equatorial field modules (SRC, PRC, and tail current) were replaced by a single regular expansion that had no predefined azimuthal or radial structure derived from the general magnetic vector potential solution of a thin current sheet in the cylindrical coordinate system [21] taking the form: where B (s) 0n , B (o) mn , and B (e) mn are the basis functions with azimuthally symmetric, odd (sine) symmetry, and even (cosine) symmetry, respectively; while a (s) mn , a (o) mn , and a (e) mn are the corresponding amplitude coefficients determined in the fitting procedure. M represents the number of azimuthal harmonics (odd/even pairs) and N determines the number of radial (Bessel functions) harmonics used in the expansion. The thickness of the current sheet comes about by substituting z with ζ = √ z 2 + D 2 in the magnetic vector potential solution, introducing D as the characteristic half-thickness parameter. The SST19 model expanded upon this approach by including two such systems, one for the thick current sheet and one for the TCS, giving: where D TCS is the half-thickness for the TCS. This TCS system is key for reconstructing the enhancement of a thin cross tail current sheet which acts to thin and stretch the magnetotail during the substorm growth phase [15]. Later storm investigations also found the TCS facilitated in reconstructing the eastward ring current during quiet and weak storm times [35]. The FAC module is similarly mathematically described using a regular expansion, in this case a Fourier series in local time [39] which is duplicated and initialized to different latitudes to mimic an expansion [34] in both local time and latitude. The SST19 configuration is used here, which employs four different latitudes with the first four Fourier harmonics, totaling 16 total basis functions which describe the FAC field structure [15]. This factor of four increase over the original TS07D configuration proved critical in reconstructing a realistic FAC morphology associated with substorms (e.g., [40]). The FAC current sheets are bent to flow along approximately dipolar field lines [41] and are allowed to expand and contract by introducing two global rescaling factors κ R1 and κ R2 . κ R1 applies to the 8 basis functions at higher latitudes (region-1 or R1 FACs) while the κ R2 corresponds to the 8 at lower latitudes (region-2 or R2 FACs).
Each current system along with B int is given a complementary shielding field together represented as B MP which acts to contain B tot within the magnetopause boundary [20]: B tot ·n| S = 0, where S is the modeled magnetopause boundary [42].

kNN Method Application
The second key element of the TS07D model is the DM approach. DM identifies time intervals when the state of the magnetosphere was in a similar global configuration as the moment of interest. These time intervals are intersected with the historical magnetometer database to form a subset of data that is then used to fit the model's a i s and β j s. The procedure is repeated for each step in time giving the a i s and β j s their time dependence, thus fulfilling the motivating principle that the data should dictate the model's structure and dynamic evolution. The DM algorithm employed is k-nearest neighbors (kNN), where the global configuration of the magnetosphere and its dynamics are assumed to be represented by some finite dimensional state-space constructed from global parameters, such as geomagnetic indices or solar wind conditions and their time derivatives [43].
As the magnetosphere evolves in time, it traces curves, represented by the state-vector G(t), within state-space. Here, as with previous studies, the state-vector is discretized to a 5 min cadence, forming a cloud of points. Similar dynamical events (such as storms and substorms) will trace similar curves in this state-space, meaning for a moment of interest G(t = t ′ ) there will be other points in proximity to it, termed nearest-neighbors (NNs). The number of NNs used (K NN ) is much larger than unity but much smaller than the total number of points in the statespace database (K DB ): 1 ≪ K NN ≪ K DB . The distance between an NN point G (i) and the moment of interest G(t = t ′ ) = G (q) is determined using the standard Euclidean distance metric: Frontiers in Physics | www.frontiersin.org where each component is standardized by dividing by its standard deviation σ G k . Here, the combined storm-substorm 5D state-space from SST19 is utilized: (8) where Sym-H * and AL are common indices used to measure storm and substorm intensities, respectively [3,4,44]. The correction Sym-H * = A · Sym-H − B · P dyn is performed to approximately isolate the ring current contributions to this index (by removing magnetopause and induction fields) [45].  [46] and substorm [47] activity. The integrals notated by . . . | in Equations (4), (6), and (8) act to smooth the inputs. The storm intensity parameter Sym-H * | and the substorm intensity parameter AL| use half-cosine smoothing windows sized to the characteristic time scales of storms and substorms with st = 12 h and sst = 2 h, respectively. vB IMF s is instead smoothed using the more responsive exponential smoothing window with τ 0 = 0.5 h and τ ∞ = 6τ 0 , which better captures the substorm growth phase. Also included in the state-space are the smoothed time derivatives of Sym-H * and AL constructed using derivative windows (5) and (7) represented by the notation D . . . |/Dt [48]. It is critical to include the derivatives in the state-space as they act to differentiate between storm/substorm phases (main/expansion vs. recovery phases). Each of the components are standardized by dividing them by their standard deviations computed over the entire state-space.
Not all the hand tailored elements were removed from the TS07D model, in particular, as with earlier models, some dynamical features are explicitly built into the model structure, including the contraction/expansion of the magnetosphere in response to the changes in the solar wind dynamic pressure P dyn and the warping of the current sheets due to dipole tilt angle effects [21]. The contraction/expansion can readily be modeled by assuming all the current systems change in a self-similar way, that is, by using a simple spatial rescaling: r ′ ∝ P Dyn −κ r (e.g., [20]). Meanwhile, near the Earth (r 4R E ), the geometry of the current systems tends to be oriented with respect to the geodipole axis, while further down the tail (r 8R E ), the current geometry is controlled by the solar wind flow direction. These dipole tilt angle effects are accounted for in the model by application of the general deformation technique [49,50]. For example the flat current sheet described by Equation (1) is warped to account for the "bowl-shaped" deformation when the dipole tilt angle is non-zero [51], introducing three additional non-linear parameters (the hinge distance R H , the warping parameter G, and the twisting parameter TW), yielding B (eq) (ρ, φ, z) [35].
It is important to mention here that the kNN DM technique is an instance-based machine learning method [52] which drastically differs both from model-based ML methods, such as ANNs [53] and classical Tsyganenko models [17,19]. In the case of Tsyganenko models, tuning the parameters other than the linear regression coefficients and selected non-linear parameters is redundant, because the resulting model is universal and its architecture is custom-made and fixed from the outset. In case of ANNs, it is usually possible to split the model into training and validation sets and to use the latter to further optimize the model architecture [parameters, like K NN , or M and N in Equation (1) that are termed hyper-parameters]. This can also be done when kNN is used for data classification or when the binning and fitting spaces are the same. However, this is almost impossible in our case, when the binning and fitting procedures are made in different spaces, the global index state-space (4)- (8) and the real space, which is extremely sparsely filled with data, with no more than a dozen of probes available for validation at any given moment.
At the same time, the selection of the kNN binning space (4)-(8) can be optimized in our DM method using important physics constraints and dynamic systems theory. First, we take explicitly into account the storm and substorm states of the magnetosphere that are known to be assessed by indices Sym-H [3] and AL [4], whereas their trends are described by the corresponding time derivatives (5) and (7). The latter can be extended to higher time derivatives following the idea of the time delay embedding in the non-linear time series analysis [48]. The averaging time scales in (4)-(8) are also physics-based and they are consistent with the observed characteristic times for storms [19,46,54] and substorms [55]. This physics-based optimization makes our kNN method similar to gray-box models (e.g., [56]).
Further selection of the hyper-parameters is done as follows. For a given set of the magnetic field model complexity (parameters M, N, and the number of the FAC modules NFAC) there is usually an optimal range of K NN values where the reconstruction is stable (no overfitting) and yet resolves important storm or substorm features, such as the westward and eastward ring currents, thinning and dipolarization of the magnetotail. Then validation tests are performed within that optimal range (e.g., SI in [30] to quantify the reconstruction fidelity). These tests are discussed below in section 3.4. At the same time, the uncertainty of the kNN binning is independently quantified by comparing the original binning parameters (4)-(8) with their means and standard deviations within the bin as is also detailed in section 3.4. The optimal choice of K NN minimizes the bias and standard deviation of the NN means but at the same time avoids overfitting with the model structure parameters (M, N, NFAC, and others), i.e., K NN still needs to be large enough to resolve the corresponding storm and substorm structures and dynamics.

Magnetometer Database
The spacecraft magnetometer database used in this study is also the same as [30] spanning the years 1995 through 2018. It covers three distinct regions with different densities of magnetometer measurements as shown by the radial histogram in Figure 1. The inner magnetosphere region r ≤ 12R E has ample spatial coverage with the THEMIS [57] (five probes) and Van Allen [58] (two probes) missions sampling the inner equatorial magnetosphere, including the vicinity of the eastward current system (2R E ≤ r ≤ 4R E ) which is crucial to resolving the peak in the plasma pressure. The geosynchronous orbiting (r ∼ 6.6R E ) GOES (08, 09, 10, and 12) spacecraft reside within the ring current region. The THEMIS orbits have changed over the years but one of the primary apogees for the three inner probes has been r = 12R E , providing good coverage throughout the Inner Region.
The Near-tail Region 12R E ≤ r ≤ 22R E has a noticeable drop in data coverage, as only two outer THEMIS probes were ever located here and only for about 2 years before then were moved into a lunar orbit becoming the ARTEMIS mission [59]. With an apogee of r ≈ 18R E , the Cluster mission (four probes) helps populate this region, however, as a polar orbiting spacecraft, they spent a limited amount of time in the equatorial region.
Beyond 22R E the data density drops off by nearly an order of magnitude as the only spacecraft in the database that spent a considerable time in this region was Geotail. The near-earth reconnection sites are expected to be located here [60]. This motivated including the 2016-2017 MMS data. With an apogee of r ∼ 26R E , it nearly doubled the amount of data between 22R E ≤ r ≤ 26R E . The only other spacecraft included in the database is IMP-8, however, it comprises a relatively small amount of the total dataset. For each time step, the model is fit to the identified subset of magnetometer data by minimizing the root mean square of the difference between the model B (mod) and the observed B j,obs magnetic field vectors: where S NN is the number of data points in the magnetometer subset identified through the kNN technique. The model is evaluated at the spacecraft's position r. Two weight factors are incorporated into the objective function, one based on the data point's position in physical space w (0) (r) and the other on its position in state-space w j . The first, w (0) (r), lowers the weight factor in regions of the magnetosphere with a high density of data and was introduced to limit the bias of the fit toward these regions, in particular, to decrease the influence of the GOES satellites which are all located at the same radial distance [21]. The other weighting factor w j gives higher weights toward observations which correspond to NNs that are closer to the moment of interest in state-space and will be described in the next section. The set of non-linear parameters (D, D TCS , κ R1 , κ R2 , R H , G, and TW) are found by minimizing Equation (9) using the downhill simplex method [61], while the amplitude coefficients are solved using the singular value decomposition (SVD) pseudo-inversion method [62,63] also by minimizing (9).

Model Resolution Dilemma
The kNN approach has shortcomings caused by both the disparate density of NNs in the state-space and also the disparate density of magnetometer observations within the magnetosphere. The cause of the former is that, like many observed outputs of complex natural systems, geomagnetic indices, such as the substorm (AL) and storm (e.g., Sym-H * ) indices, tend to follow lognormal distributions [64][65][66], meaning the distribution contains many weaker storms and substorms than stronger ones. The result is that the NNs are inhomogeneous within the state-space, biasing the kNN method toward weaker events, which is especially problematic for strong and extreme events. Decreasing K NN helps resolve this problem but introduces another, overfitting. This problem was addressed in a pair of studies [35,36] that demonstrated introducing a simple distanceweighting of the NNs could significantly reduce this bias while keeping the K NN large enough to temper overfitting. The NNs are weighted using a Gaussian function of the form where R q is the distance of each NN to the query point and R NN is the radius of the NN n-sphere, which is the maximum R q (the distance between the query point and the most distant NN). σ determines the narrowness of the Gaussian and the value of σ = 0.3 is utilized here. These weights are then attached to the magnetometer datapoints when the model is fit by minimizing Equation (9). However, the second problem that was left unaddressed in those studies was the disparate data density in real/physical space caused by the distribution of spacecraft as shown in Figure 1.
The choice of K NN in the kNN approach is a tradeoff. A small number (K NN ∼ 1) is akin to event-oriented modeling but requires a similarly small number of degrees of freedom in the model, i.e., the combined number of scaling coefficients a i and non-linear parameters β j are also on the order of unity [25,26]. Increasing K NN permits more complex models with a larger number of degrees of freedom (more a i s and β j s), but when K NN approaches the size of the database the model becomes a universal statistical fit similar to classical Tsyganenko models (e.g., [17,19,38]) with a weaker sensitivity to storm and substorm phases. During the development of the SST19 model, it was found that an equatorial model resolution of (M, N) = (6, 8) using K NN = 32, 000 was sufficient for reconstructing the primary substorm configuration of the magnetotail while avoiding overfitting. For this study, the model described in section 2 using these SST19 values of (M, N) = (6, 8) and K NN = 32, 000 will be labeled the Tail Model. Figure 2 (right panels) displays the Tail Model's 2D equatorial distributions of the magnetic, electric current, and pressure fields during the late growth phase of a substorm (described below in section 3.2). For the sake of simplicity, these equatorial slices ignore dipole tilt and twisting effects, thus aligning the magnetic equator with the equatorial plane. Current densities are determined by numerically evaluating Ampere's law j = 1 µ 0 ∇ × B. The pressure is computed by integrating j × B radially inward starting at the boundary r = 10R E in the manner detailed in [35]: where the boundary pressure P(r 0 ) is assumed to be small and set to zero. Figure 2 demonstrates how the Tail Model performs well throughout the near-tail and midtail regions and thus is suitable for substorm studies. However, within the predominant THEMIS inner probes apogee (gray circle) the current densities appear underresolved. Specifically, the eastward currents are almost entirely absent ( Figure 2D). Without an eastward current, the integrand [−j × B] r monotonically increases earthward from the boundary r = 10R E meaning the pressure also increases monotonically. Thus, the location of the pressure peak is absent and the pressure radial gradients are also missing ( Figure 2F). In contrast, storm-time empirical reconstructions of the inner region (r ≤ 12R E ) were able to reconstruct the eastward current and thus the pressure, its peak, and gradient with the essential difference in the configuration of these models being a larger radial expansion number N = 20 used in the equatorial current module [33,35,36,67]. Furthermore, the inclusion of the TCS particularly helped in resolving the eastward current when the storm state of the magnetosphere was relatively quiet [35]. However, the addition of the TCS in the model resulted in reconstructed pressures with rather large azimuthal gradients. This was mitigated by halving their number from M = 6 to M = 3, yielding a smoother pressure distribution, resulting in an equatorial resolution of (M, N) = (3,20). Thus, the Inner Model is defined using this resolution of (M, N) = (3,20) with all other model configurations being the same as the Tail Model. The analogous reconstruction of the substorm growth phase using this Inner Model are plotted in the lefthand panels of Figure 2 demonstrating that it does indeed capture the eastward ring current ( Figure 2C) and appropriate pressure distributions ( Figure 2E). However, the magnetic and current density reconstructions now contain numerous artifacts throughout the near-tail region (Figures 2A,C) indicating overfitting there.
Herein lies the dilemma; no single resolution is capable of adequately reconstructing both the inner and near-tail regions using K NN = 32, 000 (the amount required to reconstruct the near-tail). Higher equatorial resolutions are necessary to describe the eastward current systems but will overfit the near-tail due to the lesser amount of data there and vice versa; lower equatorial resolutions perform well in the neartail but miss the eastward currents in the inner region. This results in a bifurcation of models: the Inner Model (higher resolution) is suitable for storm spatial-scale reconstructions and the Tail Model (lower resolution) is applicable to substorm scales. A potential solution is to dramatically increase the value of K NN , however, due to the disparate density of data in state-space, this would weaken the model's sensitivity to the event of interest and would begin to resemble statistical modeling instead of the DM approach sought. A simple solution would utilize a piecewise field, that is, to evaluate the Inner Model in the inner magnetosphere region and the Tail Model in the near-tail region. However, the equatorial current sheet described by Equations (1) and (2) ensures a divergenceless B and j fields. Such a piecewise field would introduce discontinuities which would violate these conditions and would also introduce infinitely thin current sheets. The resultant question is how to transition between these regions in a way that maintains ∇ · B = 0 and ∇ · j = 0? The next section 3 presents a simple resolution to this dilemma that smoothly transitions between the two regions while maintaining divergenceless fields.
It must be stressed that this resolution dilemma is not just a technical issue related to the sparse data distribution in space and the need to improve it is some regions. It reflects different physics processes associated with magnetic storms in the inner magnetosphere and substorms that reconfigure the magnetotail and create new FAC systems. In the earlier DM algorithms both storm and substorm descriptions used fleets of synthetic probes mined in the history of the magnetosphere. But these descriptions did not exchange the information gained from inner and outer magnetosphere description. So the fitting challenge resembles the first-principles model problem of the concurrent description of magnetic storms using the dedicated ring current models [68][69][70] and global MHD models [71][72][73] that often properly describe the outer magnetosphere. The resolution of that physics-based problem is similarly offered in the form of coupled models [74][75][76][77], when the information on the plasma conditions at some boundary separating inner and outer magnetosphere is transferred from MHD to ring current models. At the same time, the latter are used to adjust the equation of state in MHD taking storm effects into account.

Merged Model Algorithm
In this section, a simple method for concurrent reconstruction of both substorm and storm spatial scales is proposed, thus addressing the model resolution dilemma described above. The general approach is to introduce a third model, termed here the Merged Model, constructed using an equatorial resolution of (M, N) = (6, 20), the lowest common resolution between the Tail (M, N) = (6,8) and Inner (M, N) = (3, 20) models. However, instead of fitting the Merged Model to actual spacecraft magnetometer data points, it is instead fit to virtual data, which are simulated by randomly sampling the equatorial regions of the Tail and Inner models outside and inside the merging boundary, respectively. Because the Tail and Inner models can be sampled at any location, the introduction of these virtual datapoints allows the Merged Model to be fit with an arbitrary density of points. The result is a Merged model (given a proper distribution of virtual data points) that tends to reflect the Inner Model within the merging boundary and the Tail Model beyond it with a smooth transition between the two in the vicinity of the boundary. Also, because the Merged Model's equatorial field is still described by Equations (1) and (2), by construction ∇ · B = 0 and ∇ · j = 0 is ensured.
For simplicity, a cylindrical boundary is used identified by the merging radius R M . Initially, R M was defined to be commensurate with the demarcation between the inner and near-tail regions as characterized by the predominant THEMIS inner probes apogee (r = 12R E from Figure 1), though, this yielded unsatisfying results. Figure 3 shows 1D plots of the external magnetic field B ext z (Figures 3A,B) and the westward component of the current density j w = −j φ (Figures 3C,D) for the Inner (blue lines) and Tail (orange lines) models along the X axis. Two different moments are shown corresponding to the late growth (left panels) and expansion (right panels) phases, respectively. Ideally, there would exist a clear boundary where B ext z and j w from the two models intersect. Instead, the intersections (between the blue and orange lines) vary between ∼ 6R E -10R E . This justifies moving R M earthward of the THEMIS apogee of 12R E . Ultimately, the center of this region, that is, a cylindrical radius of R M = 8R E was chosen, which will be used as a constant merging boundary throughout the rest of the study. Although this is a rather arbitrary value, as shown below, this choice of R M enables the Merged Model to fully capture the deep B IMF z minimums and eastward currents of the Inner Model while also reconstructing the substorm scale features of the SST19 and Tail Model. The merging algorithm can be enhanced in future studies by incorporating a dynamical merging boundary that minimizes the differences of B ext z and j w between the Tail and Inner models. The use of elliptic cylindrical boundaries should also be pursued as Figure 3 indicates that day-night asymmetries exist in the optimal merging location.
To create the Merged Model, the non-linear parameters (D, D TCS , κ R1 , κ R2 , R H , G, and TW) are taken to be the average of their values from the Inner and Tail models and are not included in the fit; while the linear amplitude coefficients for the equatorial and FAC systems are fit in the manner described before, that is by using the SVD least-squares method to minimize Equation (9). However, instead of using real magnetometer data points, the B j,obs i terms in (9) are populated with virtual measurements constructed by evaluating the Inner and Tail models. The Inner Model was randomly sampled to achieve a data point density of 1R E −3 within the cylindrical volume 1.1R E ≤ r ≤ R M , −6R E ≤ z ≤ +6R E . The value of ±6R E was chosen as it is twice the typical value of the thick current sheet half thickness D ≈ 3R E . Meanwhile, given the lower radial resolution of the Tail Model, it was only sampled at a data point density of 0.2R E −3 within the cylindrical volume R M ≤ r ≤ 31R E , −6R E ≤ z ≤ +6R E . While this worked well for resolving the thick current sheet, due to the small height scale size of the TCS, additional samples were needed from the region −1R E ≤ z ≤ +1R E , where ±1R E is twice the typical TCS half-thickness D TCS ≈ 0.5R E now using six times the densities from before, 6R E −3 and 1.2R E −3 for the Inner Model and Tail Model, respectively. The reason for using a random distribution of points in contrast to a regular grid, is that fitting a regular grid of points with a magnetic field described by the regular expansion in Equation (1) (Figures 3A,B) and the eastward currents (Figures 3C,D), both of which are key characteristics needed to reconstruct the storm-time dynamics of the inner magnetosphere [33]. Tailward of R M the Merged Model closely matches the Tail Model (orange line). This indicates that the Merged Model should also reconstruct the primary substorm scale features, such as the thinning and stretching of the magnetotail during the growth phase and its rapid dipolarization during the expansion phase, as will be shown in the next section 3.2. In the ∼ ±2R E region bounding R M , the Merged Model smoothly transitions between the Inner and Tail models.
In order to qualitatively compare the Merged Model to the Inner and Tail models, 2D equatorial slices of the magnetic and current density fields are plotted for the late growth phase in Figure 4. Inside the merging radius R M = 8R E (yellow circles) the Merged Model qualitatively resembles the Inner Model, and importantly for the evaluation of the pressure, resolves the eastward currents (orangish-red region r 5R E in Figures 4D,F). Outside of the R M , the Merged Model (Figures 4C,F) is nearly indistinguishable from the Tail Model (Figures 4B,E). Notably, the Merged Model resolves the B z minimum formation in the tail at 10R E r 13R E and the enhancement of the cross tail current at 7R E r 16R E . The distribution of the virtual spacecraft projected into the equatorial plane used to fit the Merged Model is overplotted in Figure 4C, showing how random sampling of the other two models evenly fills in data gaps present in the near-tail and midtail regions, particularly beyond the Cluster apogee (r ≥ 19R E ). A note, as can be seen in Figure 4, the Inner Model is now only being fit using data within the predominant THEMIS inner probes apogee of r ≤ 12R E . Including data beyond this serves no purpose because the Inner Model is only being sampled within r ≤ R M when creating the Merged Model and their inclusion may bias its reconstructions.
The differences between the Inner/Tail and Merged models from the panels in Figure 4 are displayed in Figure 5. Within geosynchronous orbit, the differences between Inner and Merged models are relatively small (Figures 5A,C), with maximum values of B = 1.9 nT and j = 3.6 nA/sq.m and mean differences of B = 0.66 nT and j = 0.9 nA/sq.m. In this same region, the differences between the Tail and Merged models are about two to three times larger, with maximum values of B = 3.9 nT and j = 6.8 nA/sq.m and means of B = 1.8 nT and j = 2.5 nA/sq.m. The comparison between the Tail and Merged models shows negligible differences beyond ∼ 10R E (Figures 5B,D), while the equivalents for the InnerModel display large differences there. This confirms that the Merged Model largely mimics the Inner Model in the inner magnetosphere and the Tail Model in the near-tail. However, of interest are the differences in the vicinity of the merging boundary. Within the ±2R E region bounding R M , the mean differences are B ≈ 1 nT and j ≈ 2 nA/sq.m but show rather large maximum differences; B = 5 nT, j = 10 nA/sq.m for the Inner Model and B = 3 nT and j = 6.5 nA/sq.m for the Tail Model. The cause of this being the relatively large mismatch between the Inner and Tail models as is evident in the Figures 3, 4. This further supports additional investigation into a more optimal merging boundary in future works.

Merged Model Reconstruction of 8 March 2009 Substorms
Global reconstructions for the second 8 March 2009 substorm are displayed in Figures 6, 7, the first corresponding to the late growth phase at approximately the time of the substorm onset and the second is 25 min later during the expansion phase as indicated by the AL index (Figure 6a). These reconstructions highly resemble the reconstructions of the same event using the SST model [15] (Figures 3, 4 from that work). This confirms that the Merged Model is mostly analogous to the SST19 model throughout the near-tail region as was expected from the analysis shown in Figures 3, 4. However, the Merged Model now reconstructs the eastward currents as can be seen near the planet in Figures 6e, 7e.
During the growth phase, the magnetotail is stretched (Figure 6e), particularly in the region between ∼ 8R E and 15R E corresponding to the B z minimum region (Figure 6d). The TCS is especially intense on the night side in this same region (Figure 6e). The 66 and 68 • field lines cross the magnetic equator at 12R E and 18R E , respectively while the 68 • line is open. The height integrated thick current sheet is plotted (Figure 6f) as it is a better proxy for the content of the ring current compared to the equatorial slice of the current density (as is shown in Figures 2, 4) because the current density of the TCS tends to be the dominant source of current density along the magnetic equator. This panel shows a modest westward ring current within geosynchronous orbit r 6.6R E as well as a strong outflow to the magnetopause at the dusk terminator. The pressure during the growth phase is relatively azimuthally symmetric, with a pressure peak at 4-5R E .
During the expansion phase, which takes only 25 min, the global configuration of the magnetosphere is drastically altered (Figure 7). The magnetotail becomes much more dipolar. The 66 • field line (Figure 7e) now crosses the magnetic equator at 7.5R E (compared to 12R E in the growth phase) and the 68 • line crosses at 10R E (compared to 18R E ). The previously open 70 • field line now crosses at 14R E , indicating the conversion of open to closed flux presumably from reconnection. This dipolarization is congruous with the formation of the magnetic flux pileup in this region (Figure 7d). There is a strong enhancement of an westwardly directed thick current all across the night side taking the appearance of a PRC (Figure 7f). However, the SST19 analysis [15] along with magnetohydrodynamic (MHD) simulations [78] indicated that the enhancement of this PRC is associated with closure through the substorm current wedge. The substorm current wedge manifests as a eastwardly directed TCS in the vicinity of the magnetic equator, but its closure out of the plane is via a westwardly directed thick current, some of which closes through the ionosphere as R2 currents. Another point of agreement between the Merged Model and the SST19 model is the strengthening of the dayside thick current sheet within geosynchronous orbit (r ≤ 6.6R E ), which was interpreted as an intensification of the symmetric ring current.  Figure 4. The format is the same as Figure 4 but the range on the color bars has been halved to emphasize the differences.

Merged Picture vs. Earlier Substorm Reconstruction
The quantitative analysis performed using the SST19 model for the 8 March 2009 substorms [15] is now recreated using the Merged Model, in particular their Figure 6 is recreated and is shown as Figure 8. Compared to that substorm-focused reconstruction, the Merged Model reveals more profound variations of the R2 FAC amplitudes (Figure 8d), stronger buildups and decays of the TCS integrated current (orange line in Figure 8f) with significant negative values implying bifurcated current structures. One such bifurcated current structure is well seen in the meridional current distribution presented in Figure 7e.
The merged resolution picture reveals stronger variations of the B z magnetic field at the Van Allen Probes inner probes apogee (red line in Figure 8g). These stronger variations reflect the impact of distance-weighting the NNs along with a better resolution of the inner magnetosphere region resulting eventually in the resolution of the eastward current and the plasma pressure shown in Figures 6g, 7g. At the same time, it resolves the formation of the near-Earth X-line around the substorm onset (red and orange lines in Figures 8h,i), which is also seen in Figure 6d. Its resolution becomes possible not only due to the solved overfitting problem for the tail region, but also due to the use of the MMS data and the distance-weighted kNN algorithm as is elaborated in detail in [30].
The comparison also reveals some pitfalls of the present merging algorithm. In particular, unlike the SST19 model, the merged resolution model no longer has a clear correlation between the substorm phase and the amplitude of the R1 currents (Figure 8c). The probable explanation, is that because the R1 currents do not close through the equatorial plane (in contrast to the R2 currents), relatively few of the virtual spacecraft points sample them. To rectify this in future reconstructions using the Merged Model, the amplitude coefficients for the FAC systems should be set in a similar fashion to the non-linear parameters and not included in the fit. This is reasonable since the focus of    the merged modeling approach is to combine two models with different equatorial resolutions.
The primary enhancement of this Merged Model over the previous SST19 model is the resolution of the pressure in the inner magnetosphere. In order to quantify the change in pressure during the course of a substorm, the 2D pressure distributions from Figures 6g, 7g are averaged over all local times for a given radial distance r. The resultant average pressures P are plotted as a function of radial distance r for the second 8 March 2009 substorm in Figures 9D,E for several moments during the growth and expansion phases (the solar wind values and geomagnetic indices are included for context in Figures 9A,B). During the growth phase, the average pressure stays relatively stable, with a max value between ∼ 5-6 nPa located at 4.25-4.5R E . In contrast, during the expansion phase P rapidly increases reaching values as high as 10 nPa while the location of the peak stays relatively stable. The time series for P at 4.5R E is then plotted in Figure 9C for both substorms, showing that during quiet and growth phases the pressure is stable and low but then rapidly increases during the expansion phase. P then decreases during the recovery phase, although for the second stronger substorm, it does not return to the nominal value until more than 3 h after the substorm onset. Such an enhancement of inner magnetosphere pressure is consistent with in-situ observations [12] of particle injections as well as statistical analyses [13] of substorm expansion phases.

Uncertainty Quantification
An important aspect of empirical model development is its validation via comparison of the model output to the data that went into its construction. In this section, three different analyses are performed to test the model's fidelity. First, the DM approach is analyzed to ensure that the identified NNs characterize the 8 March 2009 substorm event and next the in-situ spacecraft magnetometer data is compared to the reconstructed field. Lastly, the three different models are statistically cross-validated using ∼28 days of model reconstructions to quantitatively assess the model error as a function of activity mode and radial distance.
The DM approach assumes that there are enough similar events in the state-space that the NNs collectively match the event of interest. This assumption can be tested by using the state-space analysis performed in earlier studies [30,35,67] in which the time evolution of the state-space components G i (t) (the binning parameters) for the studied event are compared to the mean of those parameters over the NNs G (NN) i (t). Ideally, G (NN) i (t) should match G i (t). Figure 10 shows the evolution of the five components of the combined storm and substorm state-space (black lines) during the two 8 March 2009 substorms corresponding to the quantities described by Equations (4)- (8). Panels a, c, and e demonstrate how the integral convolutions act to smooth the original data (pink lines). It should be stressed that the smoothing windows are asymmetric in time, that is, they are only computed over previous and not future data points [note the zeros on the upper bounds in the integrals (4)- (8)]. This was introduced to prevent the smoothed parameters from changing prior to the start of an event [48]. For example, without this asymmetry, the model would begin reconstructing storm activity prior to the arrival of the southward IMF in the solar wind. This causes the smoothed values to lag the original (black lines lag behind the pink lines). However, as demonstrated by Figure 8, this is not problematic as the time lag is universally applied over the entire dataset. Note how the modeled magnetic field (Figure 8g) and model parameters (Figures 8d-f) are still largely correlated with the original value of AL (Figure 8a solid orange line) and not its smoothed value (dotted orange line).
The weighted average of the state-space vector over the closest K NN = 32, 000 NNs is shown by the blue lines: G (WNN) = K NN j w j G j / K NN j w j , where the weight factor w j is computed using the Gaussian distance weighting from Equation (10). In a similar manner, the weighted standard deviations of these parameters over the set of NNs is also computed and the ±1σ values are indicated by light blue envelopes that surround G (WNN) (blue line) in Figure 10. Overall, G (WNN) largely stays within 1σ for the event of interest (black line).
The most notable inconsistencies appear in Sym-H * | and D Sym-H * |, which is a result of mixing parameters of different smoothing scales. Sym-H * and its derivative are smoothed over storm time frames of st = 12 h while the AL equivalents are instead smoothed over substorm scales sst = 2 h. This also demonstrates that the pressure-corrected storm indices contain contributions from substorm current systems and is not representative of just a pure ring current as has been discussed in previous studies [79,80]. Other deviations appear in particularly when AL| it reaches its minimum. The cause of this is the inhomogeneity of the datapoints in state-space as has been extensively discussed in earlier works [35,36,67]. This biases the DM approach toward weaker events as they occur more frequently which is only exacerbated as the events of interest become more extreme [36]. Distance-weighting the NNs can drastically mitigate this issue [35], although, as is shown in the second substorm, it does not entirely correct it. For instance, the minimum value for the second substorm is AL| = −576 nT while the mean over the NNs is only AL| (WNN) = −478 nT.
Next, the reconstructed magnetic field of the Merged Model is compared to the available in-situ magnetic field measurements from the THEMIS and GOES satellites (Figure 11). For this event, all five THEMIS probes are located in the magnetotail, with probes A, D, and E being in a similar orbital configuration situated at about post-midnight and r ≈ 12R E , while probes B and C are at ∼ 21 MLT with r ≈ 25 and 16R E , respectively. As the magnetotail current sheet thins, less magnetic flux threads through it resulting in a decrease in B z (e.g., [81]), as is observed by all five THEMIS probes during the growth phases, with B z approaching values of B z ∼ 1 nT prior to substorm onset (panels b, e, h, k, and n; black lines). The model captures these thinning signatures for the inner probes A, D, and E (panels b, k, and n; brown lines) as their location places them within the B z minimum and TCS region as was shown in Figures 6d,e. The model underestimates the thinning, as the B z values only reach B z ∼ 4 nT during the growth phases. One explanation is that the vB IMF s | parameter is a suboptimal proxy for the development of the TCS in the magnetotail while another is that the DM approach smears the singularly thin nature of the TCS effectively making it thicker than in reality.
Following substorm onset and throughout the expansion phase, magnetic flux is transported earthward accumulating in the near-tail region and enhancing the value of B z there (e.g., [82]), leading to a more dipolar magnetic field configuration. The three inner probes observe this dipolarization as a ∼ 20 nT increase in B z for the first substorm and a ∼ 30 nT increase for the second (panels b, k, n; black lines). As was shown in Figure 7d, the model reconstructs this flux pileup across the whole nightside magnetotail to a distance of r ≈ 16R E . The in-situ validation demonstrates that the model nicely captures these B z increases for the inner probes (panels b, k, n; black lines). While the model effectively captures the magnitude of the flux enhancement for the first substorm, it underestimates it for the second. This is not unexpected, as the state-space analysis shown in Figure 10 demonstrated the DM approach underestimated the intensity of the second substorm (Figure 10c). However, the model misses the TCS formation and the dipolarizations at the locations of outer probes B and C (Figures 11e,h) which are further down tail and whose MLTs are significantly away from midnight. Further, while the model displays good consistency with the GOES-10 and 12 magnetometers (Figures 11p,q,v,w), there are quite large deviations for GOES-11 (Figure 11t). This particular pair of substorms appears quite global in nature, showing signs of current sheet thinning as far as r ∼ 20R E down tail at 21 MLT and large dipolarizations at geosynchronous orbit at 3 MLT. While the AL index is a good indicator of the overall substorm intensity, it has its limitations. Firstly, it is derived from a rather small (10)(11)(12) number of magnetometer stations [44] and secondly, it yields no information about the local time configuration of the substorm currents [83]. As such, the DM approach will tend to reconstruct the typical substorm as is represented by the statespace parameters but will have no specific insight into the event's MLT configuration. This may be the underlying cause of the deviations between the model and the THEMIS B and C as well as the GOES-11 satellites. Some of these issues might be addressed by using another substorm index constructed from a more robust collection of magnetometer stations, such as the SML index [84], which furthermore, is computed for differing local time bins [85].
Overall, these results largely mirror to the validation of the SST model as shown in [15] (their Figure 2). This is not unexpected as the Merged Model closely matches the Tail Model at these spacecraft's locations and given that the Tail Model is very similar to the SST19 model. A shortcoming of the spacecraft configuration for this particular event is that none of them are earthward of geosynchronous orbit and according to Figures 3,  4, the improvement of the Merged Model over the SST19 model should be most apparent in this region. This will be addressed by a quantitative statistical analysis next.
In order to preform a quantitative uncertainty estimate the three models are cross-validated. In this technique, a subset of the data is reserved, that is it is not included when the model is fit, and instead forms a validation dataset. Statistical analysis is then performed on the independent validation set, and because this set was not used while fitting the model, it represents an out-of-sample test.
For the purposes here, the validation set should include spacecraft which adequately cover both the inner magnetosphere and near tail regions while there is also storm and substorm activity. The year 2015 was chosen to be the validation time interval as both the Van Allen probes and the THEMIS missions were sampling the inner magnetosphere, while the Cluster and MMS missions observed the near-tail. Further, the maximum sunspot count of solar cycle 24 peaked in April of 2014, so solar activity was still quite high throughout 2015 resulting in an elevated occurrence of geomagnetic storms [86]. Thus, when performing the model cross-validation, the entirety of the spacecraft dataset for the year 2015 was excised from the magnetometer database when fitting the models. That is, only data from spacecraft missions spanning the years 1995-2014 and 2016-2018 are used during fitting, while only data from 2015 are used in this cross-validation. Owing to the computational expense needed to fit the model, it was not feasible in this study to model the entirety of the year 2015. Instead, only times corresponding to storm and substorm activity from the 3 months of June, September, and December were modeled. At least one strong storm (Dst −100 nT) occurred during each of these months, and their separation in time allowed the spacecraft to sample a range of different local times. The start of September corresponds to the beginning of the MMS primary science mission phase, and as such is when the MMS magnetometer data becomes available. Here, storm and substorm activity time intervals are defined as when Sym-H * | ≤ −50 nT or AL| ≤ −300 nT. To ensure that broader portions of the storm main and recovery phases and substorm growth and recovery phases were also included, the time intervals were expanded by ±3 h and ±30 min, respectively. In total, the identified storm and substorm activity time intervals across these 3 months spans 673.5 h or about 28 days, corresponding to 8,046 model fits and ≈ 68, 000 spacecraft magnetic field observations when evaluated at a 5 min cadence. Here, the out-of-sample model error is quantified using the magnitude of the difference between the model evaluated at the spacecraft location and the observed magnetic field there: This quantity is similar to the model residuals which are key to its fitting, for example, as seen in the objective function (9), but without the weighting factors. They also have been the primary metric for testing empirical magnetic field model fidelity in previous studies (e.g., [87][88][89]). Specifically, the average magnitude of the difference between the model and observation ǫ = 1 N ǫ i is reported in Table 1, where N is the number of samples.
The Inner Model performs marginally better for the Van Allen Probes data than the Tail Model as their orbit is entirely within the inner magnetosphere region. As Figures 3, 4 indicated, this marginal improvement in B is crucial for reconstructing the current density j which depends on the spatial derivatives of B. In contrast, the Tail Model has lower errors for the Cluster spacecraft which spends relatively little time in the inner equatorial magnetosphere. The THEMIS and MMS missions have similar errors for both models. All models display notably lower errors for the less active times, presumable because the global configuration is more regular. Substorm and storm performance is rather similar, although, the storms during the validation interval contain significant substorm activity based on the AL index. Further analysis could be performed to separate isolated, multiple, and storm-time substorms, and storm intervals with and without substorm activity. Importantly, the Merged Model generally matches the lower error of the two other models, statistically validating the algorithm discussed in 3.1. Indeed, the Merged Model has the smallest error across the entire validation set at ǫ = 16.228 nT. To put this in context, running the same model cross-validation for the commonly used T89 and T96 yields errors of ǫ = 20.006 and 21.777 nT, respectively, while using no model (B ext = 0) gives is 37.458 nT.
To get a better indication of how the model errors change as a function of distance, the errors were collected into different radial bins ( r = 0.25R E when r < 12R E ; r = 1.0R E when r > 12R E , where r is the size of the bin) and for the different types of activity levels as are displayed in Figure 12. The Tail Model (orange lines) has significantly lower errors tailward of r = 12R E compared to the Inner Model (blue lines), although, recall the later was not fit using data in this region, so it should not be expected to perform well here.
Meanwhile, the Inner Model (blue lines) has lower errors within 3R E ≤ r ≤ 6R E , indicating that it does yield a more accurate reconstruction of the ring current region. However, errors are still relatively high here, particularly during storm times (Figure 12D), indicating the model still potentially underestimates the ring current intensity. One cause might be, given that there were several strong storms during this validation period (one storm reached a minimum Dst ≈ −200 nT while another hit Dst ≈ −150 nT), the aforementioned bias toward weaker events. Another, is that strong storms tend to have a complex morphology (e.g., [36]), and the multitude of mesoscale features can simply not be discerned with the DM approach. Also of note, the errors significantly increase earthward of r = 2.5R E . This has been observed in in-situ comparisons in previous studies and was attributed to attitude uncertainty issues which make it difficult to distinguish between the internal and external fields (e.g., [35]) as the former is very large close to the planet. If these are indeed observational errors, then perhaps these datapoints earthward of r = 2.5R E should be excised from the model database in future studies. Importantly for the context of this study, is that the errors of the Merged Model (dashed brown line) tend to follow the smaller of the other two across all radial bins. This confirms that the merging algorithm discussed in section 3.1 works as intended.

DISCUSSION AND CONCLUSION
In this paper we presented a new method of the empirical reconstruction of the magnetospheric substorms, which not only resolves the corresponding reconfiguration of the magnetotail but also resolves both westward and eastward currents in the inner magnetosphere that reflect the associated storm-type phenomena. In particular, it becomes possible to resolve and quantitatively evaluate the buildup of the storm-time plasma pressure, which was under-resolved in the previous substorm reconstructions [15]. The key to the solution of such a combined description of the inner magnetosphere is similar to the original kNN DM method [15,22], where a swarm of many synthetic probes neighboring the event of interest in the binning parameter space was used. Here similar virtual probes were used and The equivalent errors with no external field model is given for reference. combined from different versions of the SST19 model focused on the inner magnetosphere and the tail region. This merged resolution approach is also similar to coupled first-principles models of the magnetosphere [74,75,77,90], where the kinetic description of the inner magnetosphere is combined with the magnetohydrodynamic (MHD) description of the whole magnetosphere in global MHD models. An important advantage of the present empirical method, compared to the aforementioned combinations of the first-principle models, is that the empirical reconstructions weakly depend on the location of the coupling boundary and they do not require any special description of the interaction between inner and outer magnetosphere models. This approach may be further improved by using a more optimized merging boundary instead of the simple static cylindrical boundary used here. For instance, the merging boundary can be made dynamical, either redetermined for each time step or made a function of storm and substorm activity level. Further, the azimuthally symmetric boundary used here is suboptimal as it does not account for day-night asymmetries. This can be addressed by introducing a shift to the center of the cylinder or by instead using an elliptic cylinder.
In conclusion, the merged modeling technique using virtual observations effectively reconstructs regions of the magnetosphere possessing different spatial scales. This may also have utility in other DM and machine learning applications in which disparate density of data makes it difficult to model the system using a single resolution.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
GS was the primary author writing the original draft, led conceptualization, investigation, and formal analysis including software, data curation, visualization, and validation. MS assisted in writing by reviewing and editing, developed the methodology presented here including its conceptualization, software, visualizations, validation, managed the project, including project administration, funding acquisition, resources, and supervision. All authors contributed to the article and approved the submitted version.

FUNDING
This work was funded by NASA Grants 80NSSC19K0074, 80NSSC20K1271, and 80NSSC20K1787, as well as NSF Grants AGS-1702147 and AGS-1744269.