Toward the Reconstruction of Substorm‐Related Dynamical Pattern of the Radiowave Auroral Absorption

In addition to existing empirical models describing the average distributions of energetic electron precipitation into the auroral ionosphere at different activity levels, we develop and test a semiempirical approach to construct dynamical models describing the recurrent features of spatiotemporal development of auroral absorption in the ionosphere during individual substorms. Its key ingredients are (a) usage of linear prediction filter technique to extract from riometer data the response function to the injection of unit magnitude and (b) characterization of injection parameters by midlatitude magnetic variations caused by the substorm current wedge. Using global riometer network we test the method performance for stations in the middle of auroral zone (at corrected geomagnetic latitudes of 65–67°) where generally the absorption amplitude is largest. In this paper we use the midlatitude positive bay index, recently developed by X. Chu and R. McPherron, to drive the model. We evaluate the model performance, discuss the dynamical properties of energetic electron precipitation as revealed by the linear prediction filter response function analyses, and finally, we discuss possible future improvements of this method intended for both science and applications.

Historically, one of most simple and useful monitoring tools for enhanced D region ionization were the riometers (relative ionospheric opacity meter; Little & Leinbach, 1959) which utilize the cosmic radio noise absorption (CNA) in the auroral ionosphere (also known as the auroral absorption) to evaluate the enhanced D region electron concentration caused by EEP (see, e.g., Hargreaves (1969) for a review of method and its early results). During substorms the absorption events in riometer records often have the bay-like shape and thus are frequently called as the CNA bays. Observationally, the CNA magnitude was shown to be closely related to the precipitated electron fluxes (e.g. Rodger et al., 2013). The international synoptic study by Berkey et al. (1974) used a global high-latitude riometer network in the Arctic region to construct and analyze the synoptic CNA maps during 60 substorms. Viewed at 15-min time resolution, these maps manifested some systematic dynamical changes of CNA distributions. In the absence of strong absorption during substorm growth phase, the impulsive CNA enhancement usually starts somewhere on the nightside (between 21 and 03 hr MLT) in the middle of auroral zone (roughly at 65-67°corrected geomagnetic latitude) at the substorm expansion onset, as confirmed also by Cresswell-Moorcock et al. (2013) superposed epoch study of the energetic electron precipitation which focused on latitudinal limits of precipitation. Later on the enhanced auroral absorption propagates eastward at the speed of a few kilometer per second, roughly consistent with the magnetic drift velocity of energetic electrons in the magnetosphere at around geostationary distance. After passing the dawn meridian, the intense slowly varying absorption bays are registered in the dawn-noon sector for more than an hour, forming a pronounced long-lived prenoon CNA maximum. The absorption magnitude decreases in the postnoon sector, being generally low throughout the substorm in the dusk sector. One may also observe some poleward expansion as well as occasional westward and equatorward expansion of the absorption region from its initial location. Some morphological elements like the nightside CNA maximum, pronounced late-morning maximum, and low intensity in the dusk sector can also be recognized in statistical EEP/CNA models like those published by van de Kamp et al. (2018), but important temporal dynamics is certainly missed in these statistical models.
Conjugate observations in the ionosphere and magnetosphere (most often at geostationary orbit) revealed that the CNA bays and associated X-ray bays observed by stratospheric balloons appear nearly simultaneously with the increases of energetic electron flux in the equatorial magnetosphere, with very similar intensity variations at magnetically conjugate locations (Baker et al., 1981;Parks, 1970;Spanswick et al., 2007). The multispacecraft studies at geostationary orbit established that initially the dispersionless-(inenergy) flux increase of energetic electrons (hereafter EE) appears in the near-midnight MLT sector, later on continuously expanding to the east (with larger time delays at later MLTs) with energy dispersion corresponding to the energy-dependent eastward magnetic drift of the electron cloud particles. The observation-based picture is consistent with the simulation studies (Birn et al., 1998;Li et al., 1998). Jointly they show that, as a result of reconnection-related flow bursts localized across the tail, the clouds of energetic electrons are injected into the nightside inner magnetosphere (see also Sergeev et al. (2012) and a recent summary by Gabrielse et al. (2019)), from where they drift around the Earth and gradually precipitate into the ionosphere. The injection is closely linked to the magnetic field dipolarization in terms of their amplitudes and locations, basically via the betatron acceleration (Gabrielse et al., 2019;Lopez et al., 1989;Nagai et al., 2019).
There were a few attempts to model the precipitation from the injected and drifting EE cloud. A comprehensive simulation was reported by Yu et al. (2018), who computed the drifts of different EE components in modeled magnetic and electric fields of the inner magnetosphere and included their pitch angle scattering to get the variations of D region ionization in the conjugate ionosphere. Particularly, after each substorm injection the enhanced fluxes of eastward propagating energetic electrons (10 to 100 keV) were found in the dawn sector outside the plasmapause, whose scattering by whistler mode chorus waves caused the ionization surges in the D region, propagating eastward from nightside to the noon in thẽ 5-10°wide auroral zone. However, the simulation was limited to the inner region at r < 6.6 R E and required the specification of injection source parameters at the outer simulation boundary from spacecraft observations at geosynchronous orbit (GEO). Therefore, this model could be suitable to simulate specific storm time substorms, but hardly helps to simulate moderate-and low-activity events, during which many injections are stopped outside GEO (Sergeev et al., 2012). Therefore, together with high computation cost, such comprehensive simulations are of limited help for monitoring purposes and space weather applications.
One more group of models is based on simple analytical calculations of electron drift in the dipole magnetic field with prescribed injection and precipitation rates to simulate different scenarios of precipitation dynamics and different shapes of CNA precipitation. From them, Liu et al. (2007), Liang et al. (2007), and Kabin et al. (2011) studies were focused on understanding the origins of the variability in possible riometer signatures, rather than on providing a framework for applications. Beharrell et al. (2015) developed a similar-type model describing the injection, drift, and subsequent precipitation of the energetic electrons, and attempted to calculate the substorm contributions to energetic electron precipitation into the atmosphere in actual events. For that purpose they evaluated some average characteristics of the "typical substorm injection" (neglecting differences of injection magnitudes in different events) and used the list of SuperMAG-based substorm onsets to initiate each particular event. The major problem (in this and other similar models discussed in above) is how to characterize quantitatively the substorm intensity and location of specific substorm injection and how to account for changing magnetic configuration. In addition, it is known that during substorms the dipolarization/injection region develops (expands or shifts) in several discrete steps rather than in continuous fashion (Sergeev, 1974), and that the location of initial activation, direction of expansion in the following steps, and other dynamical parameters are very individual in specific substorms. That means we need to use some monitoring tool which helps to specify the intensity and location of multiple injections in each event from the observational data.
Instead of describing analytically a such complicated chain of phenomena (including injection, electron drift in variable magnetic and electric fields, instabilities and precipitation via wave-particle interaction), in which many coefficients in the relationships should be somehow specified, in this paper we test another approach. Namely, following Shukhtina and Sergeev (1992), we attempt to describe complicated spatiotemporal behavior of precipitation by extracting the recurrent features of CNA substorm dynamics from the observational data. Moreover, our empirical approach specifically takes into account the variability and diversity of individual substorms.
One key ingredient of our approach is the characterization of injection intensity by the amplitude of magnetic field dipolarization. Indeed, computational models by Li et al. (1998), Zaharia et al. (2004), Liu et al. (2007), and Kabin et al. (2011), although being different in details, all showed quantitatively that magnetic field dipolarization of observed magnitude in the nightside inner magnetosphere is capable to provide the EE flux sufficient to explain the observed CNA response. As quantitatively demonstrated by Sergeev et al. (2011), Sergeev, Nikolaev, Kubyshkina, et al., (2014, a good proxy of dipolarization amplitude in the magnetosphere is provided by ground observations of the substorm current wedge (SCW; see Kepko et al., 2014) effects, observed as bay-like magnetic variations at midlatitude stations distributed around nightside. Using this approach we, in principle, can evaluate the SCW magnitude and location at every time step and, therefore, investigate the influence of variable spatiotemporal development of the dipolarization during the specific substorms. In this paper, we use the recently developed midlatitude positive bay (MPB) index (McPherron & Chu, 2017) as a possible SCW proxy.
A second key ingredient of our approach is the linear filtering technique (LPF; see, e.g., McPherron et al. (2015) for recent brief summary and references), which is used to extract the repeating features of substorm dynamics from the CNA records. First attempt to apply such approach for substorm-related CNA response reported by Shukhtina and Sergeev (1992) showed quite promising initial results, and here we considerably extend this model using advantage of better riometer coverage and much larger data sets available for both CNA and magnetic variations. In section 2 we describe observations and the event selection procedure and further test the amplitude relationship between the CNA bays and midlatitude magnetic bays. Then in section 3 we present the LPF method and compare its results with superposed epoch ones. Some validation tests are also presented there. The azimuthal dynamics of the auroral absorption and its interpretation in terms of the precipitation from drifting electron cloud are discussed in section 4. Finally, a brief summary of results and discussion of further steps to improve the performance of LPF method are presented in section 5.

Riometer Observations
In this paper we focus on azimuthal dynamics of the auroral absorption, that is on the modeling of dynamical CNA distribution observed by the riometers in the middle of auroral zone, at 65-67°corrected 10.1029/2019SW002385

Space Weather
SERGEEV ET AL. geomagnetic latitude. At these latitudes, roughly corresponding to the near-geostationary equatorial region, the CNA dynamics is closely related to the evolution of drifting EE clouds injected during the dipolarizations. Also, as follows from synoptic CNA maps presented by Berkey et al. (1974) and from statistical EE flux distributions in Figure 2 of van de Kamp et al. (2018), at these latitudes most of the time the energetic particle precipitation and CNA amplitude attain their maximal values. The extension of the dynamical model to include the latitudinal variations will be investigated in the following papers.
Five azimuthally separated riometers of the Canadian Space Agency's NORSTAR riometer array constitute the most valuable resource in this particular study. They include Gillam (GIL; 66.1°CGM latitude, magnetic midnight at 6.6 UT), Rabbit Lake (RAB; 67.0°, 7.5 UT), Fort Smith (FSM; 67.4°, 8.4 UT), Fort Simpson (FSI; 67.4°, 9.2 UT), and Dawson (DAW; 66.1°, 10.6 UT). CGM coordinates of the point at h = 100 km in the station zenith are calculated for the year 2010. University of Calgary Space Physics Data portal (http://data.phys. ucalgary.ca/) contains riometer data for roughly 30 years providing a possibility for massive statistical studies. These five stations are spread at roughly constant latitude over the 4-hr-wide MLT sector, allowing one to test the data quality, evaluate the inherent variability of CNA signal, and accumulate the large statistics. They are complemented by observations made at Russian Arctic station Tixie (TIX; 66.5°, 15.7 UT) as well as at Scandinavian stations Ivalo (IVA; 65.5°, 21.5 UT) and Abisko (ABI; 65.8°, 22.0 UT), allowing us to cross-check the results at widely displaced geographic locations. An example of observations during a sequence of several substorms is presented in Figure 1.
All these instruments are the standard wide-beam riometers operating at 30 MHz, with the riometer beam size being about 100 km (at 100-km altitude). The riometer observations are vulnerable to a number of problems including extraneous noise sources (both solar radio emissions and man-made ones). Also, the computation of CNA magnitudes relies upon derivation of quiet day curve baseline which may deviate from actual baseline due to various factors caused by changes in ionospheric conditions (changes in temperature gradients in the mesosphere) and technical problems (temperature conditions for the device, changes in the antenna-feeder system geometry due to icing and snow loads, interference caused by snowstorms, and domestic interference and radio noise). Therefore, the data quality control becomes an important issue in such types of studies. In the following we use CNA data from Canadian and Finnish riometers obtained using standard quiet day curve derivation procedure, which are available from the corresponding websites. In our case we use an advantage of closely spaced observations (within 1 hr MLT) made at nearly the same latitude, because the natural CNA events usually look rather similar in this case-compare, for example, records of GIL/FSM/DAW and IVA/ABI riometer pairs in Figure 1. General similarity of variations between neighboring riometer records helps to remove most of events significantly corrupted by noise interference and identify the serious baseline shifts by careful manual and automatic inspections of the data. Using the possibility to collect as many events in the database as required, we tried to clean the database by rejecting suspicious events of different kind. In addition, dealing with isolated substorms we usually benefit from having a quiet time period prior to the substorm onset in which the auroral absorption value is presumably low. A usual practice is that CNA magnitudes below 0.1-0.3 db are not treated as confident (e.g., Rodger et al., 2013), so the events with preonset absorption exceeding 0.4 db were deleted from the data set used to construct the prediction filter. The data obtained during the solar proton events (with enhanced fluxes of solar protons with energy >10 MeV exceeding 10 proton flux units according to GOES data) were also excluded. In total, due to all kinds of problems (including the absence of data) usually up to 30% of station events were lost from the long initial list of isolated substorms for the riometers used in this study. Figure 1 illustrates that every of five substorms seen in both AL and MPB data causes significant absorption increases in the auroral zone riometer records, which are seen in different MLT sectors. As known from previous studies (e.g., Berkey et al., 1974), the response usually starts on the nightside and propagates eastward to the later MLTs, its amplitude being usually lower in the dusk sector. The brief negative CNA spikes during daytime seen in 1-min traces correspond to solar radio bursts and are automatically filtered out when operating with 5-min median values. (The records with longer solar radio burst interference were removed from database.) This figure also illustrates some variability of CNA behavior, manifested by the differences of amplitudes and variation details observed at the nearby riometers.

Midlatitude MPB-Index: Event Selection Procedure
In order to quantify the power of the magnetic bay-like perturbations caused by the SCW, Chu et al. (2015) and McPherron and Chu (2017) produced a new MPB index. After removing both solar quiet and storm time variations, the sum of squares of the x and y residuals is determined at each of 35 midlatitude stations distributed in MLT. The power averaged over the nighttime sector provides the MPB index. This index was recently calculated at 1-min resolution for years 1980-2012. The multiyear list of the substorm onsets together with original MPB index data were provided in McPherron and Chu (2018), constituting the second main resource for our study, which requires many years of data to reach the goal.
The testing of MPB index for the substorm identification purpose was recently discussed by McPherron and Chu (2018). Because the onset identification is sensitive to a number of details (like threshold parameter values), two versions of substorm lists generated independently by X. Chu and by R. McPherron showed some difference in the numbers of identified events and, sometimes, in their onset times. Similarly, these results somewhat differ from the substorm onset lists generated based on auroral zone SML index in the SuperMAG project (Gjerloev, 2012;Newell & Gjerloev, 2011). Most difficulties actually arise during continuous activity and weak substorms, whereas less problems exist for isolated events. This is clearly seen in Figure 1, in which the MPB-based and AL-based onsets reasonably agree for isolated events #2 to #5. To be maximally safe concerning the onset time, which is a critical parameter in our study, we limit ourselves to isolated substorms. In addition, to minimize the influence of method details on the final result, we chose those events which appear in both MPB and SML onset lists and also ignored those events, in which the onset time differences between the lists were larger than 15 min. When selecting the events, we generally required that no substorm activity (T > 3 nT; see below) was observed at least 3 hr before the onset; we also avoided the situations with the long-lasting activity after the substorm.
In the following we suggest that intensities of dipolarization and injection are linearly related to each other. As the MPB-based proxy of the injection intensity we use a square root from MPB values (MPB 1/2 ). Their

Space Weather
SERGEEV ET AL. onset-to-peak differences ΔT will serve as the magnitude of the midlatitude magnetic bay. Weak events having ΔT <3 nT were discarded. Finally, the substorm onset time was taken as the first time point of bay-like variation in which the difference MPB 1/2 (t + 5 min) − MPB 1/2 (t) exceeded +0.5 nT 1/2 . In total, with severe restrictions described in this and previous subsections, only about 220 isolated substorm events were selected during years 2007-2012, with roughly 20 events per 1-hr UT bin available for the statistical study. Following the restrictions described in above, in Figure 1 only event #3 satisfied all criteria and was used in constructing the CNA response model.

Comparison of CNA Peak Amplitude to the MPB and AL Injection Proxies
Here we briefly discuss the observational justification of our choice how to drive the CNA model. Savinov et al. (1986) noticed that long events with continuing intense AL-index activity (like the steady convection events) may be accompanied by rather modest riometer absorption values, which are smaller compared to those recorded during substorm-related magnetic bays having similar AL magnitude. They also suggested to use the peak amplitude of CNA bay-like perturbations in the late morning sector as the estimators of injection strength and showed their fair correlation with the magnitudes of midlatitude magnetic bays. Here we follow the same approach. The choice of late morning sector is supported by the slowly varying character of CNA bays in this region (see, e.g., Figure 1), as expected at large distance from the injection region, where the injected electron clouds are interfered and smoothed due to the drift dispersion effects. This differs from a more impulsive appearance of CNA variations on the nightside, where the spikes and substantial differences of absorption variations at closely spaced stations are not uncommon (e.g., Spanswick et al., 2007).
In Figure 2a we show the relationship of CNA late-morning peak amplitudes (Ap) for each of three riometers RAB, FSM, and FSI, located at nearly the same magnetic latitude, with the total amplitude of midlatitude magnetic variation (ΔT). Quantities Ap and ΔT show a positive correlation with CC = 0.55 for 94 points included. Similar comparison with the peak magnitude of SML index (SuperMAG version of AL index), which potentially may also serve as a measure of substorm current wedge intensity during isolated substorms, shows a lower correlation (CC = 0.42) in Figure 2b. This illustrates that MPB-based midlatitude variation is a more perspective proxy of the dipolarization (and injection) magnitude than the AL-based product is. Relatively low values of correlation coefficients in Figures 2a and 2b are partly due to the low dynamical range of events in this data set, in which the medium-scale activity dominate, namely, events with SML~300-500 nT and ΔT~5-15 nT. The data set in Savinov et al. (1986)  In addition, Figure 2 indicates two other reasons of large scatter. The first is that, besides the dipolarization, there exists another factor which also influences the resulting energetic electron flux. Sergeev et al. (2015) showed that the temperature of ambient plasma sheet electrons (T e ) strongly depends on the solar wind state, so that statistically T e in the magnetosphere is considerably larger for fast and tenuous solar wind compared to the slow and dense wind. As the fluxes in the energetic tail of electron distribution should increase with the T e increase, we may expect the growth of >30-keV electron precipitated fluxes during fast tenuous solar wind conditions. The abovementioned differences are distinctly manifested in Figures 2a and 2b where systematically the red points (fast/tenuous wind) appear above the blue points (slow/dense wind). This is especially pronounced for Ap/ΔT correlation in which the linear regression slope for the red points is about 50% higher compared to the slope for the blue ones. This certainly introduces a significant scatter in the combined data set.
Another source of significant scatter is illustrated by Figure 2c, which compares the peak absorption values in the same event observed at the pair of closely spaced riometers and shows that they are rarely equal to each other. Whereas some influence of artificial factors (like the determination of actual quiet day curve) cannot be fully excluded, systematic appearance of these differences rather suggests their probable natural origin, that is, inhomogeneity of either the electron precipitation or the ionospheric/atmospheric parameters. The differences are rather large; for example, differences by a factor 2 or more are rather frequent even in the slowly varying CNA events typical for the late morning sector. As a result, the absorption correlation between two neighboring riometers in Figure 2c is only CC~0.7. All this indicates that the following reconstruction of response function characterizing the regular substorm-related dynamical features will be made in conditions 10.1029/2019SW002385

Space Weather
SERGEEV ET AL.
of considerable noise imposed on data by the external factors and by the inherent variability of energetic electron precipitation and/or ionospheric absorption.

Realization of LPF Method
Like other models of auroral absorption dynamics mentioned in section 1 Kabin et al., 2011;Liu et al., 2007;Liang et al., 2007;Yu et al., 2018), we believe that drift of (energy-dispersed) electron clouds injected in the nightside dipolarized regions provides the basic spatiotemporal context for the phenomenon considered, with the precipitation from this cloud being responsible for the observed CNA effects. This suggests a few implications of this conceptual approach. First, we believe that dipolarization proxy may serve as an input to describe the injected EE flux, with the betatron acceleration being behind this link. Second, different energy-dependent drift speeds imply the intrinsic time delays of CNA effects after the injection. Third, in case of several localized dipolarizations/injections contributing to the drifting cloud, we believe that linear summation/superposition provides an adequate mathematical approach to compute their integral effect.
Therefore, a choice of LPF model looks natural for our task, which has to include both the delayed response from particular injections as well as cumulative effect of multiple injections of different intensity. The LPF approach provides a simple way to model the response of a physical system to an external driver taking into account the delayed system response. Its simplest form gives the current output of the system as a weighted sum of the current and past inputs (see recent discussions and more details about the method in McPherron et al. (2015)).
Taking into account that the MPB rise time distribution during the substorm-related magnetic bays has a peak at 16 min (most frequently occurring between 5 and 30 min; see Chu et al., 2015, Figure 7), we selected the 5-min time interval as the basic time resolution of input and output data. Correspondingly, both CNA and MPB data sets are discretized to 5-min time resolution. Whereas we use the median CNA values in 5-min bin, for MPB product we use its differences over 5-min time step (T(t) = MPB 1/2 (t + 5 min) − MPB 1/2 (t)), as illustrated in the top part of Figure 3. Negative differences T are set equal to zero, because only positive differences were considered as the characteristics of dipolarizations/injections. Following Berkey et al. (1974) results we suggest that the relaxation time (when the CNA variation returns to its preonset level) is shorter than 4 hr. Formally, we describe that every injection T i causes a 5-hr CNA response, starting 1 hr before injection and ending 4 hr after, for example, Ai ¼ ∑ 60 k¼1 Fk × Tik , where F k is the desired response function (or the linear prediction filter). From causality a zero response is expected for the hour before injection, which helps to control the accuracy of the derived response function. For each isolated substorm the first CNA value A 1 is taken 1 hr before the substorm onset time t 0 so it requires knowledge of T values for 5 hr preceding t 0 (see T 1k in Figure 3). The last value in CNA data set (T 60k ) is taken 4 hr after t 0 . In total for every substorm at each station a 10-hr-long MPB data set (t 0 − 5 hr, t 0 + 5 hr) and a 5-hr-long CNA data set (t 0 − 1 hr, t 0 + 4 hr) are required. Adding another station or another event will increase the index i by 60, and so on. This is how the final matrix T and vector A are formed. Then the matrix equation was solved using singular value decomposition method (Forsyth et al., 1977). The function F has the sense of absorption response to the unit injection and in the following it is measured in db/10 nT units.
After that, using the known F, the output can be predicted from input data as Â = TF. Its correlation with measured A as well as the prediction efficiency index can be used to evaluate the quality of the prediction. Now we show the particular filters and investigate some methodical aspects using the collection of substorm events distributed in UT, which have been carefully identified during six-year time period in 2007-2012 according to criteria described in section 2. To increase the population of intense events we also added seven intense isolated substorms from the 1999-2007 interval. Since the response function should vary with the distance from the source injection, we group the events according to their MLT at substorm onset t 0 . We note that a so obtained LPF characterizes the precipitation response at the station rotating with the Earth rather than the response at the fixed MLT. Another thing to notice is that, for the given MLT bin, the events at different stations start at different UTs, so the combination of data from different riometers helps enormously to increase the statistics of LPF derivation without accumulation and processing many additional substorm events.
In view of observed CNA differences at closely spaced stations (Figure 2c), our first question is whether the response functions derived at the same MLT (and magnetic latitude) but at different stations are similar and,

