Turbulence in a global magnetohydrodynamic simulation of the Earth’s magnetosphere during northward and southward interplanetary magnetic field

We report the results of MHD simulations of Earth’s magnetosphere for idealized steady solar wind plasma and interplanetary magnetic field (IMF) conditions. The simulations feature purely northward and southward magnetic fields and were designed to study turbulence in the magnetotail plasma sheet. We found that the power spectral densities (PSDs) for both northward and southward IMF had the characteristics of turbulent flow. In both cases, the PSDs showed the three scale ranges expected from theory: the energy-containing scale, the inertial range, and the dissipative range. The results were generally consistent with in-situ observations and theoretical predictions. While the two cases studied, northward and southward IMF, had some similar characteristics, there were significant differences as well. For southward IMF, localized reconnection was the main energy source for the turbulence. For northward IMF, remnant reconnection contributed to driving the turbulence. Boundary waves may also have contributed. In both cases, the PSD slopes had spatial distributions in the dissipative range that reflected the pattern of resistive dissipation. For southward IMF there was a trend toward steeper slopes in the dissipative range with distance down the tail. For northward IMF there was a marked dusk-dawn asymmetry with steeper slopes on the dusk side of the tail. The inertial scale PSDs had a dusk-dawn symmetry during the northward IMF interval with steeper slopes on the dawn side. This asymmetry was not found in the distribution of inertial range slopes for southward IMF. The inertial range PSD slopes were clustered around values close to the theoretical expectation for both northward and southward IMF. In the dissipative range, however, the slopes were broadly distributed and the median values were significantly different, consistent with a different distribution of resistivity.


