The dependence of magnetospheric plasma mass loading on geomagnetic activity using Cluster

Understanding changes in the magnetospheric mass density during disturbed geomagnetic conditions provides valuable insight into the dynamics and structure of the environment. The mass density plays a significant role in a variety of magnetospheric processes, such as wave propagation, magnetic reconnection rates, and radiation belt dynamics. In this study, the spatial variations of total plasma mass density are explored through the analysis of Cluster observations. Data from the WHISPER (Waves of High frequency and Sounder for Probing of Electron density by Relaxation) and CODIF (ion Composition and Distribution Function analyzer) instruments, on board the four Cluster spacecraft for a time interval spanning 2001–2012, are used to determine empirical models describing the distribution of the total plasma mass density along closed geomagnetic field lines. The region considered covers field lines within 5.9≤L < 9.5, corresponding to the outer plasmasphere, plasmatrough, and near‐Earth plasma sheet. This study extends previous work to examine and quantify spatial variations in the electron density, average ion mass, and total plasma mass density with Dst index. The results indicate that during periods of enhanced ring current strength, electron density is observed to decrease and average ion mass is observed to increase, compared with quiet geomagnetic conditions. The combination of these variations shows that although heavy ion concentration is enhanced, the decrease in plasma number density results in a general decrease in total plasma mass density during disturbed geomagnetic conditions. The observed decrease in mass density is in contrast to prevailing understanding and, due to the dependence of the Alfvén speed on mass density, has important implications for a range of plasma processes during storm time conditions (e.g., propagation of wave modes).