Space Weather
SERGEEV ET AL.
therefore, provide an objective information about the substorm effects in the electron precipitation. In Figures 4a and 4b we compare the LPFs obtained separately using data of GIL, RAB, FSM, and FSI riometers. Here the amount of points in the input data sets varies between 10 and 22, and we emphasize that data sets are mostly composed from different events, the particular event lists being slightly overlapped for stations closest in longitude. In spite of the small differences between filters, generally their shapes and amplitudes look very similar. The correlation coefficients between predicted and observed time series are between 0.6 and 0.7, and PE varies between 0.4 and 0.6. Similarity of response functions obtained for different stations suggests that statistical shape of the prediction filter provides an objective information about the precipitation response to substorms. We note that the agreement between different stations is not always as good as in Figure 4, especially for intervals with relatively small amount of events (10)(11)(12)(13)(14)(15)(16)(17)(18)(19)(20) in the data set, such differences being most probably caused by the specific (usually strong) particular events under insufficient statistics.
Encouraged by the similarity of LPFs obtained for individual stations, we used the combined data set of four Canadian riometers (GIL, RAB, FSM, and FSI) to investigate the global development of auroral absorption in the middle of auroral zone. Figure 5 shows the response functions in different MLT sectors. (More details are presented in Figure S1 and Table S1 in the supporting information.) Each LPF was constructed for all events in which some station was located in 2-hr-wide MLT bin at the substorm onset time t 0 , so now the amount of station events varies between 36 and 81 in different bins. The variations of median and quartile CNA values obtained by the superposed epoch analyses (SEA) are also shown for comparison. Generally, both methods agree in describing the basic trends known from previous studies, including initial impulsive CNA increase on the nightside from where it propagates eastward and the increase of both rise time and duration while shifting from midnight to late morning. At the same time there are also some noticeable differences between LPF and SEA results which will be discussed later in section 4.

