Time-varying magnetotail magnetic flux calculation: a test of the method

We modified the Petrinec and Russell (1996) algorithm to allow the computation of time-varying magnetotail magnetic flux based on simultaneous spacecraft measurements in the magnetotail and near-Earth solar wind. In view of many assumptions made we tested the algorithm against MHD simulation in the artificial event, which provides the input from two artificial spacecraft to compute the magnetic flux F values with our algorithm; the latter are compared with flux values, obtained by direct integration in the tail crosssection. The comparison shows similar time variations of predicted and simulated fluxes as well as their good correlation (cc>0.9) for the input taken from the tail lobe, which somewhat degrades if using the “measurements” from the central plasma sheet. The regression relationship between the predicted and computed flux values is rather stable allowing one to correct the absolute value of predicted magnetic flux. We conclude that this method is a perspective tool to monitor the tail magnetic flux which is one of the main global magnetotail parameters.


Introduction
The magnetic flux circulation in the solar windmagnetosphere system determines the dynamical regime of magnetosphere (Russell and McPherron, 1973).This regime depends on the balance between the magnetic flux, reconnected on the dayside magnetopause and coming to magnetotail from the solar wind, and the magnetic flux, reconnected in the magnetotail and transported to the dayside Correspondence to: M. A. Shukhtina (mshukht@geo.phys.spbu.ru)magnetosphere.If the dayside reconnection rate exceeds the tail reconnection rate, magnetic flux is stored in the tail and the substorm growth phase is observed; the situation when the tail reconnection prevails corresponds to substorm expansion phase; if there is a balance between dayside and nightside reconnection rates, steady magnetospheric convection is realized.Thus monitoring of the magnetotail magnetic flux value is necessary to characterize the state of the magnetospheric system.However, until recently the knowledge about this quantity was very poor.The reason is that F is a global characteristic, which is difficult to infer from local observations.Several approaches are possible to calculate the magnetic flux F .First, its value may be determined from optical observations of the polar cap (PC) area, as the PC magnetic field lines are believed to project to tail lobes.Recently routine optical observation of large polar areas from Polar and Image spacecraft made available PC observation during long periods (Brittnaher et al., 1999;Hubert et al., 2006;Milan et al., 2007;DeJong et al., 2007).However, weak luminosity near the PC boundary combined with dayglow contamination may influence the accuracy of the PC boundaries determination, the accuracy itself being difficult to estimate.
An alternative approach to magnetotail flux calculation was proposed by Petrinec and Russell (1996) (further PR96).Their method is based on tail lobe magnetic field observations combined with time-shifted simultaneous solar wind measurements.Petrinec and Russell (1996) derived a statistical formula to describe the tail radius R T as a function of the solar wind parameters and magnetotail spacecraft position, allowing to compute the flux value as F =0.5 π R 2 T B L , where B L is the measured lobe magnetic field.The analysis of this R T model showed its good correspondence with both the experiment and other models (Shue et al., 1998;Kawano et al., 1999;Yang et al., 2002).However, all these magnetopause models depend only on the external parameters and are therefore the same for different magnetospheric states.
Meanwhile it has been established long ago that the tail radius (as well as B L and F values) depends also on the state of magnetosphere: all three variables grow during the substorm growth phase and decrease after the substorm onset (e.g.Caan et al., 1973;Maezawa, 1975;Baker et al., 1994).The dependence of B L and R T values on the magnetospheric state as well as their dependence on external parameters was studied by Shukhtina et al. (2004) (further Sh04).They obtained the regression models for the F value as a function of solar wind parameters and X-coordinate in the tail for three different magnetospheric regimes: Quiet (Q), Steady Magnetospheric Convection (SMC) and Substorm Onset (SO) separately.Their results showed that the magnetic flux values for all three states are independent of the dynamic pressure Pd and almost independent of X, confirming that F can be treated as a state variable for the magnetotail.While the Sh04 results displayed a substantial flux differences for the states considered, their approach still used a fixed model for any particular state and did not allow to follow the magnetopause changes and to compute the true F variations.In the present paper we propose and test a modification of the PR96 algorithm, allowing one to calculate the actual time-dependent magnetic flux in the magnetotail.
Like the original Petrinec and Russell (1996) method, the modified algorithm (hereafter referred to as the SPR algorithm) uses a number of serious assumptions and approximations, and thus needs a very detailed testing.Since the method is based on the MHD equilibria equations, the global MHD modeling provides the best and natural possibility for its verification.In this paper we compare F values, predicted by the SPR algorithm for different locations of the artificial spacecraft in the tail, with those obtained by direct integration of the magnetic flux through the tail cross-section.To make certain that the SPR scheme really improves the original PR96 one, the computations using the PR96 algorithm are also presented.

