Robust Estimates of Spatiotemporal Variations in the Auroral Boundaries Derived From Global UV Imaging

The aurora often appears as an approximately oval shape surrounding the magnetic poles, and is a visible manifestation of the intricate coupling between the Earth's upper atmosphere and the near‐Earth space environment. While the average size of the auroral oval increases with geomagnetic activity, the instantaneous shape and size of the aurora is highly dynamic. The identification of auroral boundaries holds significant value in space physics, as it serves to define and differentiate regions within the magnetosphere connected to the aurora by magnetic field lines. In this work, we demonstrate a new method to estimate the spatiotemporal variations of the poleward and equatorward boundaries in global UV images. We apply our method, which is robust against outliers and occasional bad data, to 2.5 years of UV imagery from the Imager for Magnetopause‐to‐Aurora Global Exploration satellite. The resulting data set is compared to recently published boundaries based on the same images (Chisham et al., 2022, https://doi.org/10.1029/2022JA030622), and shown to give consistent results on average. Our data set reveals a root mean square boundary normal velocity of 149 m/s for the poleward boundary and 96 m/s for the equatorward boundary and the velocities are shown to be stronger on the nightside than on the dayside. Interestingly, our findings demonstrate an absence of correlation between the amount of open magnetic flux and the amount of flux enclosed within the auroral oval.