Validation of LPF Method
According to Figure 5 for MLT bins where the precipitation amplitude is large (between midnight and noon inclusively), the correlation is at the level CC~0.6 … 0.7, PE is positive and varies between 0.35 and 0.5. Naturally, the low variance of CNA on the duskside predetermined very low CC and PE values in this sector.
To characterize the overall quality of reconstruction we carried out two kinds of tests. First, we compared the predicted and observed values for all events in the sector 06-12 hr MLT. The predicted CNA values were computed based on MPB input and on the corresponding LPF filters demonstrated in Figure 5, and then the summary scores were calculated for the entire data set in the late morning sector. For independent The results detailed in Table S2 and Figure S2 show that peak calculated absorption values demonstrate good correlation with the observed ones (CC~0.7); however, the PE values are negative, possibly indicating the shifts in the absolute values between predictions and observations. Particularly, despite the small distance between Abisko and Ivalo (separated by 143 km) their peak values differ almost 1.5 times under fair correlation between them (CC = 0.9) which needs further investigation. The results for the total data set of 5-min values are similar, although the quality of prediction is slightly worser (CC~0.63).
We also tested whether our LPF functions derived for isolated substorms can be applied to predict the auroral absorption during disturbed period consisting of nonisolated events. For that purpose we explored the disturbed five-day period between 27 April and 1 May 2007, containing 61 substorms, which was previously considered by Beharrell et al. (2015). The results for two Canadian stations, demonstrated in Figure 6, show that the prediction quality, as given by CC and PE values, is similar to that obtained for the collection of isolated individual substorms studied in our paper. A closer inspection shows that the observed peak values (red traces) in most cases systematically exceed the predicted values (blue traces). If confirmed in a more detail study, such difference may be related to the effect of substorm clusters which, according to Rodger et al. (2016), can produce stronger EE fluxes and wave intensities than the isolated substorms can. Whereas this