Description of the modified (SPR) algorithm
The method is based on the approximation of axisymmetric flaring magnetotail and uses the equations of onedimensional pressure balance in the magnetotail (vertical balance across the current sheet) and across the flaring magnetopause: where B L is the equivalent lobe magnetic field, and α is the flaring angle.The coefficient 0.88 in Eq. ( 2) is the ratio of the magnetosheath pressure to solar wind dynamic pressure for high Mach numbers (Newtonian approximation, Spreiter et al., 1966).If the α value is known, one can calculate the tail radius R T , which is necessary to compute the tail magnetic flux.
The validity of Eq. ( 1) was tested experimentally in the midtail (Fairfield et al., 1981;Baumjohann et al., 1990;Petrukovich et al., 1999), and it was concluded that the vertical balance is approximately satisfied at distances tailward of X∼−15 R E (with some exceptions near substorm onset - Petrukovich et al., 1999), where the "tail approximation" is fulfilled.Magnetic and plasma data on the RHS of Eq. ( 1) is taken from results of the MHD simulation.
The pressure balance across the magnetopause was studied and confirmed experimentally in PR96.In Eq. ( 2) parameters on the LHS are the solar wind parameters, time-shifted to the X coordinate of the observational point.For simplicity we use the time shift T = X/Vsw, where X is the distance along X between the observational points in the tail and in the solar wind.The electron temperature Tesw is assumed equal to the ion temperature Tisw.Comparison of results based on this assumption with those for Tesw=2Tisw (which is perhaps more appropriate -see Newbury et al., 1998) showed that they are almost identical due to the relatively small contribution of the thermal pressure to the LHS of Eq. ( 2).The solar wind dynamic pressure was calculated as Pd=1.94×10−6 n SW V 2 SW (assuming 4%-helium content, e.g.Tsyganenko, 2002).Solving Eq. ( 2), we determine the flaring angle α value, necessary for the tail radius R T computation.
When calculating the R T value, the magnetosphere is assumed axisymmetric (relative to the X axis, the GSM coordinates are assumed everywhere).As tan α=dR T /dx, the R T value may be calculated as where R T 0 is the tail radius value at the terminator (X=0).
As shown in Petrinec et al. (1991), the radius at terminator is almost independent of IMF B Z value, and following their results we calculate the R T 0 value as with Pd given in nPa, and X, R T and R T 0 in R E .
Here begins the difference of the present (SPR) procedure from the PR96 approach, where (as well as in Sh04) a large set of Eq. (2) for different conditions was solved to obtain model formulas for α and R T dependence on inputs Pd, IMF B Z and spacecraft X-coordinate.After sin 2 α(X) dependence is determined, the integral in Eq. ( 3) may be calculated, giving R T and F values.One function for all situations was derived in PR96, whereas in Sh04 different formulas for different states (Q, SMC and SO) were obtained.In the present approach the α, R T , and F values are calculated from Eqs. (1-3) for every single measurement in the tail and solar wind, that is for any specific time.To use Eq. ( 3) one needs to specify the α(X) dependence.Based on previous experience, we express sin 2 α(X) as According to Sh04 this dependence fits well every of three states considered, the C values for Q, SMC and SO being 0.0749, 0.0781 and 0.0612.For further calculations we use the average of these values, C=0.0714.Substitution of Eqs. ( 4) to (3) then gives: where A is determined from Eq. ( 4).
The last thing to do is to take into account the B L and α distribution in the tail to match the observation point inside the tail with the magnetopause point where the balance is evaluated.One may expect that the B L -isocontour lines should be perpendicular to magnetopause in its vicinity, but should follow the lines X=const in the central sector of the tail (this was confirmed by ISEE2 spacecraft observations, see e.g.Fig. 3 in PR96).So the X value at the spacecraft (e.g., Geotail) position should be replaced by some new X * value at the magnetopause which has the same ) sinαcosα, X, Y, Z being the coordinates of the observation point (Fig. 1, the same as Fig. 2 in PR96).Substituting X * to Eq. ( 4), we obtain the new A* value, which is then substituted to Eq. ( 5) to get the new R T value, the new X value, etc.After 3-4 iterations the solution converges, and finally Finally, ignoring the depressed (compared to the lobe values) magnetic field inside the plasma sheet and suggesting the circular tail cross-section, we approximate the tail magnetic flux as This procedure is repeated for any time step.Thus the (modified) SPR algorithm of the magnetotail flux calculation requires a) knowledge of the B L value at some point tailward −15 R E (calculated from Eq. 1), and b) knowledge of the sin 2 α value, which, according to Eq. ( 2), needs data on solar wind parameters Pd, n SW , Tsw and Bsw.The solar wind parameters are inputs for MHD simulation, whereas plasma and magnetic pressures in the magnetotail are the results of the simulation.
Assumed pressure balance in the magnetotail and at the magnetopause (Eqs. 1, 2), the presumed axisymmetric magnetopause shape (Eq.6) and the neglected plasma sheet (Eq.7) are all the approximations that may lead to systematic errors.The MHD simulation provides an excellent opportunity to test and correct the used approach.To do it, F values, Fig. 1.The scheme, presenting the geometry of Geotail measurements.X is the Geotail position, X * is the magnetopause coordinate, corresponding to the same B L and α values, R T is the tail radius.
calculated by the SPR algorithm, are compared with the integral of the magnetic flux through a magnetotail cross-section tailward −15 R E (Flux Direct, F D ).