Introduction
The aurora is a visible manifestation of Earth's coupling to near-Earth space and is produced by charged particles precipitating into the upper atmosphere.These particles are usually confined to closed magnetic field lines that connect directly between the Northern and Southern Hemispheres.As a result, Earth's aurora often appears in an oval shape surrounding the magnetic poles.Based on photographs of the aurora, Feldstein and Starkov (1967) demonstrated how the average oval moves to lower latitudes and becomes wider as the level of geomagnetic activity increases.This average response to the level of activity has since been confirmed in multiple studies (e.g., Carbary et al., 2003;Han et al., 2020;Hu et al., 2017;Milan et al., 2010;Starkov, 1994;Y. Zhang & Paxton, 2008) and it has been reported that the radius of the oval is larger when the ring current is strong (Milan, 2009;Milan, Hutchinson, et al., 2009).Auroral boundary locations evince a large spread in these statistical studies, indicating that the size of the oval is highly variable.This dynamic nature of the shape and size of the auroral oval is observable both in individual events (e.g., Akasofu, 1964;Craven & Frank, 1987;Frank & Craven, 1988;Laundal, Østgaard, Snekvik, & Frey, 2010;Milan et al., 2003) and in superposed epoch analysis (e.g., Laundal, Østgaard, Frey, & Weygand, 2010;Mende et al., 2003;Milan, Grocott, et al., 2009).
Inside the auroral oval is a region of open magnetic field lines that extend into the solar wind.This is the polar cap and the boundary separating this region of open magnetic flux from the surrounding closed flux region is the open-closed boundary (OCB).The polar cap expands and contracts as magnetic flux is opened by dayside reconnection and closed by nightside reconnection.Net changes in polar cap magnetic flux, P, correspond to differences between opening and closure of flux, where Φ D and Φ N are the dayside and nightside reconnection rates, respectively.This is the expanding/contracting polar cap (ECPC) paradigm (Siscoe & Huang, 1985).If the instantaneous location of the entire OCB is known, P can be calculated using a model of the Earth's main magnetic field such as the International Geomagnetic Reference Field (Alken et al., 2021).The net difference between Φ D and Φ N can then be inferred by considering the temporal evolution.
A major advantage of global auroral images is the possibility to identify the entire OCB with high temporal resolution, since the OCB is approximately located at the poleward boundary (PB) of the auroral oval (J.B. Baker et al., 2000;Boakes et al., 2008;Carbary et al., 2003;Hubert et al., 2010;Kauristie et al., 1999;Longden et al., 2010).The method is most reliable during active periods with bright emissions, and less certain when the emissions are weak or during periods with emissions on open field lines within the polar cap, such as High-Latitude Dayside Aurora (Frey et al., 2003(Frey et al., , 2004;;Q.-H. Zhang et al., 2021) and proton aurora (Frey et al., 2002).Using subsequent PB locations obtained from auroral images to estimate dP/dt, Equation 1 can be used to give quantitative estimates of the global flux transport.Milan et al. (2012) identified intervals with no significant nightside reconnection (Φ N ≈ 0) and estimated Φ D by quantifying the rate of polar cap expansion.
Other studies have used empirical coupling functions (e.g., Milan et al., 2012;Newell et al., 2007) to estimate Φ D from solar wind observations and then used Equation 1 to find Φ N (Milan et al., 2007;Ohma et al., 2018).It is also possible to derive local reconnection rates by quantifying the amount of magnetic flux crossing the OCB in any given region based on observations of ionospheric plasma convection (Blanchard et al., 1996(Blanchard et al., , 1997;;Boudouridis et al., 2021;de la Beaujardiere et al., 1991;Hubert et al., 2006;Østgaard et al., 2005).Such estimates must be done in the reference frame of the OCB, which means that the boundary normal velocity must be identified.These quantitative approaches thus need reliable estimates of the change in the OCB location, not only its location at any given time.
Auroral images also provide the location of the equatorward boundary (EB) of the oval, and it has been shown that this boundary corresponds well with the equatorial boundary identified by particle detectors in low-Earth orbit (Chisham et al., 2022;Kauristie et al., 1999).This boundary is often associated with the earthward edge of the main plasma sheet (Galperin & Feldstein, 1991;Kauristie et al., 1999).Since particle precipitation significantly affects ionospheric conductivity, more accurate identification of the location of the auroral oval will significantly improve our ability to resolve ionospheric electrodynamics using data assimilation techniques (e.g., Laundal et al., 2022).Zhu et al. (2020) has demonstrated that binning precipitation measurements based on dynamic boundary locations yields more confined and intense electron energy flux and mean energy, which in turn affect conductance and convection estimates.
If the full PB and EB can be identified, the magnetic flux threaded by the auroral oval A can be found in the same way as P. A naive interpretation of the ECPC paradigm suggests that A increases with Φ N and decreases with Φ D : During nightside reconnection, open (dim) field lines are converted to closed (bright) field lines, and during dayside reconnection closed (bright) field lines are converted to open (dim) field lines.A corresponding interpretation in terms of magnetospheric dynamics is that the plasma sheet grows (in terms of magnetic flux) during nightside reconnection, followed by sunward transport and eventually destruction as the closed field lines of the dayside extension of the plasma sheet are opened by reconnection with the IMF.Decotte et al. (2023) found that this simplified fluid model could partially explain large-scale features of the auroral oval.However, the ECPC does not explain variations in the EB, which can be significant (e.g., Gjerloev et al., 2008;Mende et al., 2003).It also neglects several processes that excite or deplete the population of precipitating particles producing the auroral emissions, such as particle drifts, particle acceleration above the auroral zone, and wave-particle interactions (Coumans et al., 2002;Newell et al., 2009;Ni et al., 2016).
A number of approaches have been applied to identify the boundaries of the aurora.Gjerloev et al. (2007) identified the boundaries visually, while Hubert et al. (2006Hubert et al. ( , 2010) ) fitted a truncated Fourier series to the contours of the proton aurora.Several studies have made latitudinal intensity profiles of the auroral intensity through binning or curve fitting (typically Gaussian or double Gaussian) and determined boundary locations using global threshold values (e.g., Brittnacher et al., 1999;Frank & Craven, 1988;Mende et al., 2003;Ohma et al., 2018), a fixed fraction of the maximum intensity (e.g., Kauristie et al., 1999), variable thresholds based on noise levels (e.g., Gjerloev et al., 2008) or by using the width of the fitted functions (e.g., Boakes et al., 2008;Carbary et al., 2003;Chisham et al., 2022;Longden et al., 2010).More recently, image segmentation techniques based on clustering methods or deep learning have been applied to identify auroral boundaries (Ding et al., 2017;Han et al., 2020;Hu et al., 2017;Tian et al., 2020).
Besides global auroral imagers, there are several other ways to infer the location of auroral boundaries.The most accurate estimates are provided by particle measurements from low-Earth orbiting spacecraft (Kilcommons et al., 2017;Newell et al., 1991Newell et al., , 1996;;Sotirelis & Newell, 2000).The drawback is that only a local point estimate of the different boundaries is provided each time the spacecraft traverses the oval, which means that longitudinal or temporal coverage is limited.Space-based magnetometers can be used in a similar manner, as the frequency of magnetic perturbations is much higher within the oval compared to outside the oval (Xiong et al., 2014).Remote sensing instruments such as radars, photometers, and all-sky cameras can provide regional estimates of the OCB within the instrument's field-of-view (Aikio et al., 2006;K. B. Baker et al., 1995;Blanchard et al., 1995Blanchard et al., , 2001;;Chen et al., 2017;Chisham & Freeman, 2003;Chisham et al., 2007).In some regions, the ionospheric convection reversal boundary can serve as a proxy for the OCB (Hoque & Fenrich, 2018;Sotirelis et al., 2005), but this proxy can be misleading as there are also regions where the plasma flow is similar at either side of the OCB (Reistad et al., 2021).Instantaneous maps of the global field-aligned current systems can also be used to infer the OCB (Burrell et al., 2020;Clausen et al., 2013;Milan et al., 2015).However, a major drawback with this method is the inability to observe the immediate poleward contraction during substorms, which means that Φ N is severely underestimated during the unloading phase of substorms.Of available instruments, global images from high-Earth orbiting spacecraft are therefore the best option to infer the spatiotemporal evolution of magnetic flux in the auroral oval and polar cap.
In this paper, we present a novel method for creating a smooth and differentiable model of the auroral boundaries based on global auroral images.We use the model to investigate statistical distributions of the boundary normal velocities.We further address how the size of the auroral oval, measured in magnetic flux, and its temporal rate of change relate to the polar cap flux, solar wind forcing, and geomagnetic indices.We then explore a simple, dynamic model of the auroral area's evolution.
We describe the data sets in the next section.In Section 3, we outline how we detect auroral boundaries and how we use these boundaries to construct smooth models of the EB and PB variations in time and space.Finally, we describe how we create a data set of auroral boundaries including normal velocities.In Section 4, we present statistical distributions of the boundary locations and velocities and we explore how the size of the auroral oval relates to other geophysical quantities.We summarize our main results in Section 5.

Data
The Wideband Imaging Camera (WIC) on the Imager for Magnetopause-to-Aurora Global Exploration (IMAGE) mission provided images of the Lyman-Birge-Hopfield (LBH) emission band in the far ultraviolet range from 2000 to 2005 (Burch, 2000;Mende et al., 2020aMende et al., , 2000b)).During the first years of the mission, WIC had the entire auroral zone in the northern hemisphere within its field of view for up to 8 hours per orbit, providing images with 123 s cadence.We use images from the period 2000 to 2002 in this study.The images have been processed to account for sensitivity differences across the detector and over the lifespan of the mission and geo-located on Earth in Modified Apex coordinates (Richmond, 1995).We have also removed background emissions like dayglow based on solar zenith angle and viewing geometry (Ohma, Madelaire, Laundal, Reistad, Hatch, et al., 2023).Solar wind plasma and magnetic field data are from the OMNI 1-min data set, time shifted to Earth's bow shock (King & Papitashvili, 2005).This data set also includes commonly used geomagnetic indices, of which we will use the AE and SYM-H indices.To remove small data gaps and reduce the influence of spikes, we apply a centered, 10 min rolling mean to the OMNI data before we resample to the times of the auroral images by cubic spline interpolation.

Auroral Boundary Detection
To identify auroral boundaries, we use a two-step method.As an initial step we create latitudinal intensity profiles and identify auroral boundaries when the intensities reach different thresholds.This yields a set of boundaries as well as a measure of uncertainty, that we use in a second step, to construct smooth and differentiable spline representations of the boundaries' variation in time and space.The full approach is described in the following subsections.

Boundary Detection
As outlined in the introduction, a common approach to detect auroral boundaries from images is to create latitudinal intensity profiles and then infer boundaries from these profiles.We use a similar approach to this well tested method.We define 24 latitudinal meridians with 1 hr separation.We define the intensity at each point along the meridians as the median of all pixels located within a 300-km radius of that point.Three example profiles are shown in Figures 1a-1c and the solid latitudinal lines in Figure 1d indicate the MLT of these profiles.The black ring in Figure 1d illustrates the circle within which the median pixel intensity is calculated for the location indicated by the black dot.This location is also indicated by a black dot in Figure 1c.The shaded region at lower latitudes in Figure 1d are not considered when detecting boundaries.The aurora rarely extends into this shaded region (Milan et al., 2010;Weygand et al., 2023).
To determine auroral boundaries, we identify the first and last point where the intensity crosses specific thresholds.From the intensity profiles in Figure 1, we see that the transition from non-aurora to aurora can be sharp or vague, and that multiple peaks can exist.When there is a single peak and the transition is sharp, the boundary is well defined and only weakly dependent on the threshold level.When several peaks exist and/or the transition is more gradual, the boundary is ambiguous and strongly dependent on the threshold level.To quantify the uncertainty, we use thresholds from 50 to 200 counts with increments of 5 when we determine the initial boundary locations.Sharp auroral boundaries then result in closely clustered boundary locations, whereas more ambiguous boundaries will have a larger spread.Identified boundary locations using a threshold of 50 (light red), 125 (red), and 200 (dark red) counts are displayed in Figure 1.

Boundary Modeling
Well determined boundaries from individual images are sufficient for many applications.However, as outlined in the introduction, some key applications rely on knowledge of how the boundary locations evolve with time.This is challenging to do coherently, as both the spatial and temporal derivatives determined from images are quite noisy.The signal to noise ratio can be low, especially at the dayside, and the exact locations of the boundaries can be ambiguous due to weak intensity gradients.There can also be strong emissions outside the oval, both on open field lines within the polar cap and at sub-auroral latitudes.In addition, uncertainty when geo-locating the images causes artificial movement or wobbling between images.This is a known issue with the IMAGE satellite, as a broken antenna probably affected the spin of the spacecraft (Green & Reinisch, 2003).Differentiating the boundaries thus requires outlier removal combined with heavy spatial and temporal smoothing and interpolation.We aim to produce models of both the poleward and equatorward boundaries, based on the initially identified boundaries.These models must cover all local times and be smooth and differentiable in both space and time and be minimally influenced by outliers.To achieve this, we use different approaches for the two boundaries, described in more detail below.

Equatorward Boundary Model
To represent the co-latitude of the EB, we seek a function τ e of magnetic local time ϕ and time t that represents the smoothed boundary co-latitude, represented in terms of B-splines (de Boor, 2001).The B-spline representation along ϕ is periodic to ensure smoothness at all local times.In order to guarantee that τ e is never negative, we model it as the exponential of a sum of B-splines so that where τ e = expτ′ e and a ij are model coefficients to be determined.The B-splines B i and B j are defined by an order and by the location of knot points, and sum up to represent any spline function.We have used an order of 3. Since the EB is fairly smooth, we set the spatial knots at 0, 6, 12, and 18 MLT, corresponding to four B-splines.This limits the spatial variation allowed in our model.We place temporal knots every 10th minute.
Equation 2 can be used to construct a matrix equation that relates the set of model coefficients, comprising the vector m e , to the set of initially determined boundary locations, whose logarithms comprise the vector θ′ e : (3) The design matrix G e is the linear relationship between θ′ e and m e described by the spatiotemporal spline function, Equation 2.
The coefficient vector m e can now be found by inversion.We use iteratively re-weighted regularized least squares (e.g., Aster et al., 2013;Madelaire et al., 2022).For the (i + 1)'th iteration, the solution vector m e is where G e is the design matrix and θ′ e is a vector comprised of the logarithms of the observed equatorward boundaries.The diagonal matrix W i contains Huber weights, based on the i'th iteration, to minimize the influence of outliers (Huber, 1964;Huber & Ronchetti, 2009), and is updated in each iteration until the relative change in the model norm is less than 0.1%.The matrices L t and L ϕ provide the temporal and longitudinal gradients of the model boundary, and the regularization parameters λ t and λ ϕ determines the strength of regularization along these two dimensions.Increasing λ t and λ ϕ give a smoother model in time and space, respectively, at the expense of fitting the initial boundary estimates.The constants a t and a ϕ are simply scaling factors to make the regularization invariant with respect to the number of input boundaries.Following Laundal et al. (2022), we set a t and a ϕ equal to the median value of the diagonal of G T e G e divided by the median value of the diagonal of L T t L t and L T ϕ L ϕ , respectively.This scaling ensures that the optimal regularization parameters are always of similar magnitudes.
To determine the value of the two regularization parameters, we use L-curve analysis (Hansen, 1992).For no regularization in the longitudinal direction, we find λ t = 1 to be a reasonable compromise between the semi-norm ‖L t m e ‖ 2 and the model misfit.Using λ t = 1, we find λ ϕ = 10 2 to be a reasonable value for the regularization of longitudinal gradients.We obtain the same values if we reverse the order to find λ t and λ ϕ .These numbers are based on experiments with several time intervals and are kept constant in our data set.
When the model coefficients are determined, EB locations τ e for any MLT and universal time constrained by the model are obtained using Equation 2.

Poleward Boundary Model
The PB co-latitude θ p should be a positive value less than τ e .To enforce this constraint, we define a fraction ρ such that which should be between 0 and 1. B-spline representations are not limited to this range, so we define θ′ p as a function of ρ such that values between 0 and 1 are mapped to R: where θ′ p is represented in terms of B-splines.The model of the PB of the auroral oval should be able to capture the dynamic and structured behavior on the nightside.We, therefore, apply a 2-hr knot separation between 18 and 6 MLT while the separation over 6-18 MLT is similar to the EB model.We do not add more knots at the dayside because structures here are often caused by features on open field lines inside the polar cap, which we try to avoid.
Similar to the EB, the vector of coefficients a ij for the PB model, m p , is found by iteratively computing the least squares solution of the regularized inverse problem where θ′ p is a vector comprised by observed values of θ′ p .Again, we use L-curve analysis to determine the two regularization parameters and again find λ t = 1 and λ ϕ = 10 2 to be reasonable values.
When the model parameters are determined, we have a model of the primed PB τ′ p (ϕ,t) as a sum of B-splines.To recover the actual (unprimed) boundary colatitude, we invoke the inverse operation:

Partial Derivatives
To quantify the amount of magnetic flux crossing the auroral boundaries, the boundary normal velocities must be known.The derivation in Appendix A demonstrates how the normal velocity depends on the derivative of τ with respect to ϕ and t.The derivative of B-splines are uniquely and analytically defined by the model coefficients and knot locations, and are new B-splines with one order less than the original B-splines.We thus have analytic representations of the necessary derivatives for τ′ e and τ′ p , such that it is not necessary to approximate the derivatives via finite differences.For the EB, the transformation to unprimed coordinates gives The longitudinal and latitudinal velocity components of the boundary normal are then given by Equations A5 and A6, respectively.The boundary normal velocity is , where we define positive to be toward the magnetic pole and negative to be toward the equator.aurora well and captures most large-scale structures.However, the predicted location does not fit perfectly everywhere as expected from the spline representation and regularization discussed above.We have created boundary models for each orbit of IMAGE from May 2000 to December 2002, inclusive, yielding boundaries from roughly 391,000 different times.Here one must bear in mind that, due to regularization based on gradients in space and time, our inversion scheme is capable of generating a boundary model at all locations and times based on a single input boundary.The model can therefore be poor even though it fits the input data well, for instance, if the data coverage is limited or if the signal to noise is very low.To be included in the statistical survey presented in Section 4, we require the boundary at a given time satisfies the following criteria:

Boundary Model Quality Evaluation
• Data coverage: The polar cap must be within the field of view of WIC.We have defined this by requiring than the region poleward of colatitude θ = 30 + 10 cos ϕ, where ϕ is MLT in radians, is observed with a viewing angle less than 65°.We use a rolling window of five images centered at the time in question, to avoid discarding images when data is missing from individual images even though the combined coverage is good.• Detected boundaries: We require that initial boundaries are detected in more than 12 MLT sectors.Again we use a rolling window of five images for the same reason as above.• Auroral intensity: The mean intensity within the defined oval must be larger that the mean intensity plus two times the standard deviation in both the region poleward of the oval and equatorward of the oval.• Boundary uncertainty: The uncertainty associated with the detection threshold level should not be too large.
We have enforced this by requiring that 75% of both the poleward and EB have an uncertainty below 1.5°T hese criteria are to some degree overlapping, but in general, ensure that the full oval is observable in terms of data coverage and intensity with reasonable certainty.Applying these criteria, we are left with 65,357 time steps with good boundaries (17%).

Boundary Validation
When developing an automated method as described here, it is necessary to make choices that cannot always be determined completely objectively.To make the latitudinal intensity profiles, we take the median intensity in a 300 km radius of each point.By visual inspection of many intensity profiles, we find this to be a reasonable compromise between reducing random noise and keeping sharp intensity gradients.Using this radius the profiles also become independent equatorward of 70°latitude, a typical location of auroral emissions.We have preferred this approach over binning the images in circular sectors as it ensures that similar numbers of pixels contribute at all latitudes.The threshold values used to detect boundaries are similar to threshold values used for the WIC camera in earlier studies (Ohma et al., 2018;Østgaard et al., 2005).We also typically obtain threshold values between 100 and 150 when we identify the Otsu threshold from background-corrected images, a standard method to separate foreground from background in gray-scale images (Otsu, 1979).
When we construct the boundary models, we must decide on knot locations.We have chosen to have as few as possible while still fitting most of the large-scale structures.This conservative approach ensures that the model is not too flexible.This prevents over-fitting, but also means that we can not expect to fully resolve all structures.For events with a high signal to noise ratio, more knots could be added to fit smaller structures.With the knot location determined, the regularization is objectively identified using L-curve analysis.
We have tried to make reasonable choices in constructing our data set, but the specific choices we make are, again, to some degree heuristic.We have therefore validated our approach by comparing our boundaries to the recently updated WIC boundaries identified by the British Antarctic Survey (BAS) as described by Chisham et al. (2022), which builds on the Longden et al. (2010) method.They divide individual images into circular sectors and fit a single Gaussian or a double Gaussian to the latitudinal intensity profiles (plus a polynomial to capture background emissions).They select the best fit and set the boundary location as the full-width-at-half-maximum distance from the appropriate peak.They also apply several criteria to discard ill-determined boundaries.Finally, Chisham et al. (2022) calibrate the auroral boundaries with particle boundaries from the Defense Meteorological Satellite Program (DMSP) spacecrafts (Kilcommons et al., 2017) and provide a statistical correction.
The comparison between the two boundary data sets is displayed in Figure 4, where only conjunctions in time and MLT are included.The polar plot in Figure 4a shows the median boundary location for our data set (blue), the uncorrected BAS data set (orange), and the DMSP-corrected BAS data set (green).The gray sectors indicate regions where this correction is not supported by data.As evident from the figure, our median boundaries are in reasonable agreement with the BAS boundaries, but our oval is consistently more narrow.It is also interesting to note the good match between our boundaries and the corrected BAS boundaries in the dusk sector.This indicates that our median boundary is close to the typical particle boundary without any correction.Our PB is also closer to the corrected BAS boundary in the pre-noon sector.However, the EB is close to the uncorrected BAS boundary, indicating a misalignment here.
The locations of the BAS boundaries relative to our boundaries are displayed in Figures 4b-4g, where Figures 4b  and 4c show the dayside (6-18 MLT), Figures 4d and 4e show the nightside (18-6 MLT) and Figures 4f and 4g show the highlighted region where DMSP data is available for correcting the BAS boundaries.Our boundaries tend to be closest to the corrected BAS boundaries (green), supporting that our data set is close to the particle boundaries on average.Note that a large fraction of the difference between our data and the BAS data is due to the wobbling of the geo-located WIC images, as our spatiotemporal model greatly reduces this artificial boundary movement.
The comparison demonstrates that our method consistently aligns with the calibrated BAS boundaries, which originate from the same data set.Unlike the BAS method, however, ours also provides robust estimates of the spatiotemporal variations of the boundary, enabling the calculation of boundary motion at any local time.While global UV image-based boundary estimates have been extensively compared with more precise OCB estimates from in situ particle measurements (Boakes et al., 2008;Carbary et al., 2003;Longden et al., 2010), the validation of boundary motion estimates is still pending.Due to the lack of alternative data sets representing this motion accurately, any validation would require simulated global aurora and UV images, which falls outside the scope of this study.

Results
In the next sections, we present distributions of boundary locations based on our boundary data set.We then examine how the amount of magnetic flux in the auroral oval relates to the polar cap flux, solar wind forcing, and geomagnetic activity.Finally, we explore the performance of a very simplified fluid model of the auroral area.

Boundary Locations and Velocities
Figure 5 displays normalized distributions of the boundary locations.Figure 5a shows the MLT-dependent distribution of the PB, Figure 5b shows the MLT-dependent distribution of the EB, and Figure 5c shows the combined distributions of both boundaries.The white lines indicate the upper and lower quartiles.The median PB changes from ∼70°near midnight to ∼76°near noon while the median EB changes from ∼61°to ∼68°, with an interquartile range of ∼3.8°and ∼4.4°.The EB distribution thus has a slightly larger spread than the PB distribution.
The corresponding distributions for the boundary normal velocities are displayed in Figure 6, showing MLTdependent PB and EB velocities in Figures 6a and 6b and the combined distributions in Figure 6c.The PB velocity distribution is wider than the EB distribution and has a slight skew toward positive values.Hence, the PB typically moves slightly faster toward the magnetic pole than toward equator; and the fastest velocities are found in the midnight sector, consistent with the explosive nature of substorm unloading.Additionally, a decrease and increase of the same amount of magnetic flux will elicit different velocities as the flux density changes with latitude.Therefore, poleward velocities will statistically be larger than equatorward velocities.Furthermore, this effect will increase with latitude as the ratio between magnetic flux density poleward and equatorward of the boundary decreases with latitude.The distribution is also considerably wider at the nightside compared to the dayside, again consistent with rapid nightside dynamics.The EB velocity distribution has a similar spread at all MLTs, but a small peak near 20 MLT.The Root Mean Squared speeds of the PB and EB when considering all MLTs are 149 and 96 m/s, respectively.For the PB, the maximum occurs at midnight (201 m/s) and the minimum at 10 hr MLT (124 m/s).For the EB, the maximum occurs at 19 hr MLT (111 m/s) and the minimum at 3 hr MLT (80 m/s).

Magnetic Flux in the Auroral Oval
It is common to measure the polar cap area in terms of magnetic flux, as changes are then directly related to the dayside and nightside reconnection rates (magnetic flux per time).It is therefore convenient to also consider the auroral oval's area in terms of magnetic flux, which we refer to as auroral magnetic flux.Since our boundaries are defined in Modified Apex coordinates, we can treat the Earth's magnetic field as a dipole.According to Richmond (1995), the magnetic flux per unit interval of modified apex latitude λ and longitude ϕ is where R = R E + h R is the modified apex reference radius,  sin I = and B e 3 = |B|/ D. |B| is the magnetic field strength and D is a geometric factor that is 1 at the reference height h = h R for dipole magnetic fields, but otherwise varies depending on the distortion introduced by non-dipole terms.
dF, ϕ, and λ are by definition constant along magnetic field lines, which implies that B e 3 is also constant along magnetic field lines.Points at high latitudes map along magnetic field lines to large distances, where the Earth's magnetic field is dipolar; and if a dipole can be used at high altitudes, it can also be used closer to Earth, since the magnetic flux per unit interval of λ and ϕ is constant along magnetic field lines regardless of the field structure.Equation 11 thus essentially translates between a realistic field configuration and a dipole, so that we can use a simple equation and still take the full IGRF into account by virtue of the definitions of ϕ and λ.
Based on this, we use dipole equations to replace B e 3 (evaluated at radius R): where B 0 is a reference magnetic field, equal to the norm of a vector composed of the dipole gauss coefficients in a spherical harmonic representation of the magnetic field (e.g., Laundal & Richmond, 2017).Inserting this in Equation 11we get Integrating this above a boundary defined by the curve λ(ϕ), we get The polar cap magnetic flux P is found by setting λ(ϕ) = π/2 τ p (ϕ) in Equation 15, and the auroral oval magnetic flux A is found by setting λ(ϕ) = π/2 τ e (ϕ) and subtracting P.
In Figure 7, we compare A and the rate of change dA/dt to P and dP/dt.Figure 7a displays the distribution of A (green) and P (orange) from concurrent images.While the two distributions peak in the same region, A has a significantly larger spread.The size of the auroral oval is therefore more variable than the size of the polar cap.The relation between A and P is explored in Figure 7b.Surprisingly, the correlation between A and P is only 0.09.
Figure 7c displays the distributions of dA/dt (green) and dP/dt (orange) and Figure 7d displays the relation between the two.Both distributions are centered near zero with the auroral distribution having the larger spread.The auroral distribution is also slightly skewed toward positive values, consistent with the rapid expansion known to occur during substorm expansion.As seen in the latter figure, the two values are anti-correlated with r = 0.65.This is expected, as contraction of the polar cap area increases the auroral area.Note that as both A and P depend on the PB, both the real and erroneous movement of this boundary contribute to the correlation.We have tried to minimize erroneous movement through regularization, but the correlation coefficient should anyway be considered as a lower limit (the anti-correlation is probably smaller).Also note that the slope in Figure 7d is about 2, much larger than the one-to-one relation that would exist if the EB was held fixed.We will return to this in Section 4.3.
In Figures 8a and 8b we compare A and dA/dt to Φ D , estimated from solar wind data using the empirical coupling function provided by Milan et al. (2012).We see that the auroral oval area is to some extent controlled by the dayside reconnection rate, but the correlation coefficient is only 0.45.The rate of change of A is practically uncorrelated with Φ D .Figures 8c and 8e compare A and to the geomagnetic indices AE and SYM-H, while Figures 8d and 8f compare their derivatives.Figure 8c shows that there is a strong relation between the AE index (which quantifies the combined strength of the eastward and westward electrojets) and A, which implies that the AE index is a good predictor of the auroral size, in addition to being a good predictor of the auroral power (Newell & Gjerloev, 2011).This result is consistent with earlier studies showing increased average auroral area for higher geomagnetic activity (Feldstein & Starkov, 1967;Starkov, 1994;Y. Zhang & Paxton, 2008).However, the rate of change of the AE index and dA/dt only have a weak correlation of about 0.3 (Figure 8d).The AE index thus has limited predictive power with regard to changes in the oval size.
Figure 8e shows that the SYM-H index (quantifying the strength of the ring current) is also clearly correlated with the auroral magnetic flux, but with a considerably lower correlation coefficient.The temporal derivatives (Figure 8f) appear to be completely uncorrelated.

A Simple Model of Auroral Magnetic Flux
The ECPC paradigm tells us that changes in polar cap flux P are equal to the difference between flux opened by dayside reconnection and flux closed by nightside reconnection.Figure 9 illustrates this concept, based on a similar sketch presented by Siscoe and Huang (1985), where we have added the auroral oval (green) outside the PB (red).Figure 9a shows a sketch of how an impulse of nightside reconnection leads to a contraction of the polar cap and an associated expansion of the auroral oval.Expansion of the auroral oval occurs during the unloading phase of substorms and has been an integral part of all substorm models since first described by Akasofu (1964).Based on this conception, it is natural to assume that the change in auroral magnetic flux dA/dt varies with Φ N .Figure 9b illustrates how an impulse of dayside reconnection leads to an expansion of the polar cap on the dayside.If this region contained aurora prior to the impulse, the trapped particles producing the aurora might be lost as the field lines change from a closed to an open configuration.The auroral magnetic flux would then decrease with a rate proportional to the dayside reconnection rate Φ D (Decotte et al., 2023).
In Figure 10 we explore these ideas in relation to our data set.Figure 10a displays dA/dt versus Φ N = Φ D dP/dt (Equation 1), where Φ D is estimated using solar wind data and the Milan et al. (2012) equation.There is, as expected, a clear relation between the two.Interestingly, the correlation is lower than the correlation found between dA/dt and dP/dt directly (Figure 7d).By including Φ D , we would have expected to account for polar cap  expansion not directly affecting the auroral size and thus increasing the correlation, but this does not happen.It is also interesting to note that we still find a slope that is larger than one.This indicates that nightside reconnection excites precipitation in a larger region than the auroral bulge directly associated with the reconfiguration of field lines.
We explore the influence of dayside reconnection on dA/dt in Figure 10b.Here we only consider the auroral magnetic flux in a 0.1-hr wide MLT sector at noon to isolate dayside effects.We find no apparent correlation, consistent with the result for the full auroral oval displayed in Figure 8b.Although appealing from a conceptual perspective, we must therefore conclude that dayside reconnection has no direct influence on the auroral area on a systematic level, at least for aurora of the energy range and brightness observed by WIC.
If we assume that nightside reconnection is the only process that produces auroral magnetic flux, we can formulate a simple model of the evolution of the auroral area based on Figure 10a: where L is some loss term.It has been shown earlier that the strength of the field-aligned Birkeland currents have an exponential decay when the solar wind forcing is low, with a decay time of about 1 hr (Moretto et al., 2018).A mean decay time of about 14 ± 4 hr have been found for the strength of the ring current following geomagnetic storms (Dasso et al., 2002).To explore if the aurora follows a similar trend, we plot L versus A in Figure 10c.The two quantities are correlated, but there is also considerable spread.A simple linear fit yields an exponential decay time of 3.6 hr.
By combining the above results, we obtain a simple dynamic model of the auroral magnetic flux Correlating the two sides of this equation we obtain a linear correlation of 0.64.While this is better than the correlation obtained using only geomagnetic indices, it is not better than the direct correlation in Figure 7d.Thus, this model does not perform any better than just considering the motion of the poleward aurora boundary directly.The model is also rather impractical as both Φ N and A must be known.Since we can only explain about 40% of the variation, the model appears to be too simple.
To illustrate the temporal evolution more clearly, Figure 11 displays the evolution of dA/dt (green) and Φ N (orange) in the left column and L (green) and A/3.6 (purple) in the right column for three consecutive orbits (28-29 August 2000).These orbits occurred during a period of high solar wind speed and high geomagnetic activity, with a visible oval and only minor gaps in the solar wind observations.From the left column, we see that the two curves often follow each other, but there are also clear deviations.From the right column, we observe that the exponential decay (purple) is not able to pick up the rapid fluctuations observed in the loss term (green).At best, it is able to describe some slow, average decay.

Discussion and Summary
In this work, we have presented a novel method to create a spatiotemporal model of the poleward and equatorward boundaries of the auroral oval.The model is smooth and differentiable in both longitude and time.We have used our own method to detect initial boundaries from auroral images to construct the models, but any boundary data set from any source can be used (the data can be regular or irregular in both space and time).Both the processed data set and the code is publicly available (Ohma et al., 2023a(Ohma et al., , 2023b(Ohma et al., , 2023c)).Since the model approach presented here provides the boundary normal velocity directly, it is ideally suited to be used in conjunction with ionospheric data assimilation techniques such as the Assimilative Mapping of Ionospheric Electrodynamics procedure (Lu, 2017;Richmond & Kamide, 1988) and the Local mapping of polar ionospheric electrodynamics (Lompe) technique (Hovland et al., 2022;Laundal et al., 2022) to estimate regional or global reconnection rates.Using the PB model and estimates of the dayside reconnection rate, we also obtain estimates of the instantaneous nightside reconnection rate.This could be used in future studies to construct scaled coupling functions of Φ N from geomagnetic indices to be used when direct estimates are unavailable.
We have used B-splines as basis functions in both space and time when constructing boundary models.However, truncated Fourier series have commonly been applied to get a mathematical representation of the average aurora (e.g., Gjerloev et al., 2008;Holzworth & Meng, 1975;Hubert et al., 2006).It is straightforward to replace the spatial B-spline B i in Equation 2with a Fourier series and this approach is also provided in the code repository.
We have preferred the B-spline approach as it is easier to regularize since the model coefficients have similar amplitudes, as opposed to the Fourier series where the magnitude decreases with increasing order.However, both methods are capable of fitting the auroral oval well.
We presented distributions of the boundary locations and velocities in Section 4.1.The boundary locations are in reasonable agreement with earlier results (e.g., Hu et al., 2017), but an exact match is not expected as different methods and time periods favor different geomagnetic activity levels.We are not aware of earlier studies showing velocity distributions as displayed in Figure 6.However, several studies have examined the boundary velocities during substorms.Craven and Frank (1987) examined two isolated substorms and found an average poleward motion of 230 m/s with a maximum of 1,000 m/s during an intensification.Gjerloev et al. (2007) combined 116 isolated substorms and report a poleward motion of 660 m/s, but expansion rates above 1 km/s for the east and west ends of the bulge.Mende et al. (2003) found an average of ∼4°poleward movement at the onset location during the first 10 min of the expansion phase, corresponding to about 740 m/s, but considerably slower changes at other locations and epoch times of the averaged substorms.While the averages found by Craven and Frank (1987) agree well with the tail of our PB velocity distribution, the most extreme speeds found by Gjerloev et al. ( 2007) and Mende et al. (2003) are rare in our statistics.There are several reasons why our long-term distributions usually show smaller magnitudes: First, the occurrence frequency of such high speeds is low, as they only occur just after an onset or intensification.Second, Gjerloev et al. (2007) and Mende et al. (2003) bin their boundaries relative to the substorm location, but we bin the velocities in a fixed MLT grid.High, substorm related velocities are thus smeared out on the nightside.Third, our temporal knot separation of 10 min combined with the regularization needed to make our model robust and unaffected by wobbling means that we do not have the time resolution to resolve these very high, short-lived boundary velocities.Our distributions should therefore be considered as slightly conservative estimates of the real distributions and the tails are likely slightly heavier.Nevertheless, Figure 3, which shows an example event that took place during active conditions (SYM-H index ≈ 50 nT and AL fluctuating between 100 and 1,000 nT) shows poleward velocities in excess of 500 m/s, in good agreement with the event studies cited above.
In Sections 4.2 and 4.3 we explored how different quantities could predict the temporal evolution of the auroral magnetic flux, dA/dt.Specifically, we investigated a naive extension of the Expanding-Contracting Polar Cap paradigm, in which the auroral area is expanded via nightside reconnection and reduced via dayside reconnection (Decotte et al., 2023).For the aurora observed with the IMAGE WIC camera, our analysis shows that this idea is too simple, indicating that kinetic effects, like acceleration and wave-particle interactions, dominate the rate of change.Based on the limited statistics investigated here, the best model for the temporal evolution of dA/dt is 2 ⋅ dP/dt (Figure 7d), indicating that nightside reconnection causes auroral emissions also in regions outside the bulge region.This is in qualitative agreement with the superposed epoch analysis of the PB and EB, which demonstrated that the average EB moves equatorward as the average PB moves poleward in the expansion phase (Gjerloev et al., 2008;Mende et al., 2003).Gjerloev et al. (2008) provide a mathematical representation of the average PB and EB relative to substorm onset.By assuming a dipole field, we find dA/dt ≈ 1.6 ⋅ dP/dt based on these boundaries, providing quantitative support for our result.Gjerloev et al. (2008) also point out that the auroral intensity in the main oval and the bulge region decay at different rates in the recovery phase, indicating that several decay rates exist.
We found no evidence for dayside reconnection reducing the spatial extent of the oval.However, it is important to recognize that boundary detection is most uncertain around noon (Chisham et al., 2022).This means that there is a large degree of smoothing and interpolation in this specific region, which might mask such an effect.In addition, WIC is not sensitive to very soft precipitation, so we do not see the full diffuse aurora.It is therefore possible that this effect exists, but we find no support for it in this analysis.
The auroral oval can also expand during prolonged periods of northward-directed magnetic field in the solar wind.This is called the horse collar aurora (Hones et al., 1989;Murphree et al., 1982), and it has been suggested to be the result of dual lobe reconnection (Milan et al., 2020).This means that reconnection at the magnetopause could directly expand the oval, which we have not included in our simple model.However, the horse-collar aurora is usually below the signal to noise ratio in WIC, and is therefore not expected to be included in our analysis.
The aurora is the most prominent and visually striking manifestation of the complex interplay between the solar wind and Earth's magnetosphere.It is within the auroral oval that the bulk of energy dissipation resulting from this coupling occurs, making it a focal point for investigating hazardous space weather effects.However, despite its paramount importance, the spatial extent of the auroral oval exhibits significant variability, and the underlying mechanisms driving this variability remain poorly understood.The lack of recent high-quality global observations hampers our progress, as the most up-to-date data set available, used here, is over two decades old.Although the upcoming ESA-CAS SMILE mission with its UV camera holds promise for temporary improvement, the need for dedicated missions targeting global auroral observations persists (Branduardi-Raymont et al., 2021).The methodology proposed in this study and in the work by Ohma, Madelaire, Laundal, Reistad, Hatch, et al. (2023) can be readily applied to similar data sets, offering valuable insights into auroral dynamics.
where θ and ϕ are the colatitude and longitude at a time t.The material derivative of f(θ, ϕ, t) is where v is the velocity of some point on the surface.
Let us require that v always be in the direction normal to the surface, where the normal vector (without normalization) is Then v = 〈0,v θ ,v ϕ 〉 , and

Figure 1 .
Figure 1.Illustration of the boundary detection algorithm.Panels (a-c) show intensity profiles at 7.5, 13.5, and 20.5 hr MLT and panel (d) displays the MLT locations of the three profiles.The intensity at each point along the meridians is the median intensity of all pixels located within a 300-km radius of that point.The black ring in panel (d) indicates the spatial extent using this radius.Boundaries when using thresholds of 50 counts (light red), 125 counts (red), and 200 counts (dark red) are indicated in all panels.The shaded regions are not considered when detecting boundaries.

Figure 2
Figure 2 displays 16 images from 21:38:49 on 28 August 2000 to 05:20:26 on 29 August 2000.The pink and cyan dots indicate the initially detected PB and EB and the red and blue lines indicate the PB and EB models.The shaded regions around the two lines indicate the model uncertainty related to threshold levels in the initial detection.This uncertainty is estimated by creating an individual boundary model for each threshold level separately and taking the standard deviation of these boundary models.The boundary model generally fits the

Figure 2 .
Figure 2. Sixteen auroral images from 28 to 29 August 2000.The pink and cyan dots indicate the initially detected poleward boundary (PB) and equatorward boundary (EB), respectively, and the solid red and blue lines indicate the modeled PB and EB, respectively.The shaded regions indicate the uncertainty of the boundary models.

Figure 3
Figure 3 further displays the model performance, based on all images from 21:18 on 28 August 2000 to 08:04 on 29 August 2000 (full orbital pass).The medians of the initially detected PB are shown as color in Figure 3a and the PB model is shown in Figure 3b.The region between the two horizontal red lines has global data coverage in the northern polar region, to be defined below.The figure demonstrates how the model smooths and interpolates the input boundaries.The uncertainty related to threshold values is indicated in Figure3cand is generally low when the full auroral region is seen by the camera for this specific event.Figure3ddisplays the boundary normal velocity extracted from the model, which matches the gradients in the boundary latitude.

Figure 3 .
Figure 3. Evolution of the poleward boundary (PB) locations from 28 to 29 August 2000.(a) Median location of the initially detected PB.(b) Model location of the PB.(c) Uncertainty related to threshold levels and (d) boundary normal velocity.The region between the two horizontal red lines in all panels are identified as having global coverage.

Figure 4 .
Figure 4. Comparison between our modeled boundaries (blue), British Antarctic Survey (BAS) boundaries (orange), and Defense Meteorological Satellite Programcorrected BAS boundaries (green).Panel (a) displays the median boundary locations.The top row (panels (b, d, and f)) show the comparisons for the equatorward boundary, while the bottom row (panels (c, e, and g)) show the comparisons for the poleward boundary.The different columns correspond to different sectors.In each panel, the x-axis shows the relative difference between boundaries from our model with corresponding BAS boundaries.

Figure 5 .
Figure 5. Distributions of boundary locations.(a) Normalized poleward boundary (PB) locations on a logarithmic scale and (b) distributions of equatorward boundary (EB) locations on a logarithmic scale.The white lines indicate the upper and lower quartiles.(c) Distributions of PB (red) and EB (blue) locations for all MLTs.

Figure 6 .
Figure 6.Distributions of boundary velocities.(a) Normalized poleward boundary (PB) velocities on a logarithmic scale and (b) distributions of equatorward boundary (EB) velocities on a logarithmic scale.The white lines indicate the upper and lower quartiles.(c) Distributions of PB (red) and EB (blue) velocities for all MLTs.

Figure 7 .
Figure 7. (a) Distributions of polar cap flux P (orange) and auroral magnetic flux A (green) from concurrent images.(b) Polar cap flux versus auroral magnetic flux, where r is the Pearson correlation coefficient.Distributions of dP/dt (orange) and dA/dt (green) from concurrent images.(b) dP/dt versus dA/dt, where r is the Pearson correlation coefficient.

Figure 8 .
Figure 8.(a) Auroral magnetic flux A versus dayside reconnection Φ D .(b) dA/dt versus Φ D .(c) A versus the AE index and (d) dA/dt versus the temporal change in the AE index.(e) A versus the SYM-H index and (d) dA/dt versus the temporal change in the SYM-H index.The Pearson correlation coefficient r is indicated in each panel.

Figure 9 .
Figure 9. Conceptual sketch of how the expanding/contracting polar cap could influence the surrounding auroral oval.(a) An impulse of nightside reconnection contracts the polar cap area and expands the auroral area.(b) An impulse of dayside reconnection expands the polar cap area and contracts the auroral area.

Figure 10 .
Figure 10.(a) Relation between the temporal change in auroral magnetic flux dA/dt and the nightside reconnection rate Φ N .The dashed, pink line indicates the slope.(b) Relation between the temporal change in auroral magnetic flux in a 0.1-hr MLT bins at noon dA 12 /dt and the dayside reconnection rate Φ D .(c) Relation between the loss term L = 1.7ΦN dA/dt and the amount of auroral magnetic flux A. The dashed, pink line indicates the slope.The Pearson correlation coefficient r is indicated in all panels.

Figure 11 .
Figure 11.The left column displays the evolution of dA/dt (green) and Φ N (orange) during three orbits on 28-29 August 2000.The right column displays the observed loss L (green) and potential exponential decay related to the amount of auroral magnetic flux (purple) for the same orbits.
We then insert Equation A5 into Equation A4 and solve for v θ to obtain