Space Weather
interesting aspect needs a special study, the comparable results for isolated/nonisolated substorms encourage us to extend the application of LPF method for any kinds of substorm activity.

Global Dynamical Properties of Electron Precipitation as Revealed by the LPF Analyses
Twenty-four LPF response functions (obtained in 2-hr-wide MLT bins centered at 00, 01, … 23 hr MLT; all original LPFs are presented in the Figure S1) have been calculated using the joint data set of four Canadian riometers. They are used to generate the map of response function dynamics in the substorm time-MLT coordinates-see Figure 7 (left). In similar way, Figure 7 (right) illustrates the azimuthal dynamical pattern of the median absorption values. These two panels provide a view of azimuthal substorm dynamics of the energetic electron precipitation (auroral absorption).
Whereas two patterns are generally similar, there are also some important differences between LPF and SEA results, which are obvious from comparison of two plots and from inspection of original LPFs obtained at 1-hr resolution ( Figure S1). The first one is that both the rise time and duration are shorter for LPFs compared to the simple averaging (SEA), showing the differences between the effect of individual injection (LPF) and the cumulative effect of multiple injections, as provided by the SEA curves.
The second big difference is that, compared to SEA, the LPF results allow to resolve the appearance of a strong and short (about 10 min long) impulsive peak in the presumed injection region, between 23 and 02 hr MLT (see also Figure S1). This peak is completely smeared out in SEA median results. In agreement with previous studies  we  interpret it as a manifestation of precipitation from the newly injected electron cloud. The valley between the nightside peak and more smooth late-morning peak centered at around 08 hr MLT is not seen in the median SEA curves, although it is still can be recognized in the upper quartile results shown in Figure 5. Third, the late morning maximum looks like a part of gradual large-scale evolution of precipitation from the eastward drifting electron cloud. It includes a gradual LPF increase with MLT from its minimal values at around 02-03 hr MLT to the maximal values at around 08 hr MLT, as well as the gradual LPF peak decrease from 08 to 16 hr MLT. As previously argued by Shukhtina and Sergeev (1996), these gradual changes in the energetic particle fluxes can mostly be contributed by the electron acceleration/deceleration in the large-scale convection electric field, although these changes of precipitation can also be mediated by the corresponding variations of the associated chorus wave activity.
The eastward shift of the precipitation region looks as a very ordered feature in Figure 7. To evaluate its velocity and associated changes of the response function, in the bottom panel of Figure 8 we replot the prediction filters obtained in the midnight-dawn-noon sector. For each filter curve we also identify the substorm time and MLT coordinates of its half-width points (based on the 0.5 of peak value in each MLT bin), which are presented in the top panel. Two rear-side points are shown for the filter at 11 hr MLT because the half-value level crosses twice the LPF curve (at 01:45 and 02:35 of substorm time), which has a kind of plateau in that region. Also, two rear-side points are shown for the first filter depending on whether we use sharp peak (only one point) or the main peak in the analyses.
These figures illustrate a couple of interesting features. First, although visually the filters and SEA curves in Figure 5 look wider with the increasing local time as expected from the energy-dependent drift dispersion of the electron cloud, the formally identified half-width in Figure 8 (and green LPF isocontour line in Figure 7 a) does not support this expectation. Actually, one may notice that, immediately after the peak, the CNA decreases fast but then the decrease becomes slower for smaller values of CNA (see the curves for 05, 07, and 09 MLT). This change of the relaxation time may be related to different processes. For example, it may occur due to the flux change associated with either the atmospheric losses of precipitated particles or with the intensity variation in the drift-dispersing cloud. Alternatively, this difference can be due to the specifics of the precipitation mechanism. Particularly, for the cyclotron instability the electron loss time is shorter for stronger filling of the loss cone which increases with the increase of available particle flux (under equal other conditions). In reality the diffusion rate has a nonlinear relationship with the trapped electron flux, in part because chorus waves are strongly enhanced by substorms (Meredith et al., 2001). 10.1029/2019SW002385