Direct magnetic flux calculation
We have run the simulation at the Community Coordinated Modeling Center (CCMC), operating at NASA GSFC.Standard global MHD code OpenGGCM was used, which solves the MHD equations with additional dissipation (see, e.g., Raeder, 2003) in the simulation domain [24, −350 Stretched Cartesian grid with step size in X,Y and Z, changing from (0.25, 0.4, 0.25) R E in the plasma sheet to (0.25, 0.75, 1.1) R E at the magnetopause at X distances from 10 R E to −30 R E was used.
As our goal was tail magnetic flux F calculation, the most important thing was the accurate magnetopause identification.Different approaches to magnetopause determination have been discussed in the literature.We tried three different methods, based on density gradient, current density peak (Garcia and Hughes, 2007), and fluopause (Palmroth et al., 2003) determination.Fluopause is defined as a family of plasma streamlines, starting in the solar wind and passing most close to the x-axis on the nightside -see Fig. 2c  and d.The streamlines are started from X=12 R E with a step size 0.5 R E in Y and Z. the fluopause method is used hereafter as a proxy of the magnetopause.
To calculate the flux value, the YZ cross-sections of magnetosphere inside the fluopause at fixed X were considered.The area of each square cell 0.5 R E ×0.5 R E was multiplied by the B X value in the cell center.For cells, crossing the magnetopause, the area of the part inside the fluopause was taken.The integral of these elementary flux values gives the total magnetic flux value through the given cross-section F D , the flux calculated by the direct integration.It is used for testing the F values, predicted by the SPR algorithm, which is described in the previous section.The F values were com-puted separately for Northern and Southern lobes, the difference between them indicating the accuracy of F D calculations.In addition, using direct integration, we could evaluate the relative contributions to flux values from the lobes and the plasma sheet, which is ignored in the SPR.

Results
The flux values, predicted by SPR and by original PR96 algorithm, were compared with results of direct flux calculation for the CCMC run "Victor Sergeev 010907 1".shown in Figs. 3, 4. The North-South-North IMF turning sequence allows to simulate the effects of a substorm, whereas the V Z variations should change the position of the magnetotail neutral sheet (originally this simulation was undertaken to study the effect of V Z variations on the neutral sheet position -see Tsyganenko and Fairfield, 2004).Two magnetotail cross-sections, at X=−15 R E and X=−25 R E , were chosen.Observation points in these cross-sections (locations of artificial spacecraft in SPR) were taken in the lobes at Y =4 R E , Z=±10 R E (Fig. 3a and c). the initial value −4 nT to +4 nT.The subsequent 51-min interval of B Z =−4 nT (Vsw=600 km/s) was concluded by F D decrease from ∼0.9 GWb to ∼0.65 GWb following the IMF northward turning.During the negative B Z interval (growth phase) F D increased from 0.78 to 0.91 GWb (not monotonously, with a dip to 0.87 GWb, perhaps a pseudobreakup).V Z variations, switched on after the simulated substorm, are accompanied by low-amplitude F D oscillations.
The predicted F demonstrates all kinds of variation, displayed by direct flux calculation.However, the F values, given by SPR, exceed the F D ones by 0.1-0.2GWb, probably due to the fact, that we ignored the existence of the plasma sheet, threaded by a smaller magnetic flux.The PR96 algorithm uses the mean model R T (B Z ) dependence, this is why the predicted F values during the growth phase are lower than those provided by both F D and SPR calculations.
Figure 3c displays flux variations at X=−25 R E .According to all methods the flux values at X=−25 R E are very close to those at X=−15 R E , being slightly (by ∼10%) lower due to partial magnetic flux closure across the equatorial current sheet.All methods demonstrate periodic F oscillations (synphase in both hemispheres) during the interval of V z variations.Though the origin of these variations is unclear, they do not affect our testing of the methods.
Figure 3b and d demonstrates the relationship between F values, calculated by either SPR or PR96 algorithms, and F D values (here we use the average of the North and South lobe fluxes).The correlation is better for the SPR method, than for PR96 one (cc=0.99versus 0.92 for X=−15 R E , and 0.93 versus 0.77 for X=−25 R E ).Averaging over two crosssections gives the regression equations F =0.9 F D +0.2 (SPR) and F =0.3 F D +0.5 (PR96) -Fig.3b, d.Interestingly in both figures the slope of the PR96 regression line is close to that of SPR one below F D ∼0.65-0.7 GWb, flattening at higher flux values; so in this simulation the PR96 algorithm is not appropriate for high flux values.In whole the SPR regression is more realistic, its slope being closer to 1 with smaller free term.From Fig. 3 we conclude, that for artificial spacecraft in the tail lobes (at |Z|=10 R E , where plasma beta<0.01) the SPR method demonstrates a good correspondence with direct flux calculations, giving higher correlation and more realistic regression equations compared to PR96.
However, in many cases the magnetotail spacecraft, such as the Geotail spacecraft in recent years, is situated in the plasma sheet.Therefore we repeated the calculation for the same run, but putting the artificial spacecraft in the plasma sheet.In the top of Fig. 4 we present results of calculations at X=−15 R E , Y =4 R E for the "measurements" in the neutral sheet (a) and in the plasma sheet, where plasma beta values are 1, 0.5 and 0.1 (c) (all points are within 5 R E from the neutral sheet).Now the correspondence with F D is worse than in Fig. 3.The flux oscillations, predicted by empirical algorithms, exceed those computed directly (F D ).Note, that the PR96 model was created for lobe spacecraft, underestimating the difference between the neutral sheet space-craft X coordinate and magnetopause X* value at the same B L isoline compared to SPR; this is probably the reason of PR96 flux values exceeding the SPR ones in Fig. 4 contrary to Fig. 3 (lobe spacecraft).The cc values for the neutral sheet are 0.57 and 0.68 for PR96 and SPR algorithms correspondingly (Fig. 4b).In the plasma sheet the SPR algorithm gives cc=0.70 for beta=1 and 0.5, and cc=0.79 for beta=0.1 (Fig. 4d).In spite of worse correlation the regression equations for all spacecraft positions are similar to those obtained for the lobes.Particularly, from Figs. 3, 4 we conclude, that for beta <1 the SPR algorithm gives a rather stable regression equation: F =0.8 F D +0.2 GWb (F =0.7 F D +0.2 GWb in the neutral sheet).

Discussion
We propose a modification of the PR96 method, allowing one to compute the time-varying magnetic flux in the magnetotail.The modification is as follows: whereas in PR96 the tail radius R T is calculated from statistically obtained model formulas, which ignore the magnetotail magnetic flux change during substorms, in our approach (SPR) the R T value is calculated from the pressure balance on magnetopause for every time step.The lobe magnetic field in both cases is taken from measurements (though in our calculations we also tried B L estimates obtained from the measurements in the plasma sheet).The method is based on simultaneous measurements in the tail and solar wind at every time step; it uses many assumptions to make the problem tractable (neglects the y−z magnetopause asymmetry and plasma sheet existence, assumes a specific shape of B L isolines, uses the simplified formula for the tail radius value at terminator, etc.) and, thus, requires a thorough testing.It is also necessary to compare both SPR and PR96 results with independent magnetic flux measurements to understand if our modification of PR96 really gives a serious improvement.
Such independent magnetic flux "measurements" were taken from global magnetospheric MHD modeling, which allows to calculate for the given solar wind input all magnetospheric parameters of interest, in particular, the value of the magnetotail magnetic flux.We predicted the magnetic flux values, using different positions of artificial spacecraft, and compared them with the calculated F D values for two tail cross-sections, at X=−15 R E and at X=−25 R E .When the artificial spacecraft is situated in the tail lobes (at Z=10 R E , where beta<0.01), it gives the best prediction of tail magnetic flux values (Fig. 3), with high correlation coefficient (cc>0.9), the slope of the regression line being ∼0.8-0.9 with a small free term in the regression equation.For original PR96 the correlation is lower and the regression coefficients are much smaller, about 0.3-0.4,with a large free term (0.5 GWb against 0.2 GWb for SPR).That means that the PR96 algorithm predicts only 30-40% of the real F changes.When the artificial spacecraft "moves" to the plasma sheet, www.ann-geophys.net/27/1583/2009/Ann.Geophys., 27, 1583Geophys., 27, -1591Geophys., 27, , 2009 the correlation coefficients decrease, reaching minimum in the neutral sheet (0.7 and 0.6 for the SPR and PR96 correspondingly, Fig. 4).According to both algorithms the F variations are much larger, than the F D ones.However for all spacecraft positions the SPR algorithm gives nearly the same regression equation F (GWb) = 0.8F D +0.2 GW.The nonzero free term probably results from neglecting the plasma sheet existence in the algorithm; it is also probably the reason for the SPR values being on average larger (by ∼10-20%), than the F D ones.According to MHD simulation the average plasma sheet (inside beta>1) contribution to the tail magnetic flux is about 10%.The correlation coefficients grow with beta decrease from cc=0.7 in the neutral sheet to cc=0.8 for beta=0.1 and to cc>0.9 for beta <0.01.
Corresponding standard deviations are ∼0.1 GWb (>10%) in the neutral sheet and ∼0.01 GWb (∼1%) in the lobes.The stable regression equation makes possible to correct the predictions as: F corr (GWb)=1.25F-0.25 in future applications.
The dependence of the calculated F values on the spacecraft position needs a special study; it may result e.g. from the deviations from 1-dimensional geometry in the presence of the magnetic flux closure through the plasma sheet (e.g.BBFs, plasmoids etc).We plan to explore it in the future studies.
In the present simulation a limited interval of solar wind/IMF parameters is considered.To do more general conclusions we plan to study a wide range of input parameter variations.

Conclusions
An algorithm, suitable to compute the tail magnetic flux at any specific time, which is necessary for magnetospheric dynamics monitoring, is proposed.The predictions of both modified (SPR) and original PR96 algorithms are compared with independent flux estimates F D , obtained by direct integration of magnetic flux values through the tail cross-section in MHD-simulated magnetosphere.The test, based on MHD simulation of a 5-h interval with changing solar wind/IMF parameters, showed a good predictive efficiency of the SPR algorithm compared to the PR96 one, especially for the observational point, taken in the tail lobe (cc>0.9 for SPR against 0.8-0.9 for PR96 at Z=10 R E , where plasma beta is <0.01).The correlation coefficients are weaker in the neutral sheet (0.7 and 0.6 for SPR and PR96 correspondingly), but the regression equation appears to be almost independent of beta value (F =0.8 F D +0.2 GWb against F =0.4 F D +0.5 GWb), giving the opportunity to correct the SPR predictions.

Fig. 2 .
Fig. 2. The positions of magnetopause surfaces, corresponding to different definitions, in the XZ cross-section for Y =0 R E (a) and YZ cross-section for X=−16 R E (b); plasma streamlines (c) and fluopause surface (d).

Fig. 3 .Fig. 4 .
Fig. 3. Top: Magnetic flux values, calculated by different methods at X=−15 R E (a) and X=−25 R E (c) for "observational points" in the tail lobes (Z=±10 R E ).Solar wind B Z and V Z variations, shifted to the observational point, are also shown.Bottom: comparison of flux values, calculated by different methods for X=−15 R E (b) and X=−25 R E (d).
Results for the North and South lobes are shown by a solid and dashed line correspondingly.The flux values in two lobes are very close: mean relative difference between North and South F D values is 0.4% and 3% at X=−15 R E and X=−25 R E correspondingly; the SPR algorithm gives 0.3% and 0.2% correspondingly.Consider first the F D variations at X=−15 R E , Fig. 3a.Initially (till B Z southward turning) the F D value slightly (by ∼0.1 GWb) decreased, following the B Z change from Ann. Geophys., 27, 1583-1591, 2009 www.ann-geophys.net/27/1583/2009/