Introduction
Several studies have found evidence for turbulence in the plasma sheet (Consolini et al., 1996(Consolini et al., , 1998Borovsky et al., 1997;Borovsky and Funsten, 2003;Lui, 2001Lui, , 2002Weygand et al., 2005Weygand et al., , 2006Weygand et al., , 2007. Borovsky et al. (1997) and Borovsky and Funsten (2003) examined the temporal fluctuations of the plasma sheet magnetic field and flow velocity. They found power spectral indices in the range of −0.8 to −2.0 for the flow and −1.6 to −3.0 for the magnetic field. These ranges are quite large and include both the Kolmogorov (1941) and Kraichnan (1965) theoretical values. However, Weygand et al. (2005) argued that power spectral indices do not unambiguously establish the presence of turbulence within the plasma sheet. They also used probability distribution functions (PDFs) to examine the possibility of intermittent turbulence, which results from the spatial concentration of dissipation (Biskamp, 1993). Weygand et al. (2005) found that the range of magnetic field spectral indices was quite broad, averaging −2.0 ± 0.4. The transverse magnetic field component (approximately in the y-direction) had a spectral index of −1.56 ± 0.04, very close to the theoretical value of Kraichnan (1965). Evidence for turbulence is not limited to the plasma sheet. Nykyri et al. (2006a) have used data from the Cluster spacecraft in the cusp to construct magnetic field power spectra. They found slopes between −2.7 and −1 in the inertial range and between −5 and −3 in the dissipative range. They concluded that the cusp was turbulent too.
Several researchers have constructed theoretical models of the plasma sheet with turbulence included (e.g., Antonova, 2000;Antonova and Ovchinnikov, 1999;Borovsky et al., 1998), which suggest that eddy diffusion can fundamentally influence the plasma sheet configuration. Borovsky et  Borovsky (2004) suggested that large-scale shear flows at the magnetopause are a source for turbulence. Both local (Nykyri and Otto, 2001;Nykyri et al., 2006b) and global MHD simulation studies have shown vortices forming at the magnetopause (e.g., Walker et al., 1998Walker et al., , 2006Ashour-Abdalla et al., 1999, 2002White et al., 2001;Slinker et al., 2003;Hasegawa et al., 2004;Fairfield et al., 2000;Collado-Vega et al., 2007;Claudepierre et al., 2008;Hwang et al., 2011). Vortices at the magnetopause have been found for both northward and southward IMF and have been interpreted as Kelvin-Helmholtz waves (Walker et al., 1998;Nykyri and Otto, 2001;Slinker et al., 2003;Nykyri et al., 2006b;Fairfield et al., 2007;Collado-Vega et al., 2007;Claudepierre et al., 2008). Other simulations have featured vortices within the plasma sheet that were not associated with boundary oscillations (El-Alaoui, 2001;White et al., 2001;Ashour-Abdalla et al., 2002;Walker et al., 2006;El-Alaoui et al., 2009. Recently El-Alaoui et al. (2010) presented high resolution global MHD simulation results that showed turbulent magnetotail flows for southward IMF conditions. The study found that a vigorous spectrum of fluctuations and eddies occurred under steady southward IMF driving conditions. They computed PSDs and PDFs and found results consistent with insitu observations and with theory. They concluded that localized reconnection was the major process driving the turbulence. The state of the magnetosphere is different during northward and southward IMF, and both the driving and dissipative processes for turbulence are expected to be different. In this paper we investigate these differences.

Simulation model
The coupled magnetosphere-ionosphere, three-dimensional global MHD code used for this study is based on a onefluid description of the interaction between the solar wind and Earth's magnetosphere (Raeder et al., , 1998(Raeder et al., , 2001Frank et al., 1995;El-Alaoui, 2001;El-Alaoui et al., 2004). The numerical resistivity in the code is so low that an anomalous resistivity model must be introduced (Raeder et al., 2001). The ionospheric part of the model takes into account three sources of ionospheric conductance: solar EUV ionization, diffuse auroral precipitation, and the electron precipitation associated with upward field-aligned currents. A detailed description of the MHD model can be found in Raeder et al. (2001). For the results presented in this paper a combination of low resistivity and small grid spacing was used to achieve a high magnetic Reynolds number that allowed turbulence to develop. A viscous interaction at the magnetospheric boundary in the MHD simulation yields an electric potential of about 30 kV, which is consistent with the observed value (e.g., Reiff et al., 1981), even though it is purely numerical.
In the MHD simulation, the total electric field includes convective and resistive terms, and is given by: where v is the bulk flow velocity, B is the magnetic field, J is the current density, and η is the resistivity. Explicit resistivity is necessary in our code for reconnection to occur (Raeder et al., 2001). The resistivity η in the code is determined by applying the equations below.
In this equation j is the local current density; B is the local magnetic field, is the grid spacing and ε is a small number added to avoid dividing by zero. The coefficient α is determined empirically and is much less than 1. Similar models based on current-driven instabilities have been used successfully in local MHD simulations (Sato and Hayashi, 1979). To avoid spurious dissipation there is a minimum current density threshold for resistivity controlled by the parameter δ (a value of 0.65 was used for these results) (Raeder et al., 2001). This threshold is calibrated such that explicit resistivity is switched on only at a very few grid points in intense current sheets. Further assessment of anomalous resistivity on substorm development can be found in Raeder et al. (2001). We set α = 0.05, and the resulting maximum value for η was about 10 6 × m. Usually the MHD simulation results for the entire system are written out at intervals of 20 to 60 s. To obtain high time resolution time series of the simulation results for use in calculating power spectra and probability distribution functions this was supplemented by writing out all MHD variables at all grid points in the equatorial plane every second.
The simulation used a grid that had very high spatial resolution in the plasma sheet, as discussed in El-Alaoui et al. (2010). The simulation used more than 52 million grid points: 540 in the x-direction, 324 in the y-direction and 300 in the z-direction. The grid distribution is illustrated in Fig. 1. In the x-direction the minimum grid spacing was 0.11 R E (∼700 km) and the maximum grid spacing was 3.84 R E in the distant tail near a simulation outflow boundary. The grid spacing in the y-direction had a minimum of 0.13 R E (∼828 km) and a maximum of 2.2 R E . The grid spacing in the z-direction had a minimum of 0.13 R E (∼828 km) and a maximum of 2.48 R E because there were 24 fewer grid points in the y-direction. The regions of high grid resolution were broad and encompassed the near Earth and midmagnetotail.
We first imposed a steady southward IMF with a magnitude of 5 nT at the upstream simulation boundary for more than 4 h. The solar wind number density was 20 cm −3 , the thermal pressure was 20 pPa (temperature of 15.6 eV), and the velocity was 500 km s −1 in the x direction. These conditions ensured that there was no turbulence in the solar wind and all fluctuations were of internal origin. Thus, we only addressed internally generated turbulence.
A typical value for the solar wind temperature is around 10 eV (Kivelson and Russell, 1995), thus our solar wind temperature is higher than is typical. Because the thermal energy density of the solar wind is so much less than the kinetic energy density, the higher temperature should not make a significant difference. The density and temperature were set to ensure that the region of interest (plasma sheet) was confined to the maximum grid resolution area (Fig. 1). After the four hours of southward IMF, the IMF direction was changed to northward with the same magnitude (5 nT). The IMF was kept at this level for 6 h. In the paper, we will refer to the time of the northward turning as T N . Times after T N (240 min into the run) will be referred to by comparison with T N ; that is, five hours after the start of the run will be T N + 60 min. Our MHD simulations are almost always initiated with a southward IMF condition for one hour at least to remove features due the initialization process. A purely northward condition at the start dissipates these features less efficiently.

Overview of simulation results
The simulation results showed a complex pattern of flows in the magnetotail (Fig. 2), as can be seen in snapshots taken at 60 and T N + 180 min after the shift to northward IMF. The top two panels of this figure show color contours of the northsouth component of the magnetic field at the center of the current sheet (z = 0) overlain with unit vectors of the flow in the plane. The white isocontours give the locations of the last closed field lines. These isocontours were determined by tracing field lines from each MHD grid point in the equatorial plane and plotting isocontours at the boundary of the region of closed field lines. At T N + 60 min of northward IMF, there is still enough residual reconnection tailward of 40 R E to create a complicated open-closed boundary (small white islands). Although the solar wind and IMF were uniform the tail was highly variable and the B z distribution shows numerous meso-scale structures on the order of few Earth radii. We identified reconnection from the existence of strong velocity reversals and outflow jets, as well as B z ≈ 0. In addition, reconnection is often taking place at the open-closed field line boundary (Fig. 2). Finally, we can use the instantaneous value of the magnitude of the electric field due to resistivity (not shown). When the IMF turned northward, reconnection continued in the magnetotail though at a diminishing rate and the open-closed field line boundary became more orderly. Away from the boundaries, there are roughly x-aligned filaments where B z was large and positive. In the bottom two panels part of the plane, at T N + 240 min, is expanded in two stages to reveal nested vortices. Vortices and boundary waves are visible on the dawn side in the expanded regions. Similar vortices exist on the dusk flank.

Turbulent power spectra
Theoretical analyses indicate that a turbulent state features three ranges in the power spectral density (PSD), a driving (or energy-containing) scale at the longest wavelengths (timescales) at which energy is supplied to turbulent eddies from the large-scale flows, an inertial range at intermediate wavelengths, and a dissipative scale at still shorter wavelengths. We note that in a magnetofluid, there may be more ranges organized by the ion and electron scales as described in, for example, Sahraoui et al. (2010). Within the inertial and dissipative ranges the turbulent PSDs should be power law distributions with a steeper slope in the latter. To investigate the nature of the fluctuations during northward IMF in the magnetotail we calculated PSDs at different locations (Fig. 3). PSDs and other features of the turbulent spectra during southward IMF were shown in El-Alaoui et al. (2010). Figure 3 shows PSDs of the north-south magnetic field component taken from T N + 30 to T N + 360 min, all at z = 0 and y = −6R E , with x values of −12 R E , −16 R E , −20 R E and −24 R E . The three expected ranges for the slopes can be seen, in addition, at high frequencies, a flat part of the spectrum appears as the MHD time-step size is approached. It is probably related to the one second re-sampling from the variable MHD time steps (in roughly the 0.01 to 0.1 s time range). The break points between the three ranges were similar for the spectra at all locations in the plasma sheet at which they were computed. After inspecting the PSDs, we took the inertial range to be from 2.5 to 30 mHz and the dissipative range from 30 mHz to 200 mHz for this analysis. Least squares fits were performed to obtain the PSD slopes for the inertial and dissipative ranges. The predicted slope (dashed lines) for the inertial range of −5/3 (Kolmogorov, 1941) is added for comparison. We also display a line with a slope of −3 for comparison in the dissipative range, which is similar to the slope seen in the solar wind (e.g., Alexandrova et al., 2009;Chen et al., 2010;Sahraoui et al., 2010) and in the magnetosphere (Nykyri et al., 2006a;Vörös et al., 2005;Matthaeus et al., 2008).
To display the overall distribution of the PSD slopes, we made histograms of the PSD slopes at different points in the magnetotail z = 0 plane for the inertial and dissipative ranges. This was done for the southward IMF conditions from 60 to 240 min from the start of the simulation and for northward IMF from T N + 30 to T N + 360 min. These PSD slopes were taken at a rectangular grid of points. We used 337 locations in the equatorial magnetotail at z = 0 for −50 R E < x < −10 R E and −16 R E < y < 16 R E . Within these intervals we evaluated the slopes on a 17 by 19 point grid with 20 points near or outside the magnetopause removed. The grid intervals in x and y were 2 R E . The These are for northward IMF in the upper panels and for southward IMF in the lower panels. The left hand panels show the inertial range and the right hand panels show the dissipative range. Magenta lines have been drawn that corresponds to the edge of a region on high average resistivity on the dusk side (see Fig. 6).
resulting collection of PSD slopes, divided between the northward and southward intervals and inertial and dissipative ranges, were used to construct histograms. These histograms are shown in Fig. 4. The upper panels are for northern IMF while the lower panels are for southward IMF. The inertial range histograms are shown on the left and the dissipative range histograms are shown on the right. For northward IMF we found that the PSDs have slopes that were almost all between −2.4 to −1.3 in the inertial range with a median value of −1.90 and a standard deviation of 0.30. During southward IMF the PSDs have slopes mainly between −2.3 to −1.3 in the inertial range with a median value of −1.76 and a standard deviation of 0.27. The inertial range values overall are for both IMF directions and they are similar to the theoretical value of −5/3 (−1.67). It is likely that this is because the behavior in the inertial range is the result of basic fluid mechanical considerations resulting from the nonlinear term in the momentum equation in both MHD and ordinary fluids (Kolmogorov, 1941). As the standard deviations imply, the peak was sharper for southward IMF, possibly as a result of the different sources driving the turbulence. The inertial range results were consistent with the results of Weygand et al. (2005) who found a peak in the distribution of slopes at a spectral index of −2.0 and a weak secondary peak at a slope of −1.6. Weygand et al. (2005) did not distinguish different solar wind/IMF conditions. While the southward IMF and northward IMF PSD slope distributions were similar in the inertial range, the slope distributions in the dissipative range were different in the two cases. In the northward case the slopes are between −2 to −5 with a median value of −3.73 and a standard deviation of 0.76. For the southward IMF the primary range was −3.4 to −5.5 with some points reaching to about −2. The median value is −4.3 and the standard deviation is 0.75. As these numbers imply, the northward IMF has a higher median in the dissipative range. This may be due to the typically thinner current sheet in the southward IMF case that is associated with a more intense resistive electric field that leads to greater dissipation.
To illustrate the spatial distribution of the PSD slopes, Fig. 5 shows a color coded plot of the PSD slopes at the same locations in the equatorial magnetotail as in Fig. 4. Some points are added that were not included in the histograms. These points are located near and outside the magnetopause and are on the sunward edges of the plots on the dusk and dawn sides. These points can be most clearly distinguished by their shallow slopes (red and yellow) in the dissipative range for northward IMF. They are located in the top left and Nonlin. Processes Geophys., 19, 165-175, 2012 www.nonlin-processes-geophys.net/19/165/2012/ bottom left corners of the plot. During northward IMF, there is a trend toward more negative slopes on the dawn side of the tail than on the dusk side (above the magenta line) in the inertial range. In the dissipative range the pattern is reversed and the steeper slopes are on the dusk side. In the inertial range for southward IMF, the slopes are fairly uniform in the midtail region, but there is a trend toward more negative slopes for more tailward locations in the dissipative range. The level of turbulence in the inertial range throughout the tail is similar, on average, even though there is more dissipation in the distant tail. This suggests that the turbulent fluctuations in the inertial range are transporting energy tailward.
To relate the PSD slopes to the explicit resistivity in the MHD simulations we averaged the magnitude of this quantity in the equatorial plane over the northward and southward IMF intervals (Fig. 6). Northward IMF is shown in the upper panel and southward IMF in the lower panel. For north-ward IMF the overall level is lower than for southward IMF. For southward IMF the higher values are roughly consistent with the locations of thin current sheets and reconnection in the simulation. Snapshots of resistivity (not shown) indicate that it is high at reconnection sites where the current density is high. There are also regions of high dissipation near the magnetopause, in particular on the dusk side. The northward IMF result showed a diffuse region of relatively high resistivity skewed toward the dusk side creating dawn-dusk asymmetry in Fig. 6. To highlight the correspondence between the resistive dissipation pattern and the PSD slopes, we plotted a line (magenta) along the edge of the higher resistivity region on the dusk side in Fig. 6 like the ones in Fig. 5. The correspondence seems quite clear. Though we have not included the viscous interaction in the analysis, the PSD slope patterns seem to be consistent with the resistivity patterns.

Discussion and conclusions
To examine plasma sheet turbulence we used an MHD simulation of the magnetosphere. This simulation used idealized purely southward and northward IMF conditions to eliminate the effect of solar wind and IMF variations. To compare the results with observationally based studies of turbulence we applied similar analyses to the simulation.
We found that the fluctuations in the simulation had the characteristic properties of turbulence. We calculated PSDs in the central plasma sheet and found not only that the results were consistent with observations, but also that they had interesting physical implications. The MHD simulation exhibited nested vortices on multiple scales, with the largest scales associated with reconnection outflows and the diversion of high speed flows in the near-Earth region for southward IMF (El-Alaoui et al., 2010). During northward IMF, there are vortices near the flank and waves on the boundary that may be Kelvin-Helmholtz waves. The significant duskdawn asymmetries may reflect the influence of the ionosphere. A flow vortex generates a field aligned current (Vasyliunas, 1970) which will be, to lowest order, of opposite sign for the dusk and dawn sides, changing the ionospheric conductivity and presumably putting a drag on the flow that is asymmetric. The vortices in the tail show the structure expected of turbulence -driving range, inertial range and dissipative range. It is clear from the patterns of resistivity, flows and magnetic field that localized reconnection is occurring during southward IMF as we have seen in event studies (Ashour-Abdalla et al., 1999;El-Alaoui et al., 2009). The importance of localized reconnection is also supported by the PSD slope distributions. The northward IMF interval has turbulent vortices that are superficially similar to southward IMF but the PSD analyses and resistivity pattern reveal important differences. A further analytical method to determine whether fluctuations actually represent turbulence is to compute Probability Distribution Functions (PDFs). El-Alaoui et al. (2010) found that the PDFs from the MHD simulation had similar properties to those computed from observations (Weygand et al., 2005(Weygand et al., , 2006. Intermittency in magnetic turbulence can be inferred (Castaing et al., 1990;Sorriso-Valvo et al., 1999) from a non-Gaussian scaling of PDFs of the magnetic field fluctuations over time periods corresponding to frequencies that fall near the transition from the dissipative to the inertial range. We computed PDFs for both northward and southward IMF and for both B z and v x and found expected characteristics of intermittent turbulence (Fig. 7). Intermittent turbulence corresponds to localization of dissipation (Biskamp, 1993), which is exactly what is observed in the distribution of resistivity in the simulation.
We used the same technique as Weygand et al. (2005) to construct PDFs of the magnetic field fluctuations. We apply the function: where B i (t) is a component of the magnetic field at time t and B i (t +τ ) its value at a fixed time τ later (Sorriso-Valvo et al., 1999). First, a running average of 1000 s was removed from the magnetic field time series in order to eliminate longer period fluctuations in the magnetic field. As discussed by Weygand et al. (2005) the values of τ corresponds to spatial lengths determined by the plasma speed. Small values of τ (τ ≤ 10 s) correspond to distances within the dissipative range for the simulation and large values of τ (τ ≥ 900 s) are associated with physical dimensions on the order of the driving scale of the plasma. Intermediate values of τ correspond to the inertial range. The top two rows of Fig. 7 show results from the northward IMF interval and the bottom two rows show southward IMF. Shown are PDFs of the z component of the magnetic field (second and fourth rows) and the x component of the flow (first and third rows) in the simulation calculated by using (1) for four different values of τ . The dashed curves show Gaussian distributions for comparison that were fit by eye to the cores of the PDFs at x = −16 R E and y = 6 R E and z = 0 R E . If the relationship between a measurement at given time is random with respect to one at later time the result is a purely Gaussian distribution. The PDFs show significant departures from a Gaussian distribution with large deviations from the average that can be described as "wings". The PDFs for the northward and southward cases are similar. For τ = 900 s the PDFs of the fluctuations in B z are the most similar to Gaussian distributions. The PDFs for v x have prominent "wings". These "wings" in the v x PDFs are probably due to intermittent high-speed flows related to reconnection (Borovsky et al., 1997). The southward IMF inertial range and the northward IMF inertial range for points on the dusk side are generally consistent with the theoretical value (Kolmogorov, 1941) while the dawn side for northward IMF has a steeper slope, possibly related to the ionospheric dissipation of the vortices on the dawn flank. For southward IMF less asymmetry is expected because the localized reconnection that drives the turbulence happens at "random" locations. We note that there are steeper (more negative) slopes in the more distant tail for southward IMF, and the average resistivity is much higher in some parts of the magnetotail than in others. Nevertheless, the slopes in the inertial range do not show any clear pattern. The frequencies of the breaks (between the inertial and dissipative ranges) are also fairly uniform. If the dissipation, as seen in both the PSD slopes and the average resistivity, is taking place in local regions, it might be expected that this would steepen the slopes in the inertial range where there is more dissipation. That this is not seen suggests that the energy of the turbulent eddies are being effectively transported across the magnetotail. For northward IMF, the overall pattern of the inertial range slopes is anticorrelated with those in the dissipative range. This cannot be easily explained by transport. Some other factor, perhaps drag from the ionosphere, which has an asymmetric conductivity distribution, may be at work. In both the northward and southward cases the standard deviations are much lower in the inertial range than in the dissipative range. We know that the inertial range slopes can be derived from basic scaling arguments (Kolmogorov, 1941) and it probably takes a large effect to change the slope significantly. Dissipation, though, varies widely. The turbulent plasma sheet system may be responding to variations in dissipation by the transport of energy to the dissipative regions.