Space Weather
As concerns the observed azimuthal propagation velocity, which is about 12 hr MLT in roughly 1 hr (Figure 8, top), we notice that this is equal to the drift period of 50-keV electron on the dipole drift shell L = 7.5. The choice of this L value is based on the equatorial B z value (73 nT at L = 7.5), as well as on the results of , who found that the inner boundary of the energetic particle injection to geostationary orbit statistically occurs for conditions when B z~7 0 nT (see Figure 12 in that paper). The electron energy 50 keV corresponds to the results by Kellerman et al. (2015) study, who used two different methods (correlation between CNA magnitudes and trapped energetic electron fluxes, as well as study of the drift delays at the riometer network) to evaluate the most probable energy of energetic electrons contributing to the auroral absorption. The observed azimuthal propagation velocity is, therefore, consistent with the predictions of electron drift cloud model.

Discussion and Concluding Remarks
Here we demonstrate a new empirical approach to reconstruct the substorm dynamics of auroral absorption which potentially may help in future operational predictions of substorm-related D region ionization enhancements in the auroral zone. It is based on linear prediction filter technique, which (in this version) exploits a large database of riometer observations in the auroral zone combined with midlatitude MPB index designed to quantify the substorm dipolarization effects. Our approach is based on our current physical understanding of (a) close relationship between magnetic field dipolarization in the near magnetotail with the flux increase (injection) of energetic electrons and (b) drift dynamics of injected EEs whose precipitation into the ionosphere provides the observed auroral absorption. Correspondingly, we expect that response functions (LPF) of auroral absorption reconstructed at different MLTs are normalized to the dipolarization magnitude and, taken together, provide the quantitative global dynamical image of the precipitating electron clouds drifting after injection, as shown in Figures 5, 7, and 8. Two abovementioned foundations of the method receive further support in the results of our study. The correlation between peak absorption and amplitudes of MPB (and SML) indexes shown in Figure 2 confirms the close relationship between the dipolarization and injection amplitudes, in addition to other facts summarized by Gabrielse et al. (2019) and Nagai et al. (2019). On the other hand, very regular eastward propagation of the CNA pattern at the velocities of electron drift in the magnetosphere (see previous section), manifested in Figures 7 and 8, supports the concept of precipitation from a drifting electron cloud. Also, the general pattern of enhanced (or depressed) precipitated fluxes at dawn (or dusk) is consistent with acceleration (or deceleration) of energetic electron fluxes by the dawn-dusk directed convection electric field as previously argued by Shukhtina and Sergeev (1996).
Compared to previous studies of global CNA dynamics the LPF approach qualitatively reproduces all its basic features derived so far from global synoptic studies of particular events and from statistical riometer studies (Berkey et al., 1974;Kellerman et al., 2015;Kellerman & Makarevich, 2011, etc.). According to them, initially the CNA starts impulsively on the nightside (Kellerman & Makarevich, 2011). Then it propagates eastward, consistent with the EE magnetic drift (Kellerman et al., 2015). The peak CNA values are achieved, first, in the injection region and, second, in the broad late-morning maximum centered at 08 MLT which occurs roughly 0.5-1 hr after the injection (Berkey et al., 1974, etc.). This is where the EEs gain maximal energy and flux drifting in the dawn-dusk convection electric field, which stimulates the excitation of strongest whistler waves scattering electrons into the atmospheric loss cone (Meredith et al., 2001). The CNA amplitude significantly decreases in the afternoon sector, fading to low values in the dusk quadrant. Statistical pattern of CNA occurrence (e.g., Hargreaves, 1969) has two absorption maxima near midnight and in the late morning sector, in close agreement with our LPF results shown in Figure 7. Being consistent with previous studies, our results, we believe, provide an improved and refined quantitative global image of the precipitation dynamics.
Whereas this paper demonstrates the potential of the LPF approach for investigation and possible predictions of substorm-controlled dynamics of the EE precipitation, this is certainly the first step. There is a number of issues to be investigated this way. On the one hand, the modeled and observed CNA variations are similar but the prediction efficiency is not that high, CC = 0.64 and PE < 0 in the validation effort described in section 3.3. Here we note that the prediction accuracy has natural limits. For example, from narrow-beam imaging riometer studies it is known (see Makarevitch et al., 2004, and references therein) that intense absorption in the late morning sector consists of patchy structures at tens of kilometer scale which show a large variability in the patch shapes and velocities. It may originate from cold plasma structures influencing the chorus generation process or/and from the limited scale size of instability regions, which are much smaller than riometer beam diagram size as well as distance between riometers in our case. This is our explanation of the large intrinsic variability of riometer absorption demonstrated by comparison of peak values seen at closely spaced riometers in Figure 2c. We remind that correlation of CNA peak values simultaneously observed at closely spaced riometers in our data set was only CC~0.7 in Figure 2c, which clearly puts a limit for the accuracy of possible predictions.
On the other hand, an obvious resource exists to reduce the scatter and improve the predictability of the LPF method by including the solar wind control of ambient electron parameters in the conjugate magnetosphere, which has been demonstrated in Figures 2a and 2b. Previously significant influence of solar wind velocity on the CNA magnitude was statistically demonstrated by Kavanagh et al. (2004). Since the temperature of ambient plasma sheet electrons (T e ) strongly depends on the solar wind state with T e being considerably larger for fast and tenuous solar wind compared to the slow and dense wind (Sergeev et al., 2015), the amount of seed electrons which can be accelerated to >30-keV energy behaves similarly to T e . The exact way how to describe this influence (choice of the optimal function of solar wind parameters and of required time delays, etc.) are the subject of special study which is underway.
Another possibility to improve the method performance concerns the proxy parameters used to describe the injections. Due to the large database available, the MPB index is an excellent resource for various studies of the LPF method; however, with the MPB index we only take into account the variable intensities of multiple dipolarizations/injections but not their variable geometrical parameters. Meanwhile, according to modeling results by Liu et al. (2007), Liang et al. (2007), etc., the width, location, and their changes with the substorm time are also important factors when considering the evolution of particle flux and energy spectra. The magnetogram inversion methods (e.g., Chu et al., 2015;Sergeev et al., 2011) can provide such information, although their usage currently requires significant processing efforts and is still far from being fully automatized. Still this continues to be an important research area.
Last but not least is that in this study we only addressed the longitudinal dynamics of EE precipitation. The inclusion of latitudinal dynamics, which takes into account both the substorm expansion features and the variable size of auroral oval (which is especially important during strong activities, like the magnetic storms), will certainly be a next important, though nontrivial, big step in the adaptation of the proposed methodology for the application purposes.