Introduction
The magnetospheric mass density, and its variations under different conditions, is a fundamental parameter of the magnetosphere. The mass density is a significant factor for the propagation of wave modes concerned with radiation belt energization and decay [Li et al., 1993;Elkington et al., 1999;Meredith et al., 2003;. Mass density also has implications for the coupling of the solar wind to the magnetosphere, through its role in determining dayside reconnection rates [Borovsky and Denton, 2006]. It is known that the magnetosphere experiences significant distortions and global disturbances during geomagnetic storms, including large enhancements of the ring current [Akasofu and Chapman, 1961;Akasofu et al., 1963;Gonzalez et al., 1994]. Therefore, an understanding of how the mass density varies with geomagnetic activity is of importance and will provide useful information on the dynamical processes and morphology of the magnetosphere during disturbed conditions. A key motivation of this study is the application of a mass density model to understand the response time of the magnetosphere to both internal and external perturbations, through variations in the frequencies of magnetospheric ultralow-frequency waves.
This study is concerned with examining mass density changes during geomagnetically disturbed periods, by determining an empirical model of the spatial distribution of total mass density in the closed magnetosphere (specifically the outer plasmasphere, plasmatrough, and near-Earth plasma sheet) that considers variations with geomagnetic activity. Previous models describing the dependence of mass density with geomagnetic activity exist [Menk et al., 1999;Gallagher et al., 2000;Takahashi et al., 2002;Berube et al., 2005;Denton et al., 2006;Takahashi et al., 2006;Maeda et al., 2009;Takahashi et al., 2010], although they are not without

10.1002/2017JA024171
Key Points: • Empirical models describing how spatial variations of electron density, average ion mass, and total mass density vary with D st index • For disturbed conditions, number density is decreased and average ion mass is increased, which combines to result in decreased mass density • Statistical study of WHISPER and CODIF data from the Cluster spacecraft between 2001 and 2012 for 4.5 ≤ L < 9.5 Supporting Information: • Supporting Information S1 • Figure S1 • Figure S2 • Figure S3 Journal of Geophysical Research: Space Physics 10.1002/2017JA024171 their limitations. Most of these studies indirectly infer the mass density on a field line from observations of harmonic standing Alfvén wave frequencies, either from spacecraft measurements [Takahashi et al., 2002;Huang et al., 2004;Denton et al., 2006;Takahashi et al., 2006;Tu et al., 2006;Takahashi et al., 2010] or ground-based magnetometer arrays [Menk et al., 1999;Berube et al., 2005;Maeda et al., 2009]. This method requires the assumption of a functional form describing the field-aligned distribution of mass density. Many of the mentioned studies assumed a power law form to describe the field-aligned variations [Menk et al., 1999;Takahashi et al., 2002;Huang et al., 2004;Berube et al., 2005;Takahashi et al., 2006;Tu et al., 2006;Maeda et al., 2009;Takahashi et al., 2010]; however, previous observations indicate that this form does not account for a local peak in mass density at the magnetic equator [Takahashi et al., 2004;Denton et al., 2006;Takahashi and Denton, 2007;Denton et al., 2009;Sandhu et al., 2016a]. Therefore, this choice of functional form may have introduced some inaccuracies to the results. In contrast, a study conducted by Denton et al. [2006] accounted for the peaked field-aligned distribution of mass density, although the CRRES (Combined Release and Radiation Effects Satellite) data set used in the study provide poor coverage over magnetic local time (MLT). The spatial coverage is such that only 11% of the data lie outside of the 12-18 MLT sector.
All the studies mentioned previously used an indirect inversion technique to infer the mass density of magnetospheric plasma. The Global Core Plasma Model [Gallagher et al., 2000], a combination of previous models, provides the spatial distribution of density for the key magnetospheric regions, including variations with geomagnetic activity. This model is based on in situ observations of mass density from the Retarding Ion Mass Spectrometer instrument on board DE 1 [Gallagher et al., 1988] and electron density observations from measurements of plasma characteristic frequencies by the Sweep Frequency Receiver instrument on board ISEE 1 [Carpenter and Anderson, 1992]. Although the model provides good spatial coverage, the key limitation is that the model mainly focuses on electron density distributions, with restricted detail on the contribution of ion composition to the total mass density.
This study is motivated by the lack of models providing a sufficiently accurate description of magnetospheric mass density distribution and the variations with geomagnetic activity while maintaining good spatial and temporal coverage. The data sets used in this study are relatively large and cover a sufficiently long time interval, as well as using direct measurements of the plasma. Therefore, the resulting empirical model provides a statistically reliable description of the distribution of mass density on closed field lines and how the distribution varies with geomagnetic activity.
In this study, we develop a previously determined empirical model of the total mass density in the closed magnetosphere [Sandhu et al., 2016a]. The previous study will now be briefly discussed, highlighting the key features of the model and reviewing the approach used. Full details are presented by Sandhu et al. [2016aSandhu et al. [ , 2016b. By using observations of the electron density and the average ion mass obtained by Cluster, specifically the WHISPER (Waves of High frequency and Sounder for Probing of Electron density by Relaxation) and CODIF (ion Composition and Distribution Function analyzer) instruments, the field-aligned distributions of these values were modeled. These empirical models of the electron density and average ion mass represented the spatial distribution of plasma for the regions between 4.5 ≤ L < 9.5 and 5.9 ≤ L < 9.5, respectively. The L shell ranges correspond to the outer plasmasphere, plasmatrough, and near-Earth plasma sheet, and the models provided coverage over all MLT (magnetic local time) sectors. The models described how the values vary along magnetic field lines and included dependences on L shell and MLT.
The field-aligned distribution of electron density was observed to follow two separate forms, depending on the region considered. At the high-latitude region of the field lines, the electron density was modeled using a power law form function with a positive power law index [Cummings et al., 1969], where the electron density increases toward the ionospheric ends of the field lines and approaches a minimum toward the magnetic equatorial plane. Although this dependence is commonly used to represent the distribution of electron density along the whole field line [Décréau et al., 1986;Olsen et al., 1987;Olsen, 1992;Gallagher et al., 2000;Goldstein et al., 2001;Denton et al., 2002Denton et al., , 2004, the results from Sandhu et al. [2016a] indicated that the power law form is inappropriate to represent the plasma at low latitudes, close to the magnetic equatorial plane. An equatorial enhancement, that is increasingly prominent for increased L values and toward the nightside, is observed and is suggested to be a result of solar wind-magnetosphere coupling acting as an additional source of plasma in the closed magnetosphere.
In addition to the electron density, the field-aligned distribution of average ion mass was also modeled. The data showed that the average ion mass followed a power law form with a negative power law index. This form 10.1002/2017JA024171 describes the average ion mass approaching a maximum toward the magnetic equatorial point of a field line, with values decreasing off equator. This represents the effects of the centrifugal force, acting more effectively on the heavier ions and preferentially trapping them at the magnetic equator [Lemaire and Gringauz, 1998;Takahashi et al., 2004;Denton et al., 2006Denton et al., , 2009. A key feature of the average ion mass model is the enhancement at low L values in the nightside sector. This is expected to result from the mass dispersion of outflowing ionospheric ions as they convect into the plasma sheet and inner magnetosphere [Sandhu et al., 2016a, and references therein].
By combining the spatial distributions of electron density, n e , and average ion mass, m av0 , the corresponding spatial distribution of the total plasma mass density, , was inferred from This made the reasonable assumptions that the plasma is quasi-neutral and the electron mass is negligible compared to the ion mass. Furthermore, the approach assumes that the average ion mass measurements provided by the CODIF instrument are representative of the full plasma population and are not a strong function of ion energy for the region considered in this study. This assumption has been validated by Sandhu et al. [2016b], specifically referring the reader to Figures 2 and 3 of Sandhu et al. [2016b] for further details. From equation (1), it can be seen that the model included contributions from both the number density and ion composition, in determining total mass density. The present study aims to develop this model further, by including dependences on geomagnetic activity, providing information on how the spatial distribution of the total mass density in this region varies during geomagnetically disturbed intervals. This will involve using a similar statistical approach for the same data sets from the WHISPER and CODIF instruments on board the Cluster spacecraft. Empirical models describing the field-aligned distribution of electron density and average ion mass will be determined, including variations with L, MLT, and geomagnetic activity, which will then be combined to obtain the desired mass density model.

Method
This study uses measurements obtained by the Cluster spacecraft, which consists of four identical satellites (C1, C2, C3, and C4). For further details on the Cluster mission, the reader is referred to Escoubet et al. [1997] and Escoubet et al. [2001]. The electron density and average ion mass data sets required are obtained from the WHISPER and CODIF instruments respectively, provided by the Cluster Science Archive (CSA) for all four of the spacecraft, and the full data sets cover the interval from 2001 to 2012. The WHISPER instrument, a resonance sounder, measures the total electron density of the local plasma within the range 0.25-80 cm −3 in the active mode [Décréau et al., 1997;Trotignon et al., 2001Trotignon et al., , 2003. The CODIF instrument, part of the CIS (Cluster Ion Spectrometry) experiment [Rème et al., 1997], provides density measurements of the key ion species (H + , O + , and He + ) covering an energy range of 0.02-40 keV/charge, which are used to estimate the average ion mass. For the average ion mass data set, measurements from the Cluster spacecraft C2 (where the CIS instrument is nonoperational) and C3 (which has instrumental issues with the CODIF sensor) were excluded, although the coverage provided remains sufficient for reliable results. Furthermore, CODIF observations corresponding to spacecraft L values below 5.9 were rejected from the data set, due to strong background contamination by penetrating energetic radiation belt particles at lower L values. Additional details on the motivation for choice of instruments and methodology used are included in Sandhu et al. [2016aSandhu et al. [ , 2016b.

Data Reduction
The data reduction applied to the WHISPER and CODIF data sets, presented by Sandhu et al. [2016aSandhu et al. [ , 2016b, will now be briefly reviewed. The data were initially binned for position in the GSM (geocentric solar magnetospheric) coordinate system, using a bin size of 0.5 × 0.5 × 0.5 R E and taking each orbit separately. For observations in each position bin for each orbit, the average measurement value and average time of measurement were determined, from which the corresponding MLT (magnetic local time) was also found.
For each position bin, the corresponding L value was found, defining the L value as the radial distance of the field line position in the magnetic equatorial plane, according to the T96 magnetic field model [Tsyganenko, 1996]. The T96 magnetic field model is parametrized by the solar wind dynamic pressure, IMF (interplanetary magnetic field) B y and B z components, and the D st index. The parameter values corresponding to the average measurement time of each bin were obtained from the NASA/Goddard Space Flight Center OMNI data set through OMNIWeb, for the 1 min averaged solar wind parameters, and from the World Data Center for Geomagnetism (Kyoto) data set, for hourly averaged D st values.

10.1002/2017JA024171
As this study is concerned only with the closed magnetosphere, any field lines which the T96 magnetic field model traced as open were discarded from the data set. However, it was identified that some open field lines were incorrectly traced as closed, particularly close to the magnetopause. These measurements were identified using a method adapted from Clausen et al. [2009], which involved examining the HIA/CIS ion energy spectrogram data at the time corresponding to the average time of measurement of the position bin. By identifying whether the Differential Energy Flux of the ions in the region was characteristic of open or closed field lines, the position bins corresponding to open field lines were removed from the data sets. A detailed discussion of this process is included in the previous study.
By using these data reduction techniques, the electron density and average ion mass data sets, restricted to observations within the closed magnetosphere and binned for position, with corresponding MLT and L values determined, were obtained. These data sets provide spatial coverage along field lines. The WHISPER and CODIF data sets provide spatial coverage along field lines in the region spanning 4.5 ≤ L < 9.5 and 5.9 ≤ L < 9.5, respectively, representing the outer plasmasphere, plasmatrough, and near-Earth plasma sheet. The different lower limit in data coverage for the CODIF data set, compared to the WHISPER data set, arises due to the removal of data contaminated by background penetrating radiation, as previously mentioned.

Field-Aligned Distribution
The electron density and average ion mass data sets are used to examine the distribution along magnetic field lines. Whereas the previous studies [Sandhu et al., 2016b[Sandhu et al., , 2016a] compared data at different MLT and L values, this study also compares data at different levels of geomagnetic activity. To represent the level of geomagnetic activity, the hourly averaged D st index is chosen, and the D st index value is determined for each position bin. It is known that geomagnetic storms are characterized by a significant reduction in the horizontal component of the global magnetic field relative to the average level [Chapman, 1918], as a result of an enhanced ring current during the main phase [Dessler and Parker, 1959;Sckopke, 1966]. The D st index represents the magnitude of these global variations of the horizontal field close to the magnetic equator. Therefore, D st index is an appropriate proxy for the level of geomagnetic activity in the magnetosphere.
The variations along magnetic field lines are modeled for both the WHISPER electron density data set and the CODIF average ion mass data set, where details on the modeling process are included in Appendix A. In order to include dependences on the D st index, the data sets are binned for D st index. Figure 1 shows the frequency, representing the number of measurements, for a range of D st index values. The majority of data lies between −100 to 10 nT, where this range corresponds to conditions varying from moderately disturbed to typical quiet times. Therefore, only data within this range were selected, avoiding extreme values distorting the variations. The data are then divided into six bins for D st index, as indicated in Figure 1 by the different colors. It is noted that the bin widths used vary with D st index, such that bins at lower D st index values have an increased bin width. This was done to ensure that each bin contained a sufficient number of data points for the fits to remain statistically reliable, as it is apparent in Figure 1 that the frequency is significantly decreased at largely negative D st index values. Figure 1 shows that each bin contains at least 170 data points, providing good data coverage.

Electron Density
The electron density model, describing variations along magnetic field lines and including dependences on L, MLT, and D st index, is given below in equation (2), where details on the fitting are included in Appendix A. The position along the magnetic field line is represented by the normalized radial distance, R norm , which is the radial distance at which the measurement was obtained, R (R E ), divided by the L value of the field line. The empirical electron density model describes variations in the L range of 4.5 ≤ L < 9.5 and covers all MLT. In order to provide context in terms of the statistical variations, the standard error of the electron density values for a given L-MLT-D st index bin is 1.2 cm −3 . The spatial distribution described by the resulting model is shown in Figures 2a-2d, where values are mapped azimuthally to the X-Z plane and averaged separately over dayside and nightside MLT sectors. The SM (solar magnetic) coordinate system is used here, where the geomagnetic dipole axis is aligned with the Z axis, in order to clearly illustrate variations with L value and magnetic latitude. The azimuthal mapping is done by assigning the radial distance of the bin from the Z SM axis to the X SM coordinate. The panels show the distribution at four values of D st index, chosen to indicate variations between disturbed and quiet conditions. In addition, Figures 3a-3d show the distribution of electron density predicted by the model (equation (2)) in the T96 magnetic equatorial plane. Each panel corresponds to a different value of the D st index, using the same D st index values as in Figure 2. Figures 2 and 3 demonstrate some clear dependences of the electron density spatial distribution on D st index, which are examined further in section 5.

Average Ion Mass
The average ion mass distribution along field lines is also modeled, including dependences on L, MLT, and D st index; refer to Appendix A for further details. The resulting model, shown below in equation (3), covers field lines over all MLT and spanning 5.9 ≤ L < 9.5. The corresponding standard error for the average ion mass values is 0.40 amu for a given L-MLT-D st index bin. Further details on the standard error of the electron density and average ion mass values are included in the supporting information (S1).

Mass Density Model
The empirical models for the field line distribution of electron density (equation (2) discussed in section 3.1) and average ion mass (equation (3) discussed in section 3.2) are combined, using equation (1), to infer the corresponding model for the total plasma mass density including dependences on D st index. Due to the spatial coverage of the WHISPER and CODIF data sets, this mass density model represents plasma in the region covered by 5.9 ≤ L < 9.5.
This resulting model is used to examine variations in the X-Z SM plane for a range of Although the mass density model presented here is parameterized by D st index, many other studies use K p index as a proxy for the level of geomagnetic activity [Menk et al., 1999;Gallagher et al., 2000;Takahashi et al., 2002;Denton et al., 2006;Takahashi et al., 2006Takahashi et al., , 2010Takahashi et al., , 2014. Similarly to the D st index-dependent models, models describing the field-aligned distributions of electron density and average ion mass were obtained, except that variations with K p index were included as opposed to variations with D st index. Further details on this analysis are included in Appendix B. The dependences on K p index will be compared to the corresponding D st -dependent models, providing information on how using a different proxy for the level of geomagnetic activity affects the spatial distributions. This discussion is presented in section 5.

Discussion
This section will now discuss the key features of the empirical D st -dependent models for electron density and average ion mass and the inferred mass density model. As the previous study examined the dependences of the electron density, average ion mass, and total mass density on L value and MLT in detail [Sandhu et al., 2016b[Sandhu et al., , 2016a, this discussion will predominantly focus on the dependences of the field line distributions on D st index. Furthermore, the effect of quantifying variations with K p index, as opposed to D st index, is also considered. This will provide an understanding of how the spatial distributions vary with geomagnetic activity, due to changes in the mass loading processes in the closed regions of the outer plasmasphere, plasmatrough, and near-Earth plasmasheet.

Electron Density
Using the same functional form as the previous study, the field line distribution of electron density at high latitudes is described using a power law function (equation (2a)). This form gives a field line distribution where the electron density is a minimum at the magnetic equatorial point and increases toward the ionospheric ends of the field line (where the power law index is constrained to be a positive value). The power law form represents the contribution of the ionosphere as a plasma source, loading the flux tubes through photoionization by incident solar radiation. This results in a distribution where the electron density along a field line approaches a maximum toward the ionosphere and decreases moving away from this plasma source. The model parameters for the power law function (equation (2a)) are the power law index, , (equation (2d)) and the electron density at the magnetic equatorial point of the field line, n e0 (equation (2c)), which both contain information on how the power law field-aligned distribution varies with L, MLT, and D st index.
The power law index, , (equation (2d)) describes the gradient of the electron density decrease away from the ionosphere. As previously identified by Sandhu et al. [2016a], the power law index is observed to linearly increase in value with L value, which is attributed to the increase in flux tube volume and length, resulting in a decreased total electron density across the flux tube. Therefore, the decrease in electron density moving away from the ionosphere will be greater, corresponding to an increased value of .
A dependence of the power law index, , on MLT is also observed, indicated by the sinusoidal component in equation (2d), in agreement with the average model [Sandhu et al., 2016a] and previous studies . The MLT dependence is such that is peaked at approximately dawn. As discussed in Sandhu et al. [2016a], the MLT asymmetry is due to increased incident flux of solar radiation for dayside field lines, resulting in increased photoionization at the flux tube footprints.
The power law index is not observed to demonstrate any statistically significant or coherent dependences on D st index. However, results from previous studies suggest that refilling rates from the ionospheric plasma source vary with geomagnetic activity, which could imply changes in the power law index. Plasmaspheric electron density has been observed to decrease with increased K p index (representing more disturbed conditions), which was attributed to reduced refilling rates [Young et al., 1982;Denton et al., 2002;Laakso et al., 2002a;Denton et al., 2004]. This is expected to result in decreased electron density at the field line ends, reducing the power law index during disturbed conditions. In direct contrast, Su et al. [2001] observed increased refilling rates from the ionosphere during periods of high magnetic activity. In this case, the power law index would be increased during disturbed conditions. However, this view does not consider the changes in loss processes with geomagnetic activity, which would determine the electron density along the field line away

10.1002/2017JA024171
from the ionosphere and influence the value of the power law index. Therefore, a lack of a clear dependence of the power law index on D st index could reasonably be due to the combination of multiple competing effects.
The power law function (equation (2a)) is also parametrized by the electron density at the magnetic equatorial point on the field line, n e0 , (equation (2c)). This parameter is observed to show some important dependences with D st index, which will now be discussed and compared to previous studies. The mean value of n e0 over all MLT (the offset component of equation (2c)) includes dependences on both L value and D st index. The L value variation is such that the equatorial electron density decreases with increased L value, as observed in the previous study and in agreement with multiple previous observations [Carpenter and Anderson, 1992;Goldstein et al., 2001;Sheeley et al., 2001;Denton et al., 2002Denton et al., , 2004Berube et al., 2005;Ozhogin et al., 2012]. As discussed previously by Sandhu et al. [2016a], this is expected to be due to the increased flux tube volume and length with L, so the total density of plasma along the field lines will be reduced.
This study observes that with decreased values of D st index, the electron density logarithmically decreases across all L values (illustrated by a comparison of Figures 2a-2d) and that the gradient of the L profile becomes flatter, consistent with previous findings [Chappell, 1972;Young et al., 1982;Denton et al., 2002;Laakso et al., 2002b;Denton et al., 2004;Reinisch et al., 2004;Denton et al., 2014]. This is as expected, as previous studies indicate that during the main phase of a geomagnetic storm, an increased convective electric field results in an earthward motion of the plasmapause [Grebowsky, 1970;Chappell, 1972;Chen and Wolf , 1972;Berchem and Etcheto, 1981;Horwitz et al., 1984;Loto'aniu et al., 1999;Gallagher et al., 2000;Sheeley et al., 2001;Su et al., 2001;Laakso et al., 2002a;O'Brien and Moldwin, 2003;Dent et al., 2006;Grew et al., 2007;Thaller et al., 2015]. This is thought to be caused by a southward turning of the IMF (interplanetary magnetic field) at main phase commencement, which initiates an increase in dayside magnetopause reconnection and increases the strength of the convective electric field [Cowley, 1982;Echer et al., 2008;Milan et al., 2009;Yermolaev et al., 2010]. An earthward motion of the plasmapause at decreased D st values would correspond to the electron density at each L value being decreased compared to the quiet time conditions, in agreement with these observations. An increased convective electric field would also correspond to a flattened L profile and reduced electron densities in the plasmatrough, as the convective field acts to erode the plasmatrough. Furthermore, O'Brien and Moldwin [2003] have shown that the enhanced ring current during disturbed conditions can also result in the erosion of the plasmasphere through nonconvective processes. During the subsequent recovery phase of a geomagnetic storm, characterized by increasing D st values, the plasmapause moves outward to the quiet time position at higher L values, and filling of the depleted new plasmaspheric flux tubes with cold ionospheric plasma increases the electron densities to their original level [Chappell, 1972;Horwitz et al., 1984;Comfort et al., 1988;Huang et al., 2004;Dent et al., 2006;Tu et al., 2006].
Although convective erosion is expected to dominate the loss of lower energy electrons (<100 keV), studies have shown that wave-particle interactions and magnetopause shadowing are important loss mechanisms for energetic electrons (>100 keV) that act to decrease densities during storm times [Fu et al., 2011]. The injection of particles from the plasma sheet during southward IMF conditions can lead to significant particle energy anisotropies, which are unstable to pitch angle scattering by various plasma waves, such as electromagnetic ion cyclotron, plasmaspheric hiss, and whistler mode chorus waves [Kennel and Petschek, 1966;Reeves et al., 2003;Thorne et al., 2013;Millan and Thorne, 2007;Summers et al., 2007;Fu et al., 2011]. This results in loss to the atmosphere. Magnetopause shadowing becomes important during storm times, as an enhanced ring current causes electrons to drift radially outward in order to conserve the third adiabatic invariant associated with the gradient-curvature drift motion. This process results in electrons that were previously on closed drift paths moving onto open drift paths, such that they are consequently lost to the magnetopause [Li et al., 1997;Reeves et al., 2003;Shprits et al., 2006;Millan and Thorne, 2007;Fu et al., 2011;Turner et al., 2012].
An additional contribution to decreasing electron densities with decreasing D st index could be variations in the refilling rates from the ionosphere. As previously mentioned, it has been observed that the refilling rates are lower during geomagnetically active periods [Young et al., 1982;Denton et al., 2002;Laakso et al., 2002a;Denton et al., 2004]. This change would result in flux tubes being relatively depleted compared to quiet time conditions. However, it is noted that Su et al. [2001] observed a differing dependence, where refilling rates are higher during periods of increased geomagnetic activity, contradicting these results.
Equation (2c) also includes a sinusoidal component for n e0 , demonstrating observed MLT dependences. The key feature of the MLT variation is that the peak in the n e0 parameter is observed to move from approximately dusk to noon with decreased D st values. As the previous study observed [Sandhu et al., 2016a], the MLT

Journal of Geophysical Research: Space Physics
10.1002/2017JA024171 asymmetry in n e0 is thought to be due to the presence of a "plasmaspheric bulge," such that the electron density is generally increased toward the dusk sector for average quiet conditions [Carpenter, 1966;Chappell et al., 1970;Sheeley et al., 2001;O'Brien and Moldwin, 2003;Thaller et al., 2015]. Therefore, the D st dependence suggests that this duskside bulge is observed to rotate sunward from the dusk region to noon during periods of increased geomagnetic activity. As the duskside bulge is a result of the cancelation of the corotational and convective electric fields, an increased convective field during the main phase is expected to affect features of the bulge region. Chappell [1972] explains, through a consideration of the direction of the convective flows, that an increased convective field acts to erode plasmaspheric flux tubes from the afternoon region (the erosion is most efficient in this region, where the convective flow streamlines are almost perpendicular to the plasmapause). These high-density detached flux tubes are then convected sunward, increasing the electron density in the noon sector and forming a plasmaspheric plume that extends into the morning sector [Elphic et al., 1996;Ober et al., 1997;Burch et al., 2001;Sandel et al., 2001;Su et al., 2001;Goldstein and Sandel, 2005;Grew et al., 2007;Borovsky and Denton, 2008;Walsh et al., 2013;Katus et al., 2015;Thaller et al., 2015]. Therefore, this process supports the observed D st dependence of n e0 and previous observations of the plasmaspheric bulge rotating toward noon [Nishida, 1966;Carpenter, 1970;Chappell, 1972;Higel and Lei, 1984;Moldwin et al., 1994;Elphic et al., 1996;Gallagher et al., 1998Gallagher et al., , 2000Su et al., 2001;Laakso et al., 2002b;Katus et al., 2015;Thaller et al., 2015].
As identified by Sandhu et al. [2016a], a localized peak in electron density is observed for the low-latitude region of the field-aligned distribution. This feature is not accurately described by the power law function, so a Gaussian function with an offset of n e0 is used to model this region (equation (2b)). This electron density peak is thought to be a result of solar wind-magnetosphere coupling, acting as a source of plasma in the closed magnetosphere, in addition to the ionospheric source previously discussed. The height of this peak at the magnetic equatorial point on a field line relative to the background power law value, a, is defined by equation (2e), where the sinusoidal component represents variations with MLT. The amplitude and phase of the MLT dependence were found to show no observable correlations with L value or D st index. In agreement with the previous average model, the peak of a is located in the nightside MLT sector, which supports the proposition that nightside reconnection may be acting as a plasma source in the closed magnetosphere.
The average value of a over all MLT sectors includes dependences on D st index (see equation (2e)). It is noted here that unlike Sandhu et al. [2016a], no statistical trends of a with L could be identified; this is likely an effect of separating the data into further bins. However, this parameter is observed to have a logarithmic dependence on D st index, such that a decreases with decreasing D st index. As with the previously discussed n e0 parameter, the decrease could be a result of increased convective erosion in the plasmatrough acting to deplete the flux tubes of plasma. As this plasma population is not thought to be directly originating from an ionospheric source, the decrease in a could be expected to be independent of variations in the ionospheric refilling rates. However, if the source of this plasma is from closure of open tail field lines by nightside reconnection, then refilling rates may have an effect. These open field lines will contain both ionospheric and solar wind plasma, so a reduced ionospheric source during disturbed conditions [Young et al., 1982;Denton et al., 2002;Laakso et al., 2002a;Denton et al., 2004] could result in decreased densities on the tail flux tubes. Therefore, the total input of plasma on the nightside will also be reduced.
As shown by equation (2e), the a parameter represents the electron density at the magnetic equatorial point relative to the background power law value, n e0 . Therefore, the observed electron density at the magnetic equatorial plane will be equal to a + n e0 , corresponding to the distributions shown in Figures 3a-3d. For a clearer quantitative description of the total electron density in the magnetic equatorial plane, values are plotted as a function of L in Figures 4a-4d. Each panel corresponds to field lines at a different MLT value, as labeled. To compare the total electron density during quiet and disturbed conditions, the profiles are included for D st index values of 0 nT and −100 nT, plotted as blue and red, respectively. It can clearly be seen from Figures 3a-3d and 4a-4d that the equatorial electron density decreases with decreased D st index. As previously discussed, both the n e0 and a parameters are reduced during the main phase of a geomagnetic storm, so this dependence is in agreement with the expected result. Another feature of the equatorial distributions is the strong MLT variations. During quiet times, the electron density is peaked toward the dusk sector, indicating that the MLT variation is dominated by the power law distribution (n e0 maximizes at approximately dusk, as shown by equation (2c)). However, at decreased D st index values, corresponding to high levels of geomagnetically active conditions, Figures 3 and 4 show that the peak in electron density moves from dusk toward noon. This represents the sunward motion of the plasmaspheric bulge, previously discussed, resulting from the contribution of the power law n e0 parameter. In addition, a relative enhancement is apparent at approximately dawn. This is caused by the Gaussian function (peaked in the postmidnight sector, as shown by equation (2e)) contribution becoming comparable to the background power law, acting to shift the MLT peak toward dawn.

Average Ion Mass
The field-aligned distribution of average ion mass is modeled using a power law form function, where the power law index is constrained to be negative, as discussed in section 3.2. This functional form provides a distribution where the average ion mass maximizes toward the magnetic equatorial point on the field line. As discussed in the previous study [Sandhu et al., 2016b], this represents the effects of the centrifugal force, which acts more effectively on the heavier ions [Lemaire and Gringauz, 1998;Takahashi et al., 2004;Denton et al., 2006Denton et al., , 2009]. The resulting model is defined by equation (3) and is parametrized by the equatorial average ion mass, m av0 , and the power law index, .
The average ion mass at the equatorial point on a field line, m av0 , is described by the sinusoidal function in equation (3b). Results from the previous study [Sandhu et al., 2016b] observed an MLT asymmetry in the average ion mass distribution, with enhanced values on the nightside, peaking in the evening sector. In addition, an L dependence in m av0 described increasing values of average ion mass with decreasing L. These features of the distribution were attributed to the convection of enhanced ionospheric outflows of heavy ions from high latitudes in the cusp and auroral zones, through the plasma sheet, and earthward into the nightside inner magnetosphere [Sandhu et al., 2016b, and references therein]. Mass dispersion effects cause heavier ions to be convected to lower L values compared to lighter ions, resulting in the observed L gradient of average ion mass [Chappell, 1982;Horwitz et al., 1984Horwitz et al., , 1986Roberts Jr. et al., 1987;Comfort et al., 1988;Berube et al., 2005;Darrouzet et al., 2009;Nosé et al., 2011Nosé et al., , 2015. It was also noted that an additional minor contribution to the average ion mass distribution is the presence of a heavy ion torus located just outside the plasmasphere,
The dependences of the equatorial average ion mass, m av0 , as shown by equation (3b), are illustrated in Figures 3e-3h and 4e-4h. In terms of the MLT dependence of m av0 , no coherent variations with D st index could be identified from the data. The sinusoidal variations are in good agreement with the average model, with m av0 peaking in the evening sector. It is also noted that in contrast with the previous study, no discernible L dependence in the amplitude of the MLT variations was apparent in the data analysis, and consequently, a dependence was not quantified. However, it was observed that the sinusoidal offset term in equation (3b), representing m av0 averaged over MLT variations, demonstrated significant dependences with D st index. The L gradient increases with decreased D st index, to the extent that the gradient becomes positive for active geomagnetic conditions. Therefore, m av0 decreases with increasing L value for quiet times (Figure 3h), whereas m av0 increases with increasing L value for disturbed times (Figure 3e). This is clearly illustrated in Figures 4e-4h, through a comparison of the blue and red profiles, corresponding to quiet and disturbed conditions, respectively. In addition, the intercept of the linear dependence decreases with decreased D st index (equation (3b)). The overall result is a significant increase in m av0 at the higher L values of the model L range (L ≳ 6.5) and a small decrease in m av0 at the lower L values (L ≲ 6.5), over all MLT sectors. This feature is in agreement with multiple previous studies [Young et al., 1982;Takahashi et al., 2006;Maeda et al., 2009;Nosé et al., 2009;Mouikis et al., 2010;Ohtani et al., 2011;Maggiolo and Kistler, 2014] and can be attributed to the increased O + concentration of the plasma sheet population. Increased magnetic activity causes increased heating of the atmosphere and ionosphere, a result of the dissipation of auroral currents and E × B drifts, resulting in a rise in the ion and neutral scale heights. As O + ions have a relatively low scale height, they react more strongly to changes in magnetic activity compared to other ions [Young et al., 1982]. In addition, an increase in magnetic activity can also increase ionization by auroral electrons, where the precipitating electrons are deposited at altitudes where oxygen is the dominant species [Young et al., 1982]. Therefore, both changes in the scale height and ionization by auroral electrons have effects that are most significant for O + ions [Young et al., 1982;Kronberg et al., 2012], resulting in enhanced O + concentrations for ionospheric outflows in the cusp and nightside auroral regions. The outflowing plasma at high latitudes is convected through the lobes to the plasma sheet and to the nightside inner magnetosphere [Yau et al., 1985;Lennartsson and Shelley, 1986;Ebihara et al., 2006;Kistler et al., 2006;Haaland et al., 2009;Liao et al., 2010;Yau et al., 2012;Denton et al., 2014;Maggiolo and Kistler, 2014]. This process supports the observed increase in m av0 with increased geomagnetic activity, most significantly at higher L values toward the plasma sheet (Figures 3e-3h and 4e-4h).
As mentioned previously, the heavy ion torus also contributes to the m av0 dependence on D st index, although the L range of the model covers only the outer regions of the heavy ion torus, so the contribution to the average ion mass is minor compared to the plasma sheet population. Although no CODIF data are available in this region, Figure 4 includes an extrapolation of the average ion mass model to lower L values, covering the same range of L values as the electron density model, as shown by the dotted lines in Figures 4e-4h. The extrapolated profiles aim to give some indication on the relative average ion mass values between quiet and disturbed geomagnetic conditions. It has been suggested by Nosé et al. [2011] that during the storm main phase, the strong convective fields result in the erosion of the heavy ion torus, due to the formation of a plasmaspheric plume extending to the dayside magnetopause [Elphic et al., 1996;Ober et al., 1997;Sandel et al., 2001;Goldstein and Sandel, 2005;Grew et al., 2007;Borovsky and Denton, 2008;Walsh et al., 2013] and acceleration by magnetic field dipolarization, accelerating O + ions and forming the O + -rich ring current. This feature correlates with the m av0 dependence on D st index at lower L, where values decrease for disturbed conditions (Figures 3e-3h and 4e-4h). The interactions between the expanding plasmasphere and ring current in the recovery phase of a storm are thought to repopulate the heavy ion torus, a feature of the quiet inner magnetosphere .
The modeled field-aligned distribution of the average ion mass is also dependent on the power law index, , defined by equation (3c), where the magnitude represents the gradient of the decrease in average ion mass moving away from the magnetic equator along a field line. Equation (3c) shows that the parameter is described using a sinusoidal function, to include MLT dependences, with an offset defining the value of averaged over MLT variations. As for the previous model, the parameter, averaged over MLT, is observed to linearly become less negative with increasing L value, as represented by the offset terms in equation (3c). This feature was attributed to the increasing magnetic field flux tube volume and length at increased L. Equation (3c) shows that no dependence on D st index is included for the model parameter offset term, which is due to a lack of consistent variations with D st index in the data.

Journal of Geophysical Research: Space Physics
As shown by equation (3c), the power law index, , is observed to demonstrate dependences with MLT, indicated by the sinusoidal component. The phase of the sinusoidal term is such that the most negative values at a particular L value are observed at dusk, in good agreement with the average model [Sandhu et al., 2016b]. In this region, there will be reduced upwelling of ionospheric O + ions, due to reduced solar illumination [Young et al., 1982;Lennartsson, 1989;Stokholm et al., 1989], resulting in a larger gradient along the field line. Furthermore, the convection of the O + -rich plasma sheet population acts to enhance the equatorial average ion mass for nightside MLT sectors and steepens the field-aligned gradient. The phase of the sinusoidal term is observed to show no coherent dependences on D st index, implying that the MLT location of the power law index minima is independent of geomagnetic activity.
The amplitude of the MLT dependence of is shown to depend on D st index (see equation (3c)). However, it is noted that the amplitude term does not include a linear dependence on L value, unlike the average model, as no identifiable variations were present in the data. The geomagnetic activity dependence is such that the amplitude is observed to decrease and the gradient in L becomes flatter with decreased D st index values. This corresponds to the dawn-dusk asymmetry in reducing and its variation with L also reducing. As previously mentioned, during disturbed conditions, enhanced ionospheric outflows of O + ions at high latitudes convect, through the plasma sheet, to the inner magnetosphere. This increases the average ion mass across all MLT sectors (Figure 3). Therefore, this will act to reduce the diurnal variations and decrease the amplitude of the MLT dependence of .

Cusp Enhancement
The spatial distribution of average ion mass has been shown to vary strongly with D st index, as illustrated by Figures 2a-2h and 3a-3h. The features of the average ion mass model indicate that during disturbed conditions, enhanced heavy ion outflows occur at high latitudes and strongly influence the spatial distribution of average ion mass in the closed magnetosphere. By mapping the model average ion mass values close to the ionosphere, the distributions shown in Figure 5 are provided, corresponding to the noon sector (1130 ≤ MLT < 1230). Figure 5 (left) shows the average ion mass values for a D st index of 0 nT, which corresponds to quiet geomagnetic conditions, and Figure 5 (right) shows the average ion mass values for a D st index of −100 nT, representing relatively active conditions. It has been assumed that the field-aligned dependence observed in the CODIF data set can be extrapolated along the whole field line. Although this cannot be proven to be a reasonable assumption, the purpose of the distributions shown in Figure 5 is to provide Journal of Geophysical Research: Space Physics 10.1002/2017JA024171 estimates on the relative magnitudes of the heavy ion compositions, not to represent the true values of average ion mass in the region. These figures provide some detail on heavy ion outflows on the field lines corresponding to the spatial distributions shown in Figures 2 and 3, which will now be discussed.
A comparison of the panels shown in Figure 4 indicates that for disturbed conditions, the average ion mass increases, which is a feature consistent across other MLT sectors (not shown). Furthermore, it can be seen that for quiet conditions, the values decrease with increasing latitude (Figure 5, left). In contrast, for active intervals, this dependence reverses, resulting in an increase in average ion mass with increasing latitude, as illustrated in Figure 5 (right). This provides support for the expected enhancement in heavy ion outflows at high latitudes during geomagnetically active condition [Shelley et al., 1972[Shelley et al., , 1982Lockwood et al., 1985aLockwood et al., , 1985bChappell et al., 1987;Chappell, 1988;Andre and Yau, 1997;Yau and Andre, 1997;Peterson et al., 2008;Liao et al., 2010;Li et al., 2012]. These heavy ion outflows are then convected into the inner magnetosphere, providing the spatial distributions of average ion mass shown in Figures 2e-2h and 3e-3h.

Mass Density
By combining the electron density and average ion mass empirical models, using equation (1), a model describing the variations of the total mass density along closed field lines can be inferred, which includes dependences on L value, MLT, and D st index. As discussed in section 1, this model accounts for the contributions of both the number density and ion composition of the plasma. The spatial distributions predicted by this mass density model for various values of D st index are shown in Figures 2i-2l and 3i-3l, mapped to the X-Z SM plane and T96 magnetic equatorial plane, respectively, where it is clearly illustrated that the mass density varies significantly with geomagnetic activity. Variations in the mass density with L are shown in Figures 4i-4l, where the model has been extrapolated to L = 4.5, following the extrapolation of the average ion mass model. It is emphasized that the extrapolated values may not represent the true values but are informative on possible dependences on geomagnetic activity at lower L values. The spatial distribution of the mass density is now discussed, specifically focusing on the variations with D st index. For additional details on the dependences on L and MLT, we refer to the previous studies [Sandhu et al., 2016b;, 2016a].
Figures 4i-4l show variations in the total plasma mass density with L value, where each panel corresponds to a different MLT. The first feature that is apparent from Figure 4 is that the equatorial mass density decreases with decreasing D st index, except for the morning MLT sector at high L values (see Figures 4j and 4k). This indicates that the D st variations are generally dominated by changes in the plasma number density, as the electron density is observed to decrease during geomagnetically active conditions (as discussed in section 5.1). However, for high L values in the morning sector (Figures 4j and 4k), this dependence reverses, such that the mass density along the field line for all MLT is observed to slightly increase with decreasing D st index. In this case, the ion composition is the dominant factor, where the increasing relative concentration of O + ions during disturbed conditions acts to increase the mass density despite the decreased number density. This feature may be a result of the relative enhancement in electron density at approximately dawn for disturbed conditions (Figure 3a), expected to be due to the plasma sheet population, as previously discussed. As the number density decrease is not as significant in this region, compared to lower L shells and other MLT sectors, the ion composition variations dominate, resulting in the mass density increase.
The general decrease in mass density during decreased D st index is also apparent from the spatial distributions shown in Figures 2 and 3. This feature is a somewhat unexpected result, which will now be examined in further detail. Previous statistical studies have observed the mass density to increase in this region during geomagnetically disturbed periods, in particular Min et al. [2013] and Takahashi et al. [2002Takahashi et al. [ , 2006Takahashi et al. [ , 2010. These observations were interpreted as a result of the injection of heavy ions into the closed magnetosphere during the main phase of a storm, causing high relative concentrations of O + ions. The increased average ion mass was thought to dominate over the relatively weak reductions in electron density and act to increase the plasma mass density. However, it was noted by Min et al. [2013] and Takahashi et al. [2006] that the mass density enhancements were relatively short lived and weakly correlated with D st index. The studies [Takahashi et al., 2002[Takahashi et al., , 2010Min et al., 2013] all indirectly inferred the equatorial total mass density values from spacecraft observations of toroidal standing Alfvén wave harmonic frequencies. Using the MHD wave equation [see Takahashi et al., 2010, equation (1)] combined with a numerical magnetic field model, the mass density was determined by choosing the value that gives an eigenfrequency matching the observed frequency. In order to estimate the equatorial mass density using this method, a function describing the field-aligned dependence of mass density must be assumed. These studies all used a power law function to represent the distribution along magnetic field lines; however, previous studies have shown that this form does not represent an observed peak in mass density close to the magnetic equator [Takahashi et al., 2004;Denton et al., 2006;Takahashi and Denton, 2007;Denton et al., 2009;Sandhu et al., 2016a]. The localized peak is also observed in this study. Figure 6 shows the field line profiles at multiple D st index values, comparing for different MLT and L values, and it can be seen that an equatorial peak in total plasma mass density is present. However, as the mass density contribution to determining the field line frequency is mainly at the equator, where the magnetic field is weakest, the estimated equatorial mass density should not differ significantly from the true value. In addition, these inaccuracies would not be expected to reverse the dependence on geomagnetic activity.
In contrast, other studies have observed decreased mass densities associated with geomagnetically active conditions [Menk et al., 1999;Dent et al., 2006;Denton et al., 2006;Menk et al., 2014], in agreement with this model. Denton et al. [2006] used ratios of observed toroidal standing Alfvén wave frequency harmonics with a Monte Carlo fit to infer field line distributions of mass density at different D st index values. Unlike the previously mentioned studies, the equatorial peak in the mass density distribution was accounted for, assuming a polynomial form to represent field-aligned variations. Denton et al. [2006] found that at high latitudes, away from the magnetic equatorial plane, the mass density decreased with decreased D st index. This resulted in a more peaked distribution for disturbed conditions. This feature was attributed to an increased ring current and reduced plasmaspheric density (as previously discussed in section 5.1). The equatorial mass density showed little variation with D st index, unlike the results presented here. However, it should be noted that Denton et al. [2006] compare only two samples of the mass density observed at different D st index and that the inversion technique employed requires a functional form describing the field-aligned distribution to be assumed. A study conducted by Menk et al. [1999Menk et al. [ , 2014 also observed decreased mass density during disturbed conditions. The plasmaspheric mass density was inferred from ground-based observations of Field Line Resonance frequencies for different K p values, where the K p index represents the level of geomagnetic activity. Menk et al. [1999Menk et al. [ , 2014 found that the mass density was increased during low K p intervals, corresponding to quiet conditions. This feature was attributed to the refilling of the plasmasphere in the recovery phase and following storms, which caused the mass density to be increased compared to storm times, where the plasmasphere is depleted. Furthermore, Dent et al. [2006] presented measurements of the total plasma mass density, inferred Journal of Geophysical Research: Space Physics 10.1002/2017JA024171 from ground magnetometer observations of FLR frequencies, showing an overall depletion of the total plasma mass density during a geomagnetic storm period. The analysis indicated that although an enhancement in the heavy ion population occurred, when combined with a decrease in number density, the overall effect was reduced values in plasma mass density across all L shells examined (2.3 ≥ L ≤ 6.3) during disturbed conditions. Therefore, there is evidence both in agreement and in opposition to the observed decrease in mass density at low L values, but the reasons for these differences are unclear at present.
As mass density determines the Alfvén speed, the distribution of mass density along a field line is important for controlling the frequency of field line resonances. Changes in the field line distribution of mass density with D st index, as described by this model, are now assessed. Figure 6 illustrates variations in the shape of the field line profiles with geomagnetic activity. Note that the dependence of profile shape on L and MLT during quiet conditions are examined in the previous study [Sandhu et al., 2016a] and will not be discussed in detail here, instead focusing on the variations on D st index. First, the high-latitude region (left of dashed vertical lines) will be considered, where power law dependences are observed, with the electron density model following a positive power law index (equation (2)) and the average ion mass following a negative power law index (equation (3)). It has been shown that the power law indices (equations (2d) and (3c)) do not demonstrate any strong dependences on D st index, and therefore, the difference between profiles shown in Figure 6 are predominantly due to variations in the equatorial values (equations (2c) and (3b)). It is clear from Figure 6 that the high-latitude region is characterized by a decrease in mass density values during disturbed conditions, for the majority of L values and MLT sectors. As mentioned above, this is due to the decreased number density dominating over ion composition variations for disturbed conditions. The exception is for field lines at higher L values in the morning sector, referring to the upper two left panels of Figure 6, where an increase in mass density during disturbed conditions is shown. As previously discussed, this increase is due to the ion composition contribution dominating over the number density contribution.
The lower latitude region is also seen to vary with geomagnetic activity, as shown by the profiles in Figure 6. In this region, the electron density is locally peaked above the background power law distribution, where the relative peak height reduces during disturbed conditions, as described by equation (2). The average ion mass also approaches a maximum toward the magnetic equator (equation (3)), so the resulting mass density distribution is expected to be peaked in this region, as shown by the previous study and multiple others [Takahashi et al., 2004;Denton et al., 2006;Takahashi and Denton, 2007;Denton et al., 2009]. Consistent with the previous discussion, the peak decreases during disturbed conditions across all regions, except for the morning sector at higher L values (upper two left panels in Figure 6), where the peak height slightly increases. As discussed, this feature is attributed to the magnitude of electron density depletion being reduced in this region, such that the increase in average ion mass is the dominant contribution to the total plasma mass density.

Comparison to Parametrization With K p Index
Although the dependences of mass density on geomagnetic activity have been quantified using D st index as a proxy for the level of magnetic field disturbances, it was mentioned previously (section 4) that many studies choose to use K p index instead [Menk et al., 1999;Gallagher et al., 2000;Takahashi et al., 2002;Denton et al., 2006;Takahashi et al., 2006Takahashi et al., , 2010Takahashi et al., , 2014. The key difference between D st index and K p index is attributed to the location of the magnetometer stations, from which the magnetic field disturbances are observed. Whereas D st index is a result of observations taken at low latitudes, K p index is based on observations at a larger range of latitudes. Therefore, the variations in D st index are largely dependent on the ring current intensity, and the variations in K p index are controlled by the auroral current systems in addition to the ring current.
In order to understand how using K p index as a proxy for the level of geomagnetic activity orders the data, the variations in the field-aligned distributions of electron density and average ion mass were quantified, including dependences with K p index. The corresponding mass density model was inferred, with details provided in Appendix B. These can be directly compared to the D st -dependent models. Although a detailed comparison of the quantitative differences between model parameters is outside of the scope of this study, some key features will be highlighted.
The same functional form for the electron density field-aligned distributions was chosen for both the D st -dependent model (equation (2)) and the K p -dependent model (equation (B1)), and both electron density models cover 4.5 ≤ L < 9.5. A comparison of the model parameters for the D st -dependent model and the K p model exhibits the same variations with increasing geomagnetic activity, where the variations in electron Journal of Geophysical Research: Space Physics 10.1002/2017JA024171 density are attributed to increased convective flows in the closed magnetosphere with increasingly active geomagnetic conditions. However, the magnitude of the parameter variations is observed to be greater for the D st model than the K p model. As an example, for a field line at L ∼ 5 and 00 MLT under disturbed conditions, the n e0 parameter is decreased compared to quiet conditions, by a further 5 cm −3 for the D st model than for the K p model.
The average ion mass was also modeled including dependences on both D st index and K p index, covering closed field lines within 5.9 ≤ L < 9.5. The model parameters show that the dependences in plasma ion composition with geomagnetic activity are stronger with respect to D st index, compared to K p index. For example, the increase in the m av0 , from quiet to disturbed times, is greater for the D st model than the K p model, by approximately 3 amu for a field line at L ∼ 9 and 21 MLT.
Using the same method as for the D st -dependent model, the electron density (equation (B1)) and average ion mass models (equation (B1)) are combined to infer an empirical model for the total plasma mass density, including dependences on K p index. The model describes the total plasma mass density distribution for closed field lines, over all MLT sectors in the region of 5.9 ≤ L < 9.5. As for the D st model, the combination of decreased plasma number density and increased heavy ion concentration results in decreased total mass density for disturbed conditions. This indicates the dominance of number density variations compared to ion composition. It has been identified that variations in both electron density and average ion mass are stronger when parameterized with D st index than K p index. Therefore, dependences in the total plasma mass density are correspondingly stronger for D st index compared to K p index, as expected. For example, from quiet to disturbed conditions, the total mass density is decreased by a further ∼50 amu cm −3 for the D st model than the K p model, considering a field line at L ∼ 7 and 18 MLT. It can be concluded from the comparison that the D st index is a better proxy for storm-related variations and convective flows in the closed magnetosphere than the K p index. This is thought to be due to the D st index being predominantly controlled by the ring current, so it is more effective at isolating storm-related magnetic field perturbations. In contrast, the K p index includes larger contributions from other current systems, such as auroral current systems. The inclusion of additional current systems, affected by nonstorm-related processes, acts to smooth the variations in electron density and average ion mass, when using K p index as a proxy for geomagnetic activity [O'Brien and Moldwin, 2003].

Conclusions
This study has determined empirical models for the electron density and average ion mass distribution along closed magnetic field lines, including variations with L shell, MLT, and D st index, from an analysis of data from the WHISPER and CODIF instruments on board Cluster. The data span a time interval from approximately 2001-2012 and represent the region of 5.9 ≤ L < 9.5, corresponding to the outer plasmasphere, plasmatrough, and near-Earth plasma sheet. The models are combined to infer an empirical model for the total plasma mass density in this region, which considers contributions from both the number density and ion composition of the plasma. The model provides information on how plasma mass loading varies between quiet and geomagnetically disturbed conditions and how the spatial distribution changes. Although the binning of the data sets for D st index has resulted in a loss of dependence on some parameters, possibly due to small bin sizes, the overall variations with magnetic activity are encapsulated in the model. The results illustrate the effects of various storm-related magnetospheric dynamics on the spatial distributions of the electron density, average ion mass, and total mass density. The key findings of this study include the observed decrease in mass density during disturbed conditions for the majority of field lines, a result contradicting previous observations, which indicates the dominance of changes in the number density over ion composition variations. An exception to this dependence is the observed enhancement in the total plasma mass density for the morning sector at higher L values. This feature corresponds to a relatively less dramatic reduction in number density in this region, such that the increase in average ion mass dominates. Due to the dependence of the Alfvén speed on plasma mass density, our improved understanding of mass density variations during geomagnetically varied conditions provides insight into how the propagation of wave modes, as well as a range of other processes, responds to magnetospheric activity.
Future work involves using this model to examine properties of standing Alfvén waves on closed geomagnetic field lines using time-of-flight calculations. This will further the analysis of Wild et al. [2005] and provide information on how pulsation frequencies depend on geomagnetic activity due to variations in the plasma Journal of Geophysical Research: Space Physics 10.1002/2017JA024171 mass loading. In addition, the work presented here can be continued by considering variations in the spatial distributions of electron density, average ion mass, and mass density with respect to different parameters and indices, representing dependences on solar wind conditions. Another area of future investigation involves the extension of the model, to include observations from spacecraft with different epochs and different orbital configurations. Conjunctions of spacecraft at different points on a given field line will allow us to further explore the variations in plasma density and composition along field lines and how the observations compare to the model. Furthermore, the magnetoseismology technique can be incorporated, using the realistic form for variations along field lines obtained here, to estimate the total plasma mass density and compare to the model.

Appendix A: Modeling the Field-Aligned Distributions
For each data set (WHISPER observations of electron density and CODIF observations of average ion mass), the variations along magnetic field lines are modeled. The same functional forms as used in the average model [Sandhu et al., 2016a[Sandhu et al., , 2016b are used to describe the field-aligned variations. In order to quantify the distribution, a least squares hierarchical fitting technique [Clark and Gelfand, 2006;Tabachnick and Fidell, 2006] is employed. This approach fits to the multiple levels (L value, MLT, and D st index) that the function parameter is expected to vary with, representing variations in the data set as a whole.
The initial step in analyzing the field-aligned variations of each data set is to bin the data for the normalized radial distance along the magnetic field line, R norm , which is the radial distance at which the measurement was obtained, R (R E ), divided by the L value of the field line. The T96 magnetic field model is used here to determine the field line corresponding to the measurement, and a bin width of 0.05 is used for the normalized radius. The data are also binned for each level included in the multilevel fitting technique. The L value bin width used is 1, and the MLT bin width used is 3 h. The data are also binned into six D st index bins, where the bin values are indicated by the different colors in Figure 1.

A1. Electron Density
By plotting the electron density, n e , as a function of normalized radius, R norm , the field-aligned dependence is illustrated. An example plot is shown in Figure A1, for 8.5 ≤ L < 9.5 and 15 MLT, where the points indicate the average electron density value in the normalized radius, R norm , bin. The profiles shown have been smoothed using a boxcar function, with a width of 3. Each panel corresponds to a different bin of D st index, referring to Figure 1. As previously discussed, it has been shown that the field line dependence of electron density can be separated into two regions, where the boundary between these dependences is defined at R norm = 0.8 (indicated by the vertical grey dashed line in Figure A1). At high latitudes (R norm ≤ 0.8) the electron density generally increases toward the ends of the field line, corresponding to a power law distribution. At low latitudes (R norm > 0.8), a peak in electron density close to the magnetic equatorial plane is often observed, which is represented by a Gaussian distribution. Using the same functional forms as the average model, a least squares fitting method, weighted by the number of points in each bin, is used to determine the best fit function parameters for each field-aligned profile. The function parameters are the electron density at the magnetic equatorial point on the field line, n e0 , the power law index, , and the peak height above n e0 , a. The hierarchical approach includes dependences on L, MLT, and D st index in the best fit parameters, providing the required functions. The power law model is given by equation (2a), and the Gaussian model is given by equation (2b). The model parameters (n e0 , , and a) are defined by equations (2c)-(2e). The field line distribution predicted by the model is shown by the overplotted solid colored lines in Figure A1. The empirical electron density model describes variations in the L range of 4.5 ≤ L < 9.5 and covers all MLT.
The functional form used to describe the model parameters in this study are chosen to include a sinusoidal term, to represent the MLT dependence. The phase of the sinusoidal term defines the location of the peak of the parameter in degrees of MLT eastward from the midnight meridian. The sinusoidal term also includes an amplitude to represent the magnitude of the MLT dependence. The function contains an offset added to the sinusoidal term, which defines the mean value of the parameter across all MLT. The amplitude, phase, and offset parameters are all allowed to be functions of L value and D st index. The choice of functional forms used for the electron density model was such that they are the simplest forms that described the observed variations in the data, minimizing the number of free parameters in the fitting process. It is also noted that Journal of Geophysical Research: Space Physics 10.1002/2017JA024171 Figure A1. Electron density, n e , (cm −3 ) plotted as a function of normalized radius, R norm , for 8.5 ≤ L < 9.5 and 15 MLT. The panels display data for each D st index (nT) bin (as defined in Figure 1). The vertical dashed line indicates the boundary between the power law and Gaussian dependences, and the overplotted colored lines represent the best fitting functions to the overall data set. The upper and lower quartiles of the distribution of points averaged in each bin are shown by the grey line, intersected by a short horizontal line at the median value.
when fitting for the model parameters, if no dependence on L or D st was clearly observed in the data, the dependence was removed from the relevant functional form.

A2. Average Ion Mass
Using the same method as for the electron density model discussed above (section A1), the field-aligned distribution of average ion mass is now examined. Figure A2 shows an example plot for 8.5 ≤ L < 9.5 and 15 MLT, where each panel shows the average ion mass, m av , as a function of normalized radius, R norm (smoothed using a boxcar function of width 3), for each D st index bin. Note that unlike the electron density plots in Figure A1, logarithmic scales are used for both axes, which linearizes power law dependences. As discussed in section 1,

Journal of Geophysical Research: Space Physics
10.1002/2017JA024171 Figure A2. Average ion mass, m av , (amu) plotted as a function of normalized radius, R norm , for 8.5 ≤ L < 9.5 and 15 MLT. The panels display data for each D st index (nT) bin (as defined in Figure 1). The overplotted colored lines represent the best fitting functions to the overall data set. The upper and lower quartiles of the distribution of points averaged in each bin are shown by the grey line, intersected by a short horizontal line at the median value. the field line dependence of average ion mass can be represented by a power law function with a negative power law index, describing the average ion mass maximizing toward the magnetic equator and decreasing off equator. A function of this form is least squares fitted to the data, with the fitting weighted by the number of data points in each bin. Using the hierarchical fitting approach, the function parameters (the average ion mass at the magnetic equatorial point of the field line, m av0 , and the power law index, ) are obtained, with dependences on L value, MLT, and D st index included. The resulting model for the field-aligned distribution of average ion mass is given by equation (3), with the function parameters defined by equations (3b) and (3c). The model field-aligned variations are illustrated by the overplotted solid colored lines in Figure A2. The model covers all MLT values in the region of 5.9 ≤ L < 9.5.

Appendix B: Dependences on K p Index
Using the same approach as for the D st -dependent models, dependences on K p index in the electron density and average ion mass data sets were quantified. The WHISPER electron density and CODIF average ion mass data sets were binned for K p index, using a bin size of 1, where data with a K p index greater than 4 were not included due to the low number of data points available. The same hierarchical modeling technique, as detailed in sections 3.1 and 3.2, was applied to the field-aligned distributions of each data set, except that Journal of Geophysical Research: Space Physics 10.1002/2017JA024171 variations with K p index were modeled, as opposed to variations with D st index. Combining the models, using equation (1), provides the resulting mass density model, which includes dependences on L, MLT, and K p index and represents all MLT sectors for the closed magnetosphere between 5.9 ≤ L < 9.5.

B1. Field-Aligned Distribution of Electron Density
The resulting model for the field-aligned distribution of electron density, n e , is shown in equation (B1). The electron density model covers all MLT sectors in the region of 4.5 ≤ L< 9.5. The same functional form as before was used (equation (2)), where a power law distribution represents the higher-latitude region of field lines (equation (B1a)), and a Gaussian function is used to describe the equatorial region (B1b). The model parameters for the electron density model (equatorial electron density, n e0 , power law index, , and relative peak height, a) were quantified to include dependences on L, MLT, and K p index and are defined by equations (B1c)-(B1e).

B2. Field-Aligned Distribution of Average Ion Mass
The corresponding model for field-aligned variations in average ion mass, m av , using the same power law functional form as previously (equation (3)) is given by where m av0 is the average ion mass at the magnetic equatorial point of the field line and is the power law index. The model parameters (m av0 and ) include dependences on L, MLT, and K p index and are described by equations (B2a) and (B2c). The average ion mass model covers all MLT sectors in the region of 5.9 ≤ L < 9.5.