Numerical considerations in simulating the global magnetosphere

Magnetohydrodynamic (MHD) models of the global magnetosphere are very good research tools for investigating the topology and dynamics of the near-Earth space environment. While these models have obvious limitations in regions that are not well described by the MHD equations, they can typically be used (or are used) to investigate the majority of magnetosphere. Often, a secondary consideration is overlooked by researchers when utilizing global models – the effects of solving the MHD equations on a grid, instead of analytically. Any discretization unavoidably introduces numerical artifacts that affect the solution to various degrees. This paper investigates some of the consequences of the numerical schemes and grids that are used to solve the MHD equations in the global magnetosphere. Specifically, the University of Michigan’s MHD code is used to investigate the role of grid resolution, numerical schemes, limiters, inner magnetospheric density boundary conditions, and the artificial lowering of the speed of light on the strength of the ionospheric cross polar cap potential and the build up of the ring current in the inner magnetosphere. It is concluded that even with a very good solver and the highest affordable grid resolution, the inner magnetosphere is not grid converged. Artificially reducing the speed of light reduces the numerical diffusion that helps to achieve better agreement with data. It is further concluded that many numerical effects work nonlinearly to complicate the interpretation of the physics within the magnetosphere, and so simulation results should be scrutinized very carefully before a physical interpretation of the results is made. Our conclusions are not limited to the Michigan MHD code, but apply to all MHD models due to the limitations of computational resources.


Introduction
Global magnetohydrodynamic (MHD) models are extremely useful for studying the large-scale dynamics of magnetospheric systems.There are more simple models that empirically describe some characteristics of the magnetosphere, for example, the location of the magnetopause (Shue et al., 1997;Shue et al., 1998), or the magnetic field within the magnetosphere (Tsyganenko, 1989(Tsyganenko, , 1995(Tsyganenko, , 2002)).These types of models are limited to the range of the data that was used to derive them.Other models, such as reconnection models (e.g., Shay et al., 2004) or models of the inner magnetosphere (Wolf et al., 1982;Sazykin et al., 2002) can be dominated by the boundary conditions, such that their lack of self-consistent feedback onto the magnetosphere may inhibit their usefulness in understanding how the real global system evolves in time.MHD codes can describe the self-consistent magnetic field and plasma interaction with the solar wind and transport of these through the global magnetosphere.In addition, with coupling to other models, such as those that describe the ionospheric electrodynamics, they can be used to study the self-consistent coupling between these complex regions (e.g., Fedder and Lyon, 1987;Raeder et al., 2001a,b;Ridley et al., 2004Ridley et al., , 2003;;Palmroth et al., 2004;Wiltberger et al., 2004).Recently, MHD models have been coupled to inner magnetospheric models to more realistically model the higher energy particles that are not accurately modeled by a single fluid (De Zeeuw et al., 2004;Toffoletto et al., 2004).Further, research is being conducted on methodologies to make reconnection within MHD codes more physically realistic, since ideal MHD does not actually allow reconnection (e.g., Kuznetsova et al., 2007).
One of the primary methods for examining the magnetosphere has been through the use of the ionospheric potential.If it is assumed that there is essentially no potential drop along magnetic field lines, the ionospheric potential maps out to the magnetosphere (or visa-versa).The potential then describes how plasma flows through the magnetosphere.
Published by Copernicus Publications on behalf of the European Geosciences Union.
Additionally, it is thought that the potential difference along the open-closed field-line boundary in the ionosphere represents the reconnection potential (e.g., Siscoe and Huang, 1985).Therefore, measuring the cross polar cap potential gives insight into the amount of dayside reconnection that may be occurring.
There are many different measurements of (or techniques for deriving) the high-latitude ionospheric electric field and plasma flow (e.g., Kelly, 1983;Richmond and Kamide, 1988;Ruohoniemi and Greenwald, 1996).For example, the Special Sensor for Ions, Electrons, and Scintillation (SSIES) instrument on board the Defense Meteorological Satellite Program (DMSP) satellites estimates the cross-satellite-track ion drift (Rich and Hairston, 1994).The drift in the upper ionosphere is primarily in the E × B direction, and so a measure of the cross-track drift essentially specifies the electric field in the satellite track, which can then be used to estimate the ionospheric potential (e.g., Hairston et al., 2003).Many model validation studies have been conducted comparing DMSP measured cross-track drifts with ionospheric drifts produced by models (e.g., Fedder et al., 1998;Raeder et al., 1998;Ridley et al., 2002;Bekerat et al., 2005;Kihn et al., 2006), and the DMSP measurements are used as the official metric put forth by the National Space Weather Implementation Plan to validate magnetospheric codes.
One problem with using the ionospheric potential as a metric for determining how accurately the MHD codes are modeling the magnetosphere is that the ionospheric potential is not a direct product of the MHD equations.As described by Fedder and Lyon (1987); Goodman (1995) and Ridley et al. (2004), the ionospheric potential is derived from the field-aligned currents computed from the MHD solution and an ionospheric conductance pattern, which has a great deal of flexibility in its specification.Fedder and Lyon (1987) and Ridley et al. (2004) showed that the ionospheric potential strongly depends on the conductance pattern.Both the strength and shape of the potential can be strongly controlled by the specified conductance pattern.Therefore, the conductance is an intermediate variable that can be tuned, such that, even if the field-aligned currents output from the MHD code are incorrect, the potential can still more or less match the empirical potential values (e.g., Ridley et al., 2002).One justification for doing this is that the ionospheric potential is of primary importance for the ring current and plasmasphere (e.g., Carpenter et al., 1972;Sazykin et al., 2002;Liemohn et al., 2002Liemohn et al., , 2005)).Therefore, with some tuning of the ionospheric conductances, such that the potential pattern is more accurate, ring current injections and plasmaspheric erosion events can be modeled better.In addition, if the plasmasphere and the ring current are specified more precisely, then the wave-driven energization of the radiation belts, which takes place where the ring current and plasmasphere overlap, may be modeled better.
As described by Goodman (1995) and Ridley et al. (2004), the field-aligned currents within the ionosphere are derived by taking the curl of the magnetic field at some radius slightly away from the inner boundary of the MHD code.These fieldaligned currents are mapped down to the ionosphere along dipole field-lines, and are assumed to focus as the magnetic field becomes stronger.Because the FAC strength scales with the magnetic field magnitude, in order to capture all of the FAC that is generated at a distant location, the resolution needs to scale as the magnetic field also.For example, if the source region is located at 10 R E and it contains current systems with a characteristic size of 1 R E , then the region in which the ionospheric current is drawn from (at 3 R E ) needs to have a grid resolution that can resolve 1/9 R E sized features, due to the magnetic field converging towards the Earth.Even with a grid resolution of 1/32 R E the field-aligned currents generated at 10 R E will start to diffuse.Running with 1/32 R E resolution or better in the inner magnetosphere is impractical, due to both the number of cells that would result, and the incredibly small time-step that would be needed to accurately model the propagation of information across these very small cells (most codes use explicit time-stepping).For example, the explicit time-step for a 1/4 R E cell near the inner boundary may be on the order of 0.01 s.With 1/32 R E extending out to 4 R E , there would be roughly 8 million more cells, and the time-step would be reduced to 0.00125 s.Considering that many magnetospheric simulations have less than 2 million cells, the run time would increase by a factor of 40 by going to 1/32 R E even in such a small region.It should also be noted that the MHD equations break down when the resolution becomes smaller than the local ion Larmor radius.
While the above spatial resolution argument is true, it is still possible to obtain simulation results that agree surprisingly well with observations.The reason is, most likely, that the MHD model gets the global features and overall current systems correctly, and the numerical errors on the smaller case are compensated by the tuning of the magnetosphereionosphere coupling and the ionosphere model itself.
Despite the reasonable agreement with observations, the effects of the numerical schemes remain significant, and it is important to be aware of the effects and limitations of numerical schemes.This is the subject of our paper.We will show, for example, that since first order and second order numerical schemes converge at different rates, they give different results at the practically available grid resolutions.Also, the limiters that are used to ensure positivity (of density, energy, and pressure) in second and higher order schemes play an important role.Robust limiters strongly limit the slope to very small values, ensuring positivity, but add more numerical diffusion, while more accurate limiters allow larger slopes (and slope changes), but may allow negative values to form due to the non-linearity of the equations.
This study examines these types of numerical considerations, using the ionospheric cross polar cap potential (CPCP) as a diagnostic for how accurate the solutions are.While there is an empirical understanding of the expected value of the actual CPCP as a function of interplanetary magnetic field (e.g., Rich and Hairston, 1994;Weimer, 1996;Ruohoniemi and Greenwald, 1996;Boyle et al., 1997;Siscoe et al., 2002;Ridley, 2005), there is not a very good understanding of what ideal MHD (with the allowance for numerical resistivity in the reconnection sites) should give.Some codes tend to give a very large CPCP (e.g., Fedder et al., 1998), while other codes give a value that is smaller than expected (e.g., Palmroth et al., 2005).In this paper, the values of the CPCP vary significantly as a function of many different parameters.This is meant to illustrate where some of the large differences between MHD codes may originate.The study then explores the 4 May 1998 storm, showing how the lessons learned from idealized simulations apply quite nonlinearly to a storm, making the interpretation quite difficult.

Modeling the magnetosphere
There are several groups that have developed MHD models of the global magnetosphere.The model of Ogino and Walker (1984); Ogino et al. (1985) was one of the first global models of the magnetosphere.It has also been applied to other planetary objects, such as Jupiter (e.g., Ogino et al., 1998;Walker and Ogino, 2003) The Lyon-Fedder-Mobarry (LFM) MHD code (Fedder and Lyon, 1987;Fedder et al., 1998;Lyon et al., 2004) works on a distorted spherical mesh in order to resolve the important parts of the domain with higher resolution.The LFM has recently been coupled with a thermosphere-ionosphere model (Wiltberger et al., 2004;Wang et al., 2004) and the Rice Convection Model (Toffoletto et al., 2004).The Open Geospace General Circulation Model (Open GGCM) is also an MHD-based code that has been used in many scientific investigations (Raeder et al., 1996(Raeder et al., , 1997(Raeder et al., , 1998(Raeder et al., , 2001a)).The Open GGCM uses a stretched Cartesian grid, concentrating high resolution in areas of interest in the magnetosphere.This has been used to successfully model Flux Transfer Events, for example.Robert Winglee's code solves the multi-fluid MHD equations (Winglee and Elsen, 1995;Winglee, 1998), having the ability to resolve oxygen, hydrogen, and helium ions in the magnetosphere.In addition, the code has been utilized at other planetary objects (e.g., Snowden et al., 2007;Harnett and Winglee, 2007).The Integrated Space Weather Prediction Model MHD code is similar to Winglee's code, but goes a step further -it models the magnetosphere and ionosphere as one system, using a spherical grid in the inner magnetosphere and ionosphere and a Cartesian grid in the outer magnetosphere (White et al., 1998;Sonnerup et al., 2001).The first-order Grand Unified Magnetosphere Ionosphere Coupling Simulation (GUMICS-4) code is similar to BATSRUS (described below) in some ways -namely it has a cell-based adaptive grid structure (while BATSRUS is block-based) and it uses a similar numerical scheme (e.g., Janhunen, 1996;Palmroth et al., 2001Palmroth et al., , 2003Palmroth et al., , 2005, and references within), and references within).The first-order nature of GUMICS-4 is compensated for by the ability to have higher resolutions in regions in which there are large gradients.The MHD code by Tanaka (1995) is unique in that it utilizes an unstructured grid to model the 3-D magnetosphere-ionosphere system.
The University of Michigan's MHD code is the Block Adaptive Tree Solar-wind Roe-type Upwind Scheme (BAT-SRUS) (Powell et al., 1999b;Gombosi et al., 2002;Ridley et al., 2002;Gombosi et al., 2004) that can model the Earth's global magnetosphere (Kabin et al., 2003(Kabin et al., , 2004;;Vogt et al., 2004;Ridley et al., 2006;Watanabe et al., 2005;Ridley, 2007;Fairfield et al., 2007), as well as the solar corona (Manchester et al., 2004b), the inner and outer heliosphere (Opher et al., 2003;Manchester et al., 2004c,a), Mercury (Kabin et al., 2000), Venus, Mars (Ma et al., 2004), Jupiter (Kabin et al., 2001), Saturn (Hansen et al., 2000(Hansen et al., , 2005)), Uranus (Tóth et al., 2004), Titan (Ma et al., 2009), and comets (Gombosi et al., 1999).Because BATSRUS is a general MHD solver, and is not explicitly designed for solving a specific problem, it has many more options than typical magnetospheric MHD codes.Here we describe a few of the settings that BATSRUS has that we will specifically explore within this study: -Grid resolution.BATSRUS has an adaptive mesh that can resolve regions of interest without increasing the total number of grid cells to impractical values.For block adaptive grids the limited reconstruction procedure near the resolution change interface is handled in the way described by Sokolov et al. (2006).The adaptive grid can evolve during the simulation, although this feature is rarely used in the magnetosphere simulations, due to the fact that overall structure of the magnetosphere is approximately stationary, and the regions of interest (such as the bow shock, the magnetopause, the inner magnetosphere and the magnetotail) remain in a confined enough region that high resolution can be utilized to cover almost all cases that may exist (e.g., using 1/8-1/4 R E in a large block covering 4 < X < 16 R E would include the magnetopause for almost any event).The adaptive mesh can be still very useful for examining small-scale features such as Flux Transfer Events and other propagating structures.BATSRUS is typically run for the magnetosphere with an unstructured, but static, grid.The typical grid is fine near the inner boundary at 2.5 R E , the magnetopause around 10 R E on the dayside, and in the central tail.The rest of the magnetosphere may be solved at 1/2 R E resolution, except in the mid-to deep-tail, which may be solved with 1-2 R E cells.This may take something on the order of one million grid points, and may be solved for in approximately realtime on 128 processors on current super-computers.
The run time of the code is strongly controlled by the number of cells included in the solution, the size of the smallest grid cell, and the number of processors used for the simulation.solving the inner magnetosphere with 1/2 R E resolution within −10 < X,Y,Z < 10 R E , this part of the simulation would include 64 000 grid points, and would take a certain run-time of N (on a given number of processors).
If the resolution was increased to 1/8 R E cells in this region, the number of cells would jump to 4 096 000, and the time-step would have to be reduced by a factor of four, due to the cells reducing in size by a factor of four (in each direction).The code would run 256 times slower (i.e., a run-time of N × 256 on the same number of processors), simply because of the addition of cells and the shrinking of those cells limiting the time-step.Therefore the grid size is an extremely important consideration when conducting a simulation.
-Solvers.There are many different ways to discretize the MHD equations and solve them.BATSRUS utilizes finite volume schemes in which the fluxes of the primary state variables are determined at the cell faces and are added to and subtracted from the volume averaged states.Finite volume schemes conserve mass, momentum and energy, and produce correct solutions for problems involving discontinuities, like the bow shock.The fluxes at the faces are constructed with upwind type schemes that use information from the upwind direction in terms of the wave propagation.The simplest, fastest, most robust but also most diffusive is the Rusanov flux function (Rusanov, 1961), which utilizes the fastest wave speed to determine the numerical diffusion.This property is in contrast with the Roe scheme (Roe, 1981;Sokolov et al., 2008), in which the amplitude of numerical diffusion is minimum admissible (hence, different) for each eigenwave.No excessive diffusion is applied, but the low diffusivity and high accuracy of the Roe scheme is achieved at the cost of lower computational efficiency (left and right of eigenvectors of the Roe matrix are to be constructed in conjunction with the numerical diffusion for each of the eigenmodes) and robustness.The Linde flux function (Linde et al., 1998) and the Sokolov solver (Sokolov et al., 1999) are considered to be schemes of the HLL (Harten-Laxvan Leer) family and lie between the Roe and Rusanov schemes.On one hand, the coefficient of numerical diffusion in these schemes is the same for any mode of perturbation (similar to the Rusanov scheme).Although the numerical diffusion of the HLL schemes is somewhat excessive in a general case, it is lower than that for the Rusanov scheme in some cases.Specifically, for supersonic flow, the HLL schemes all reduce to the purely upwind scheme (the Courant-Isakson-Friedrichs scheme), which is known to have the lowest numerical diffusion for this particular class of motions.Summarizing, the Roe scheme is always optimal, the HLL schemes are optimal in some limiting cases, while the non-selective Rusanov scheme is never optimal.
-Limiters.Second order total variation diminishing (TVD) schemes employ non-linear slope limiters to interpolate the cell centered values to the faces.The slopelimiters play a crucial role in avoiding spurious oscillations near discontinuities while still getting accurate solutions in smooth regions.BATSRUS contains various limiters, including Koren's third-order (in smooth monotone regions) limiter (Koren, 1993), the monotonized central (MC) limiter, and the minmod limiter.
In this study we use the MC limiter with a β parameter that can vary between 1 and 2. When β = 1 the MC limiter coincides with the minmod limiter, which is the most diffusive and robust of all TVD limiters.On the other hand β = 2 corresponds to the sharpest version of the MC limiter that is much less diffusive, but also less roubust.For magnetosphere simulations BATSRUS is often run utilizing the MC limiter with β = 1.2.
-Speed of light.When the semi-relativistic MHD equations are solved for (Boris, 1970;Gombosi et al., 2002), the fast mode wave speed is limited to the speed of light, which can be artificially reduced, thereby allowing the code to take larger time-steps in explicit time integration schemes (where the time-step is limited by the cell size divided by the maximum wave propagation speed), significantly reducing the run-time of the simulation.Typical reduction factors range from 0.02 to 0.005 (corresponding to c = 6000 km/s to c = 1500 km/s).This reduction is often referred to as the "Boris correction", as we will refer to it throughout this study.The fastest waves in the (simulated) inner magnetosphere typically reach speeds of 30 000 km/s.By reducing the speed of light in a simulation in this way, the run-time of a simulation can be reduced by a factor of 5 to 20, respectfully.The lower the reduction in the speed of light, the less stable the code is and the less physical the solution can be.For example, during the 29 October 2003 Halloween storm, the solar wind velocity reached speeds of 2000 km/s, so artificially reducing the speed of light to 1500 km/s would not be very physical.
In addition to the run-time speed-up, the artificially reduced fast mode wave speed reduces the numerical diffusion in the simulation.This is because the numerical diffusion in TVD solvers is directly proportional to the local wave speed.Therefore, when the speed of light is reduced and the wave speeds are then reduced, the diffusion can be diminished by a significant amount.
-Implicit versus explicit time-stepping.Another method of allowing the code to run in a manageable amount of time is to use an implicit time-stepping scheme instead of a typically-used explicit scheme (e.g., Tóth et al., 2005b).The implicit scheme involves an iterative linear solver, and a single implicit time step requires many iterations.In a typical magnetospheric storm simulation the implicit time-stepping can be run with a time step of 2-5 s, which is three orders of magnitude larger than the time steps of a typical explicit simulation without the Boris correction.The net practical gain is an approximate 50 times speedup, which is substantial.When the Boris correction is used, the net speed up is less.One also needs to pay attention to the fact that for fast moving features/waves/shocks, the implicit code has more numerical diffusion than the explicit code.In general a feature of interest should move less than a third of a grid cell per implicit time step (Tóth et al., 2005b) to maintain the same quality of accuracy as in the explicit code.
-Inner boundary density.Currently, most magnetospheric simulations using BATSRUS use a constant mass density at the inner boundary (28 AMU/cm 3 ).This density diffuses away from the boundary (due to numerical diffusion) and is ejected into the magnetosphere due to gradients in pressure (similar to Winglee, 2000).This is an extremely low density at 2.5 R E at low latitudes where the plasmasphere should be, and is a relatively high density for the polar field-lines.Global MHD codes do not have the grid spacing to resolve the possible orders of magnitude change in mass density through the plasmapause boundary.Therefore, a lower density is chosen, so the density is better represented at distances of 5-6 R E in the equatorial plane.At high latitudes, there are a significant number of studies currently being conducted on the actual density, field-aligned flow velocities and/or particle fluxes at high altitudes (Gombosi and Nagy, 1989;Abe et al., 1993;Strangeway et al., 2000Strangeway et al., , 2005;;Moore et al., 2007;Yao et al., 2008;Nosé et al., 2009).Recently, Glocer et al. (2009) implemented a coupling between BATSRUS and a Polar Wind Outflow Model that provided a more-realistic representation of the inner magnetospheric density and structure.It is clear that there is significant structure in the high-latitudes, and that a constant value is a dramatic simplification.The results presented here simply show how changes in the density at a uniform inner boundary could possibly affect the global magnetospheric state and the numerics.
In this study, we will show how each of the options listed above can alter the solution in an idealized simulation and in a realistic simulation of a storm.
In the idealized simulation the BATSRUS code is coupled with the ionosphere model, described by Ridley et al. (2004).The ionospheric conductance has a dayside solar zenith angle dependence and an auroral component that is strongly dependent upon the field-aligned currents.In the storm simulations, BATSRUS is also coupled to the Rice Convection Model (RCM) (Wolf et al., 1982;Sazykin et al., 2002;De Zeeuw et al., 2004) of the inner magnetosphere.The coupled models run as part of the Space Weather Modeling Framework (SWMF) (Tóth et al., 2005a).

Results
In numerical models, the grid resolution matters a great deal.Discontinuities require a certain number of grid cells to be resolved.For example, at the magnetopause it takes about 3-4 cells to capture the strong velocity and magnetic field changes.If each cell were 2 R E in size the magnetopause would be modeled as 6-8 R E wide, so it would be completely smeared out and indistinguisable from the bow shock and the inner magnetosphere on the day side.Conversely, if the cell size were 1/4 R E , the magnetopause would be modeled as being 0.75-1 R E in width, still large, but reasonably separated from other features of the magnetosphere.
Figure 1 shows an example of the magnetospheric modeling results with the block structure superposed.The simulations were run for 35 000 iterations using the local timestepping method, in which the code was pushed towards a steady-state solution, as described by Ridley et al. (2002), and references within.While the solution may not be perfectly steady, as described below, the results are good enough to illustrate the differences in the numerics.While the code was running, the resolution was dynamically altered.The initial 4000 iterations were conducted with the majority of the magnetosphere being resolved with 1 R E grid cells (extremely low resolution).The resolution was changed such that the magnetosphere was resolved to 1/2 R E , and the solution evolved for 4000 more iterations.The resolution in a sphere around the inner boundary was changed to 1/4 R E , and the solution evolved for 5000 more iterations.Next, this sphere was refined to 1/8 R E and the solution was allowed to converge for 10 000 iterations.Finally, the sphere was refined to 1/16 R E and the code was run for 12 000 more iterations.With 1/16 R E resolution, the code never reached a true steady-state, but showed signs of a slightly oscillatory solution.This behavior was not observed at lower grid resolutions.
The solar wind boundary conditions at the upstream face were held constant with the IMF B x = B y = 0, B z = −10 nT, and the solar wind number density at 5 cm −3 and velocity equal to 400 km/s in the −X direction (i.e.purely antisunward).For reference, the ionospheric cross polar cap potential given by the Weimer (2001) model for the above solar wind conditions is 130 kV.
Figure 2 shows the ionospheric field-aligned currents (FACs) and potential patterns derived self-consistently with the simulation described in Fig. 1.A FAC-dependent conductance pattern (Ridley et al., 2004) is utilized to make the simulations as close to what would be used for a real event study as possible.As the grid becomes more refined, the strength of the field-aligned currents change significantly.This drives the ionospheric potential pattern to become stronger as well.While the strength of the potential and FAC pattern changes significantly as the run progresses, the shape changes very little, indicating that the shape of the source region has not changed, but that the total amount of FAC making it to the inner boundary of the MHD code is increasing with increasing resolution.Fig. 3.The top panel shows the ionospheric cross polar cap potential as a function of run-time, with first (dashed) and second (solid) order solvers.The smallest grid cell size (in R E ) is indicated at the top of the panels.The grid resolution changes are indicated by the vertical dotted lines.The bottom panel shows the ratio between the second and first order cross polar cap potentials as a function of runtime (solid line), and the ratio between the maximum CPCP in the second order scheme and the rest of the CPCPs in the second order scheme simulation.The numbers above the line indicate the value in the last iteration before the end of the particular resolution.The numbers below the line are the differences between the ratios (i.e., the numbers above the line) in the current and previous resolutions.
The observed behavior is quite expected, because of the size of the grid cells with respect to the size of the fieldaligned currents.For example, if the size of the narrow fieldaligned current in the ionosphere is 1 • in latitude, it would be approximately 1/2 R E in size at 3 R E , simply due to the geometry of the field-lines.If four grid cells were needed to resolve a current structure such as this (which is optimistic, since the curl of the magnetic field, B, or ∇ × B, has to be resolved), the resolution needs to be a minimum of 1/16 R E .At coarser resolutions, the code can not numerically resolve features as narrow as this.

First versus second order schemes
A closer analysis of the change of the CPCP is available in Fig. 3.This figure shows how the CPCP changes as a function of iteration (and resolution) for two simulations -one with a first-order scheme and one with a secondorder scheme.The CPCP converges to a constant value for each resolution, except for the 1/16 R E resolution using the second-order scheme.The second-order scheme results in a much higher CPCP than the first-order scheme, and the ratio between them depends on the resolution, as shown in the bottom plot.For the 1/16 R E resolution, the ratio between the second-and first-order schemes is almost 2.3, while, with 1 R E resolution, the ratio is only 1.5.This means that the second-order scheme is converging to its solution much faster than the first-order scheme, as one would expect.
The first and second order schemes are not approaching the same values.It is expected that, as the resolution becomes infinitely small, the solution will be exactly the same.The problem with this argument for the runs described above, is that the resolution in most of the region is remaining fixed, and only the region surrounding the body is being refined.This means that the source region of the field-aligned currents is not being refined, and it is expected that the solutions will be different between first and second order for relatively coarse resolutions.For a true convergence study, both the source region and the body would have to be refined at the same rate.This is difficult to do, since the number of cells would be prohibitive.
The convergence can be quantified by examining how the CPCP changes as the resolution changes.For example, if the code were completely grid-converged, the change between one resolution and a higher resolution would be zero.Any change between the resolutions indicates that the solution is not grid-converged.The dashed line in the bottom plot shows the ratio between the CPCP and the maximum CPCP during the simulation.The numbers above the lines indicate the ratio in the last iteration before the end of the particular resolution.The numbers below the line are the differences between the ratios in the current and previous resolutions.If the code was moving towards grid convergence, the numbers below the line would decrease as the resolution increased, indicating that the solution is changing, but the changes are decreasing in magnitude.
This figure shows that the differences in the secondorder scheme are actually increasing for all of the resolution changes except for the 1/8-1/16 R E jump, which indicates that the code is starting to converge only at the highest resolution.This means that the 2nd order MHD code, with a Rusanov scheme and a minmod limiter, is very far away from grid convergence in the inner magnetosphere and that far more resolution is needed before grid convergence is reached.The lack of grid convergence is responsible for many issues in attempting to model the inner magnetosphere using MHD.For example, at low resolutions, the potential is significantly underestimated (recall that the "true" potential is close to 130 kV), which is what is indicated in Fig. 3 (at, for example, 1/4 R E resolution).While this fact is used as a diagnosis here, it is quite important for modeling the particle motions in the inner magnetosphere.Additionally, diffusion in the code works on all of the MHD parameters.For example, the field-aligned currents are underestimated because the magnetic field is not resolved well enough.This same diffusion acts on the velocities, such that the potential applied at the inner body may diffuse partially, causing the potential to change along the field-line.This is obviously incorrect, since field-lines in the magnetospheric region should be very close to perfect conductors, except in very specific regions (such as the reconnection site and in the auroral gap region).
With a second-order scheme, utilizing the highest practical resolution, the values of the cross polar cap potential are approximately 160 kV, which is actually higher than what is expected, given that B z = −10 nT in the interplanetary magnetic field.Further, this plot is indicating that the "ideal" MHD solution (with numerical diffusion in the reconnection site), should produce a higher CPCP: it is clear from Fig. 3 that if the code were run with 1/32 R E resolution at the inner boundary, the CPCP would be higher than 160 kV.The reason for this overestimation of the CPCP is investigated below, when the dayside reconnection is resolved with higher resolution.

Numerical flux functions
Figure 4 shows the exact same results as in Fig. 3, but comparing the 2nd order Rusanov scheme and the 2nd order Linde scheme (Linde et al., 1998).The Linde scheme has less numerical diffusion than the Rusanov scheme when the flow velocity is significant compared to the fastest wave speed.Therefore one may hope that the Linde scheme will provide more accurate results than the Rusanov scheme.This is clearly not the case for the ionospheric CPCP, indicating that the Linde and Rusanov schemes are very similar for regions with small flow speeds, such as the high-latitude magnetosphere that forms the region-1 current system.
Figure 5 shows the same type of plot as Fig. 4, except it compares the Rusanov and Roe schemes (Roe, 1981).The Roe scheme is less diffusive than the Rusanov (and Linde) scheme in almost all possible circumstances.Indeed the plot confirms that the Roe scheme converges towards the final (highest resolution) CPCP value much more rapidly than the Rusanov scheme.In fact the Roe scheme may be able to achieve a grid-converged solution at a slightly higher resolution, while the Rusanov scheme does not show much sign of convergence.
Somewhat surprisingly, at the highest resolution the Roe scheme produces a smaller CPCP than the Rusanov scheme.This is a result of two mechanisms: the reduced numerical reconnection on the dayside, and a stronger pressure on the night side.Both tend to decrease the expected CPCP.The first will be discussed in more detail later, when the dayside resolution is changed.The second is much more complex, since the inner magnetospheric dynamics is driven by the ionospheric potential.As the potential becomes larger, because of the increasing resolution, the flow shears in the inner magnetosphere increase, causing a larger pressure peak to form just tailward of the Earth.In the Rusanov scheme, this pressure peak is weak and diffuse.In the Roe scheme, the pressure peak is stronger and less diffuse, allowing the formation of more significant region-2 currents.These region-2 currents close some of the Pedersen currents in the ionosphere, thereby reducing the cross polar cap potential, as described by De Zeeuw et al. ( 2004) and references within.
Figure 6 shows the pressure distribution and magnetic field topology in the magnetospheric X −Z plane for the two runs compared in Fig. 5 in the last iteration.Each plot has the same color scale of pressure, so it is clear that the run with the Roe solver has a much stronger pressure on the night side.This pressure causes a stretching of the magnetic field-lines in the tail, resulting in a much longer tail, which is more consistent with observations of the tail length (Sergeev et al., 2000;Nagai et al., 2001Nagai et al., , 2005;;Nagai, 2006;Slavin et al., 2003).The stretching of the field-lines also serves to allow the plasma to be adiabatically heated more as it is advected in from the tail, thereby increasing the pressure in the inner magnetosphere, resulting in a nonlinear feedback.
The increased pressure also drives more current into the ionosphere, leading to strong region-2 currents.Figure 6 also shows the FAC patterns produced by the two MHD simulations.The simulation driven with the Roe scheme shows the currents to be sharper than the Rusanov simulation.In addition, there are more significant region-2 currents on the night side in the Roe simulation as compared to the Rusanov simulation.
While the Roe scheme produces superior results, there are also some drawbacks.It is much more complicated than the Rusanov scheme, and therefore it is about three times slower to run.This alone is not a reason to favor the simpler schemes, since the simulation results obtained with the Roe scheme are better.A significant issue with the Roe scheme is that it cannot be used with the Boris correction due not being able to solve for the Eigen vectors of the semi-relativistic MHD equations.The ideal solution for the Earth's magnetosphere MHD modeling would appear to be to use the Roe scheme with the implicit time-stepping algorithm, as will be discussed in Sect.3.7.

Limiters
Figure 7 shows a comparison between simulations using the Rusanov second order scheme utilizing the minmod (same as the Rusanov scheme above) and the MC limiter with β = 1.2.It is quite clear that the choice of limiters is quite important for controlling the amount of diffusion that exists in the MHD code.The field-aligned currents in the ionosphere are much stronger and sharper with the more aggressive limiter.
One of the issues with using a sharper limiter is that the code becomes less stable.This is because steep gradients in the solution can cause oscillations within the solution.In regions in which there is a strong gradient from a high to a very low value in something that should be positive, the oscillations may cause the values to become negative.This typically happens in the lobes, where the thermal pressure is significantly smaller than the magnetic pressure.A common method for handling this type of code crash is to put a floor (i.e., a minimum value, such that if the solver produces a smaller value, it is replaced with the floor value) on the pressure or density, which allows the code to continue to run, but may produce spurious results.BATSRUS does not use a floor for the pressure of the density, so if these states become negative, the code crashes.
Much of the variation that occurs when the resolution is 1/16 R E in both the simulation with the MC limiter and the Roe solver, is caused by the tail actually being dynamically unstable: the code never reaches a true steady-state in these cases, and should be run in a time-accurate mode to examine plasmoids that are being launched during the simulation.This is much closer to what observations of the tail dynamics show (e.g., Slavin et al., 2003).

Artificial reduction of the speed of light
One of the main issues with running time-dependent MHD simulations of the magnetosphere of the Earth is that the time-step is very small, due to the extremely large Alfvén wave speed in the high-latitude inner magnetosphere, where the density is quite low (i.e., the boundary value is 28 AMU/cc) and the magnetic field strength is quite large (i.e., 31 100 nT).The time-step of all hydrodynamic and MHD codes that do not have an implicit time-stepping scheme is limited by the minimum of the size of each cell divided by the maximum wave propagation speed in the cell.For example, with 1/4 R E resolution in the solar wind, and a solar wind speed of 400 km/s, the maximum allowable timestep in that particular cell would be around 3.9 s.In the inner magnetosphere the Alfvén speed is about 30 000 km/s and the cell size may be as small as 1/16 R E ≈ 400 km, this time-step is typically reduced to approximately 0.01 s, meaning that in order to do a one day simulation, 8.6 million iterations have to be taken.There are only a few known ways to get around this limitation.
One method involves using subcycling, in which cells that need to take smaller time-steps may take a few iterations for each of the other cells that have a larger time-step (e.g., Janhunen, 1996).Another significantly more complex method involves using an implicit solver (Tóth et al., 2005b).A final method, which is much easier to implement, is to solve the semi-relativistic MHD equations (Boris, 1970;Gombosi et al., 2002), and artificially lower the speed of light by some fraction.This, in effect, limits the maximum wave speed within the MHD code to be less than the artificial speed of light.When done consistently, it can be shown that the steady state solution remains the same irrespective of the artificial lowering of the light speed (Gombosi et al., 2002).For time-dependent simulations, lowering the speed of light is not physically correct, but it does not appear to affect the overall structure of the magnetosphere.
If the maximum Alfvén speed within the simulated magnetosphere is approximately 0.1c, then setting the "Boris correction" to a number less than 0.1, the time-step within the MHD code can be increased.Using a factor of 0.01 "Boris correction" allows the code runs approximately 10 times faster than without it.This allows simulations to be conducted at a reasonable resolution within a reasonable amount of time, and is therefore one of primary mechanisms that modelers have of making simulation of the magnetosphere in a timely fashion.
Interestingly, there is an additional side benefit of the "Boris correction".All TVD schemes have numerical resistivity that is dependent upon the fastest wave speed within the cell (e.g., Powell et al., 1999a;Lyon et al., 2004) Fig. 8.The top figure shows the ionospheric cross polar cap potential as a function of run-time, with the base-line 2nd order scheme (dashed) and the 2nd order scheme with the semi-relativistic Boris equations, with the speed of light lowered to 1% of the true value (solid).The smallest grid cell size (in R E ) is indicated above the solid line.The bottom figure shows the ratio between the cross polar cap potentials derived from the simulations using Boris and no Boris as a function of run-time.
reducing the Alfvén wave speed within the inner magnetosphere, the amount of numerical diffusion is decreased, reducing the amount of field-aligned current that is diffused away in the inner magnetosphere, and increasing the cross polar cap potential becomes larger.Because the Boris correction only reduces the fastest wave speeds in regions in which the wave speeds approach the (reduced) speed of light, it only reduces the diffusion in the inner magnetosphere, where the magnetic field and the Alfvén speed are quite large.
Figure 8 shows how the CPCP compares when using the "Boris correction" and when not using it.These simulations were run with a factor of 0.01c (starting at 1.0c at 0 iterations, and ending up with 0.01c at 5000 iterations.)The simulation is compared against the baseline case of the 2nd order Rusanov scheme with a minmod limiter.The only difference between the simulations is the Boris correction.It is immediately apparent that artificially reducing the speed of light (and therefore the Alfvén velocity) significantly en-hances the CPCP at moderate resolutions (i.e., 50% increase at 1/4 R E ).At higher resolutions, the simulation is becoming grid-converged, which means that the numerical diffusion has little influence over the simulation, and the Boris correction has little effect on the CPCP.The Boris correction then just increases the time-step of the code.
Typically, though, the code is run with 1/4 R E grid cells, rather than 1/8 or 1/16 R E grid cells due to the dramatic increase in the number of cells when running at higher resolutions.If the maximum resolution that is considered is 1/4 R E , it is clear by comparing all of the previous plots that the Boris correction gives the potential that is closest to the actual value (about 110 kV, compared the to expected 130 kV).By combining the Boris correction with a better limiter, an even higher CPCP can be achieved at comparable resolutions.The Roe scheme can not be used for the semi-relativistic MHD equations, because the eigen vectors of the equations cannot be analytically determined, and these would be needed for the Roe solver.

Inner magnetospheric density
One of the questions that is strongly being debated within the magnetospheric community at this time is the role of ionospheric outflow in the dynamics of the magnetosphere (e.g., Glocer et al., 2009).This seems like it is a perfect place for magnetospheric models to make significant progress, but there are difficulties with this.For example, most MHD models are not multifluid codes, so the ionospheric outflow needs to be specified as a mix of oxygen and hydrogen within the mass density used as a boundary condition within the code.Once the mass leaves the boundary, it is impossible to determine how much of the fluid is oxygen or hydrogen in a single fluid simulation.While this is not a great limitation for basic investigations of the effects of increased mass density on the magnetosphere, many important aspects of the inner magnetosphere cannot be investigated this way.Global magnetospheric codes that couple to inner magnetospheric codes that can differentiate between oxygen and hydrogen (e.g., De Zeeuw et al., 2004;Toffoletto et al., 2004) can not pass the correct information on the individual species densities, so the ring current dynamics can not be modeled accurately.Multispecies and multifluid MHD models are emerging that can actually do this (Winglee, 2000;Glocer et al., 2009;Welling and Ridley, 2010b).
For studies that simply focus on MHD considerations of the outflow, lack of constituent information may not be a hindrance.For example, if one wants to investigate how the change of mass density with ion outflow changes the timing of events within the magnetosphere, or perhaps the effect of ion outflow on the reconnection rate (in resistive MHD), then including the ion outflow within the single fluid codes may help to address these issues.
Within the single-fluid version of the BATSRUS code, ion outflow can be approximated by altering the density on the inner boundary.This density is diffused away from the boundary by a couple of cells, such that it gets into the regions in which pressure gradients can cause it to flow away from the boundary.Within the open field-line region, the lobe pressure is extremely low, so the density near the boundary is advected away by field-aligned flows that develop due to large pressure gradients (Winglee, 2000).In other regions, where the pressure gradients are radially inward, such as on field-lines that map to the tail reconnection site, there is no outflow.
Figure 9 shows the ionospheric cross polar cap potential as a function of iteration number when the BATSRUS inner boundary density is increased by a factor of two.This leads us to believe that, when the (approximate) outflow rate is increased, the ionospheric potential increases, which is contradictory to investigations by Winglee (2000).What appears to be occurring within the BATSRUS code is that the Alfvén speed is decreasing in the inner magnetosphere due to the increased densities.This causes the diffusion to be reduced, as described above, allowing more current to reach the ionosphere, and the potential to be increased.This can be observed in the simulation results by examining the solid line in the bottom plot of Fig. 9.When the cells are large, and the numerical diffusion is large, the ratio between the two simulation results are the largest (approximately 20%).As the resolution increases, and the numerical diffusion is reduced, the ratio between the two simulations reduces to about 5%.So, increasing the ion outflow rate within a non-gridresolved simulation may have a mixture of numerical effects and physical effects.This is quite important to consider when attempting to understand the physical ramification of more mass within global MHD codes.

Resolution on the dayside
Figure 10 shows the MHD simulation results with an improved grid resolution at the dayside magnetopause.As the resolution near the body is improved (as is the case with all of the runs described above), the resolution at the reconnection site is also improved, matching the grid resolution near the body.This means that during the final 12 000 iterations, the resolution near the subsolar point is 1/16 R E .
During the time in which the resolution is 1/16 R E , there are magnetic islands that form near the reconnection site and propagate away from it.These are similar to flux transfer events, and have been observed in other simulations under similar grid resolutions (Dorelli, 2007;Kuznetsova et al., 2009).Additionally, it has been shown that when the grid resolution is fine enough, wave structures appear on the magnetopause (Fairfield et al., 2007).It is clear that in an ideal world, it would be good to resolve the magnetopause with extremely fine resolution all of the time.The problem with this is that it is impractical due to the lack of computational power to handle the extremely large amount of grid cells that it would take to commonly run with 1/16 R E near the magne- topause.If one were to just approximate the dayside magnetopause as a half of a sphere being some where between 8 and 12 R E in radius, it would require 10 million cells to resolve just this region of space at 1/16 R E , not counting the rest of the simulation domain.Additionally, the magnetopause can not really be approximated as a sphere -it is much closer to a paraboloid, which would require even more cells.The adaptive mesh refinement in BATSRUS can match this shape quite well.The simulation by Fairfield et al. (2007) refined the magnetopause with 1/16 R E resolution, which required approximately 20 million total grid cells.
Another issue with refining the grid to extemely small resolutions is that the MHD equations can begin to break down when the spacing becomes close to the ion gyro radius.In regions of strong magnetic field, this is not likely to occur, but in regions of weak fields, like the cusp or the tail, the fine grid resolution could approach the ion gyro radius.
While the resolution assists in creating more wave structures at the magnetopause, it also has an effect on the SWMF to describe the magnetospheric configuration at the present time.
Table 1 describes the simulations conducted here.Each simulation is a single setting different than another simulation, Fig. 10.The current into and out of the plane (i.e., J y in µA/m 2 ) plotted as a color contour plotted under the grid and traces of magnetic field in the Y = 0 plane.The grid is once again shown as the block boundaries, where each square contains 4 × 4 cells.From upper left to bottom right, the highest resolution regions are reconnection rate, and therefore the cross polar cap potential.Figure 11 shows results of the CPCP as a function of iteration number, comparing the base-line simulation with the higher dayside resolution simulation.The cross polar cap potential is reduced in this case for almost all resolutions in which the dayside is resolved better (i.e., 1/4 R E and finer).The reason for this is that ideal MHD is not supposed to be able to model reconnection, since there is no resistivity in ideal MHD.Obviously, when the equations are discretized, and are implemented in various schemes, there is numerical diffusion, which acts as a (nonphysical) resistive term in the MHD equations.Therefore, as the solution moves towards infinite resolution, the model should disallow reconnection, which should show up as a reduction in the reconnection rate.This is what is observed in this simulation, as well as the simulation utilizing the Roe scheme, described above.In the 1/16 R E simulation, the reconnection rate reduces to the point that flux builds up, forcing the rate to increase, which in turn causes periodic flux transfer events to occur.

Time-dependent simulations
In the previous section, steady-state simulations were examined to show how an idealized simulation reacts to different resolutions, solvers, and other numerical factors.Here simulations of a realistic event (4 May 1998) are described.These simulations demonstrate some of the characteristics that were discussed above.Figure 12 shows the most relevant solar wind and IMF inputs into the MHD code that were used as a time-dependent upstream boundary condition on the MHD simulation as it ran.The studies by Yu and Ridley (2008); Wang et al. (2008b) and Welling and Ridley (2010a) show many validation results of this time period that demonstrate how the code compares against other data sets not shown here.
For these simulations, three models were coupled together in the Space Weather Modeling Framework (Tóth et al., 2005a): (1) BATSRUS, which is described above; (2) an ionospheric electrodynamics model, described by Ridley et al. (2004); and (3) the Rice Convection Model, which describes the inner magnetospheric dynamics of low and medium energy particles (Wolf et al., 1982;Sazykin et al., 2002;De Zeeuw et al., 2004).These are the typical models that are used in the SWMF to describe the magnetospheric configuration at the present time.
Table 1 describes the simulations conducted here.Each simulation is a single setting different than another www.ann-geophys.net/28/1589/2010/Ann.Geophys., 28, 1589-1614, 2010 simulation, so most of them can be directly compared.The grid resolution, solver, limiter, speed of light (Boris correction), inner boundary density, and dt in the implicit scheme are varied to determine their effect on the solution, which is quantified through the normalized root mean squared (nRMS) comparison between the output and Dst, ionospheric cross polar cap potential and GOES 08 and 09 satellite magnetic field data.The nRMS error is defined as where M i is the model result at a given time and D i is the data at that particular time.The calculation is done over the first seven hours of the simulation, so that the majority of the nRMS values are computed over the exact same time period.The time column shows how much of the 12 h simulation was completed within the given 8-h queue.Each simulation was conducted on 64 Intel Xeon 2.66 GHz processors.Two of the simulations, shown in bold, crashed while running.The nRMS results are calculated only for the runs that reached seven hours.The first set of simulations (1-6) compare the effects of grid resolution on the solution.Figure 13 shows model results at about the Dst minimum for these six simulations along with the grid resolution.The low resolution grid has 1/2 R E cells at the inner boundary; the medium has 1/4 R E cells; and the high has 1/8 R E cells.The med, med2, med3, and med4 grids have the 1/4 R E grid resolution extending to a spherical radius of 4, 6, 8, and 10 R E , respectively.The high run has 1/4 R E grid resolution to 10 R E with an 1/8 R E shell around the boundary extending to 4 R E .Examining Table 1, the clearest result from this group of simulations is that the amount of time that can be simulated with an 8h queue decreases with increasing resolution.The biggest loss of simulation time comes when the resolution is changed from 1/4 R E cells (run 5, med4) to 1/8 R E cells (run 6, high).This is because the time-step in the code has to be cut in half.While the implicit time-stepping in BATSRUS is being utilized in this set of runs, the number of iterations to converge in the implicit code goes up when the cell size decreases.This is primarily caused by the ratio between the CFL timestep and the requested time-step decreasing by a factor of two.
When the nRMS values are compared, there is not a clear trend.For the ideal simulations, there was a tendency that the higher resolution simulations produced larger cross polar cap potentials, and therefore the simulations were perhaps "better".This was not actually the case, though, since the real cross polar cap potential was supposed to be around Ann. Geophys., 28,2010 www.ann-geophys.net/28/1589/2010/130 kV with B z = −10 nT.So, some of the simulations were clearly overestimating the CPCP, and there was an ideal resolution for getting the CPCP "correct", even though it was most likely for the wrong reasons.The same trend is observed here -with higher resolution, the results are not necessarily "better", as defined by the nRMS.
In Fig. 13, there is a huge amount of variation on the amount of pressure in the inner magnetosphere and the amount of stretching on both the dayside and nightside.These images are at one instant in time, and so no assumptions about equilibrium can be made; indeed, this time is at the peak of a storm, so the solutions are most likely extremely www.ann-geophys.net/28/1589/2010/Ann.Geophys., 28, 1589-1614, 2010  1.
dynamic.The lowest resolution simulation has the highest degree of stretching on the nightside, while simultaneously having the lowest inner magnetospheric pressure.By adding a very small layer of 1/4 R E cells near the inner boundary, the pressure in the inner magnetosphere is increased by a factor of two, but the stretching is reduced significantly.The stretching is increased by expanding the resolution into the outer magnetosphere, but it is decreased again once the high resolution region reaches well beyond the reconnection site.When a layer of 1/8 R E cells is added at the inner boundary, the stretching increases again.There is a large amount of stretching on the dayside as well as the nightside in the 1/8 R E simulation, such that the cusps are at very low latitudes.The tendencies of stretching match the results in Table 1, when comparing to GOES data.By examining these plots, it is clear that the magnetosphere is highly nonlinear and that the nonlinearity is driven by both the solar wind and IMF as well as the grid resolution and the different numerics in the code.
To get a more global view of the simulations, the Dst index gives a large-scale average of how much current is running through the inner magnetosphere.This is a simplification, though, since a smaller ring current closer to the Earth can give the same Dst as a larger ring current farther from the Earth.For this type of detail, the GOES satellites can be examined.During these periods, there is a large amount of stretching in the field, which indicates the presence of a large ring current at or outside of geosynchronous orbit.If the simulation were to model both Dst and the GOES satellite magnetic field measurements well, it would show that both the size and the shape and location of the ring current are correct.For example, as both Fig. 13 and Table 1 indicate, Run 1 has more stretched field-lines in the inner magnetosphere than Run 2, indicating that it most likely has a better field-line shape in the midnight region, but has a significantly worse Dst comparison and a smaller pressure peak on the nightside, indicating that the total amount of particle energy in the inner magnetosphere is most-likely incorrect.It is most likely that an enhanced tail current is causing the stretching of the field, but not enough of a ring current to cause the correct Dst enhancement.
Figure 14 shows measured and simulated Dst indices during this storm period for all of the runs described in Table 1.None of the lines are labeled, so it is impossible to tell which run is which.The point of this figure is to show that there is significant variation between the model runs.The dashed line represents the measured Dst index.Some of the simulations are shown to reproduce the depth of the Dst well during the storm, while others more accurately reproduce the rate of Dst decrease during the injection phase.Still others do neither and completely miss the majority of the storm.
Figures 15 and 16 show comparisons between the GOES magnetic field measurements (dashed lines) and the different model runs (colored lines).These are once again illustrative figures that show the diversity in the model results, some of which do much better than others.Almost all of the model results share one common feature, though -there is too little stretching of the field.This can be observed by noticing that B x is almost always underestimated, while B z is overestimated.The comparison with the GOES-08 satellite (Fig. 15) shows that the trends in B y are extremely well modeled by all of the runs.This is true of the GOES-09 satellite also, but the model results are all too low in magnitude.
Runs 7a and 7b switch to more accurate solvers.These runs are directly comparable to run 4.Both of these simulations are slower, and they do not appear to perform as well in any of the comparisons, which is unexpected, because both of the utilized solvers are superior to the Rusanov solver.This illustrates the point that, while the MHD equations may be solved for in a more-accurate fashion, the comparisons with the data can become worse in some cases.
Runs 8 and 9 change the limiter to a lower and higher value, respectively.These simulations are directly comparable to run 4, again.Interestingly, having a more diffusive limiter allows the simulation to run faster, most likely having to do with the rate of convergence in the implicit scheme.The more diffusive limiter has a significantly worse Dst and a slightly worse comparison with GOES and CPCP.The worse inner magnetospheric results are expected, since the weak limiter doesn't allow strong gradients to exist, which are needed to drive a realistic ring current.The stronger limiter simulation (run 9) crashed.Typically crashes occur because the pressure in the lobes becomes negative.Once again, BAT-SRUS does not put a floor on either the pressure or the density in any region of space, so it can crash due to negative values in primary variables that should be positive.The strong limiter demonstrated in run 9 is never actually used in modeling the magnetosphere (in BATSRUS), since it does allow very strong gradients to form, which can result in negative thermal pressures in the lobes, as it did in this case.If a floor were put on the pressure values, the simulation would not show the location of the satellite in GSM coordinates, with the diamond, star and triangle indicated in all plots to show how the orbit corresponds to the data below.The bottom three plots show, from top to bottom, the Bx, By and Bz components of the magnetic field measured by the satellite (dashed black lines) and values from the simulations (colored lines).Each color represents a different simulation described in Table 1.
crash, but the results may be questionable.In this way, BAT-SRUS may be a less robust code, but it is more self-consistent and does not need to artificially limit values for storm-time simulations.
The next two simulations (runs 10 and 11) modify the speed of light to lower values (0.01c and 0.005c, respectively), to show that this can reduce the diffusion in the code significantly.Run 10 allows much more of a ring current to form (compared to run 4), although the location and size may not be as accurate.The cross polar cap potential is also not quite as accurately simulated.The amount of time that it simulated is reduced, due to the convergence in the implicit scheme.Run 11 crashed very soon after the storm started.Running with an extremely low speed of light (< 3000 km/s) makes the code less stable and it typically will not work during a storm simulation.
Runs 12 and 13 simply change the time-step in the implicit scheme to 2.5 and 10 s, respectively.A shorter timestep gives a better temporal resolution, and typically requires fewer iterations to converge; while a longer time-step gives worse temporal resolution and requires more iterations for convergence.The resolutions must be balanced with the number of iterations for convergence, and there is typically an optimal time-step that results in the fastest run-time.The 5 s time-step allows slightly more time to be simulated (10:06) than the 10 s time-step (09:54) and significantly more that the 2.5 s time-step simulation (08:39).Interestingly, the Dst and CPCP have better comparisons for the 10 s time-step simulation (compared to run 4).
Runs 14 and 15 increase the mass density at the inner boundary.The simulation with 56 AMU/cm 3 shows a significant improvement over the simulation with 28 AMU/cm 3 as an inner boundary density.There are two possible explanations for this: (1) the increased density allows more energy to flow into the ring current, thereby causing more distortion of the field and a better Dst; and (2) the increased density reduces the Alfvén wave speed, which reduces the diffusion in the code, thereby allowing a more accurate simulation to occur.While the second reason has been tested by increasing the resolution and reducing the speed of light, the results have been mixed and have shown that simply reducing the diffusion doesn't necessarily decrease the nRMS values.By increasing the density, the amount of energy going into the ring current significantly increases.Dst is much better simulated and the stretching of the field-lines is better modeled (as shown by the comparisons to GOES).When the density is further increased to 112 AMU/cm 3 , the www.ann-geophys.net/28/1589/2010/Ann.Geophys., 28, 1589-1614, 2010 amount of energy in the ring current becomes too large, and the Dst is over-predicted.This can be observed in Fig. 14, but looking at the lowest two lines.The bottom line is the 112 AMU/cm 3 simulation, while the one directly above this line is the 56 AMU/cm 3 simulation.The 112 AMU/cm 3 simulation obviously over-predicts the Dst, while 56 AMU/cm 3 simulation matches the Dst better than any other simulation.Interestingly, the stretching of the magnetic fields in the 112 AMU/cm 3 simulation (in the comparison with the GOES satellites) is still not enough, but is better than all of the other simulations.This can also be seen in Figs. 15 and 16 -this simulation is the closest to the data.The mixed result once again points to the fact that the code can get roughly the right amount of energy into the ring current, but it can be at the wrong place.These results clearly show that the inner boundary density is extremely important in the simulation of the magnetosphere, and that a more accurate description of the ionospheric outflow is warranted, as is further described by Glocer et al. (2009).Welling and Ridley (2010b) found that, in order to more precisely match dayside density measurements during a relatively quiet time, the inner boundary density needed to be lowered significantly.This points to the fact that the mass may be quite important in the simulationsduring storms the mass density needs to be higher because of the oxygen flowing out of the ionosphere, while during quiet times, the mass density needs to be lower because of the lack of oxygen outflow.Runs 16-18 don't utilize the implicit feature of BATSRUS, and are meant to show how much of an acceleration the implicit code gives.Run 16 and run 4 are directly comparable.Run 16 only simulates 04:17 of the storm, while using the implicit code, 10:06 is simulated, for a speed up of 236%.Because so little of the storm is simulated, the nRMS values can not really be compared.Run 17 shows that the implicit code speeds up the 1/8 R E simulations also (Run 6), but this time by a smaller amount 199%.By reducing the speed of light a little bit more (run 18), the explicit and implicit times are comparable with 1/8 R E resolution.This is because the number of iterations for convergence in the implicit solver becomes larger as the cell size becomes smaller.This is a nonlinear effect, while the reduction of the timestep in the explicit code is linear as the grid resolution is increased.Therefore, a balance is reached between the implicit and explicit Boris correction methods around 1/8 R E resolution, where the implicit code is faster for lower resolutions and the explicit Boris correction scheme is faster for higher resolutions.
Run 19 utilizes the best solver that is available within BAT-SRUS, the Roe solver.While the Roe solver is the least diffusive of all of the solvers within BATSRUS, the semirelativistic MHD equations can not be solved for using the Roe solver, which means that the speed of light can not be artificially reduced.The run with the Roe solver can be directly compared to run 8.It is clear that the Roe solver does not do as well as the Rusanov solver with the speed of light reduced.Indeed, some of the worst comparisons in Figs.14-16 are from the Roe solver (this run is the bright yellow line furtherest away from the data in almost all cases).While this solver is extremely powerful, it shows that the Boris correction is crucial when simulating the inner magnetosphere.Further, it may indicate that, if the numerical dissipation is reduced in the reconnection site, simulations of realistic time-periods may be degraded, due to the significant reduction of the reconnection rate.More tests need to be conducted to determine the true cause of the disappointing results with the Roe solver.
Run 20 is directly comparable to Run 3, but moves the zero potential ionospheric boundary condition from 50 • latitude to 10 • latitude.This condition allows the potential to penetrate to low latitudes, allowing more plasma to flow closer to the Earth in the magnetosphere.While this shouldn't affect the plasma flow within the MHD domain very much at all, it does affect the plasma flow in the RCM domain, which extends much closer to the Earth than 2.5 R E .The net result of moving the boundary is that more plasma is allowed to be energized in the inner magnetosphere, strengthening the ring current, as evidenced by the decrease in the nRMS error for the Dst comparison (compare Run 3 and Run 20).The magnetic field-lines end up being more stretched also.

Discussion
Within the magnetospheric community, there have been very few direct model-to-model comparisons.These types of comparisons are typically conducted at workshops and almost always result in disappointment by the audience, since it is almost impossible to determine why the models give different results.While the comparisons range from complicated storm time periods to relatively simple input conditions, the results sometimes disagree and sometimes agree, depending on the models.Because BATSRUS and GUMICS-4 are similar in many ways in terms of numerical scheme and adaptive mesh, they tend to agree more than the other MHD codes, although this is not quantitative.
All of the results that are presented above are produced by a single MHD code, namely BATSRUS.Given that the code is solving the same set of equations (basically, the semi-relativistic MHD equations), the results can be different given the choice of different solvers, resolutions, limiters, Boris factor, implicit time-step and inner magnetospheric density.If the choice is expanded further and different codes solving differently discretized equations (for example, using a conservative versus non-conservative form of the MHD equations), the solutions can diverge even more.
The MHD equations are nonlinear.The magnetosphere is also nonlinear.This means that by examining even one simple quantity, such as the ionospheric cross polar cap potential, a wide variety of numerical and physical effects can influence it.For example, as shown above, the grid resolution in the inner magnetosphere, the Boris correction and the resolution at the magnetopause all affect the CPCP in ways that may constructively or destructively interact with each other.The resulting CPCP may be the value desired (when compared to data), but it may be for the wrong reasons: if the dayside magnetopause is not well resolved, leading to a high reconnection rate and relatively large region-1 currents, and the inner magnetosphere is not well resolved, allowing diffusion of those currents, the net current into the ionosphere can be about what is expected, resulting in a nominal CPCP.Further, it is shown that by adding more grid resolution and reducing the diffusion within the MHD code, the solution may not actually be more accurate (when compared to data).This is most likely due to missing physics within the code and fact that as numerical diffusion is decreased, this missing physics becomes the dominant source of the errors (as it should).
Researchers who are attempting to conduct scientific investigations with a global MHD code, may consider running the simulation with the same drivers over a few times varying different parameters within the code, such as grid resolution, solvers and boundary conditions to determine the robustness of the physical result.If the result of interest only shows up for a single run with a given grid and solver (for example), this should be made public, so the scientific community can judge for themselves the validity of the result.
Given that researchers would like to compare codes and understand the differences between them, it may be best to simplify the problems as much as possible (e.g.constant, nominal inputs, no anomalous resistive terms, and constant ionospheric conductances), and make sure that the models are run with as close to the same grid resolution and same Boris factors as possible.It may then be best to compare topological features, such as magnetopause location and shape.A variety of conditions could be simulated (always using constant drivers or simple step-changes), such that researchers can examine the model's dependence of the topological feature on the input conditions.Once the simple features of the code are well understood, then researchers can investigate features interior to the magnetosphere.Comparing features such as the field-aligned current strength and topology (for example) during varying input conditions may be too challenging of a problem, given the dependence on the complexity of the numerics involved.
Finally, when comparisons are made to data during different events, modelers should make it clear how they ran their model, specifying quantities such as the resolution in different regions, the Boris correction, solver, inner boundary density, etc.These should remain the same for each comparison, so the model can be evaluated fairly.

Conclusions
The BATSRUS code is versatile, having the ability to simulate many different planetary bodies and general plasma conditions.This means that BATSRUS is a useful tool when studying different regimes within the solar system.As this study shows, and other validation studies have shown, BAT-SRUS is capable of accurately modeling the magnetospheric dynamics (e.g., Wang et al., 2008b;Yu and Ridley, 2008;Taktakishvili et al., 2007;Wang et al., 2008a;Fairfield et al., 2007;Welling and Ridley, 2010a;Yu et al., 2010).Further, BATSRUS is versatile enough that it can be used to simulate the magnetosphere under a wide variety of user needs (i.e., resolution, run-time, etc.), as is shown here.Further, this flexibility of the code means that BATSRUS can be utilized to study the effects of numerics on magnetospheric simulations.Very few studies showing how realistic magnetospheric simulations are affected by the numerics of the code have been published.These types of studies are extremely important, since they show how the different numerics can interplay with the physics, sometimes confusing the interpretation.For example, this study shows that when the grid is resolved more in the inner magnetosphere, the amount of field-aligned current into the ionosphere increases, while, when the resolution at the dayside magnetopause is increased, the field-aligned current decreases.These are both numerical effects that have to be understood in order to explore magnetosphere-ionosphere coupling due to reconnection or flux transfer events.
This study shows that the inner magnetosphere is not grid converged for nominal, and even high, resolutions.This is because of the strength of the dipole and the fact that the dipole field varies as r −3 , implying that the grid resolution should increase at a similar rate, if field-aligned currents are to be fully resolved from their generation region to the inner boundary.This is not realistic in global MHD codes because of time-step limitations (which can be partially overcome with subcycling, implicit solvers and the Boris correction) and just the shear number of grid cells that would be needed to satisfy the resolution needed.This lack of grid convergence is only truly important if the science that is being investigated is missed without the convergence.For a space weather type of application, where the prediction of something happening is more important than the exact details, the lack of grid convergence is most likely perfectly fine.For studying localized field-line resonances, even the highest resolutions are not good enough since the diffusion will damp the wave structures.
It is shown that by artificially reducing the speed of light within the inner magnetosphere (also known as the Boris correction), much of the diffusion can be eliminated, and reasonable grids can be utilized.In idealized simulations, using a Boris correction of 0.01c is roughly equivalent of having double the resolution in the inner magnetosphere.The Boris correction has very little influence when the grid is www.ann-geophys.net/28/1589/2010/Ann.Geophys., 28, 1589-1614, 2010 high enough resolution that numerical diffusion is already significantly reduced, and therefore it is only really needed in models that have limited resolution (i.e., larger grid-cells than 1/8 R E ) in the inner magnetosphere.It is a common method utilized by almost all global MHD modelers, including BATSRUS, and is most likely the best way of reducing the numerical diffusion in the inner magnetosphere.Further, it greatly increases the time-step size in the simulation, thereby reducing the run-time by as much as an order of magnitude (ignoring the use of the implicit code).The density at the inner boundary is critical in determining how much energy density will enter the ring current.It is found that utilizing a constant value of 28 AMU/cm 3 gives reasonable results, while doubling the density to 56 AMU/cm 3 improves the comparisons with data dramatically.Increasing the value further tends to put too much plasma into the ring current and over-predicts the value of Dst, although the stretching of the field is more realistic.Glocer et al. (2009) found that utilizing a more self-consistent outflow significantly improved the stretching of the magnetic field-lines in the inner magnetosphere.
When the ionospheric zero potential boundary condition is moved from 50 • latitude to 10 • latitude, the ring current is shown to be more accurately modeled, having both a stronger Dst and more stretching observed at geosynchronous orbit.This is because the low latitude boundary condition allows more particles to flow deeper into the Rice Convection Model domain, which energizes them more and allows more energy density into the inner magnetosphere.
Finally, this study shows that the numerics of global MHD codes have a large, non-linear, influence on the solution.Researchers must be aware of these effects when attempting to ascribe a physical interpretation to the modeling results.

Fig. 1 .
Fig.1.Magnetospheric current density into and out of the plane (i.e.J y in µA/m 2 ) plotted as a contour with the grid structure overlayed for different iterations numbers through the simulation.The grid boxes indicate the block boundaries, which means that within each of the plotted squares, there are 4 × 4 cells.These plots show the different resolutions that are investigated in this study.

Fig. 2 .
Fig. 2. The ionospheric field-aligned currents (left, in µA/m 2 ) and potentials (right, in kV) in an MLT-magnetic latitude coordinate system, with noon at the top and dawn on the right.The outer ring is 50 • latitude.The maximum and minimum values for each plot are indicated.The rows correspond to the various grid resolutions from 1 R E (top) to 1/16 R E (bottom).

Fig. 4 .
Fig. 4. The top figure shows the ionospheric cross polar cap potential as a function of run-time, with Rusanov (dashed) and Linde (solid) second-order solvers.The smallest grid cell size (in R E ) is indicated above the solid line.The bottom figure shows the ratio between the Linde and Rusanov cross polar cap potentials as a function of run-time, and the percentage of the final CPCP for the Linde scheme as a dashed line.

Fig. 5 .
Fig. 5.The top figure shows the ionospheric cross polar cap potential as a function of run-time, with Rusanov (dashed) and Roe (solid) second-order solvers.The smallest grid cell size (in R E ) is indicated above the solid line.The bottom figure shows the ratio between the Linde and Roe cross polar cap potentials as a function of run-time, and the percentage of the final CPCP for the Roe scheme as a dashed line.

Fig. 6 .
Fig. 6.Top: plots of the pressure color contours in the Y = 0 plane, with traces of the magnetic field in this plane overplotted.The simulation results using the Rusanov (Roe) solver is on the left (right).Bottom: plots of the field-aligned current (left) and ionospheric potential (right) for the Rusanov (top) and Roe (bottom) simulation.The color contours are the same for each pair of Rusanov and Roe plots.

Fig. 7 .
Fig. 7.The top figure shows the ionospheric cross polar cap potential as a function of run-time, with a minmod limiter (dashed) and MC limiter, with β = 1.2 (solid).The smallest grid cell size (in R E ) is indicated above the solid line.The bottom figure shows the ratio between the cross polar cap potentials derived from the simulations using the MC and minmod limiters as a function of run-time.

Fig. 9 .
Fig.9.The top figure shows the ionospheric cross polar cap potential as a function of run-time, with an inner boundary density of 28 cm −3 (dashed) and 56 cm −3 (solid).The smallest grid cell size (in R E ) is indicated above the solid line.The bottom figure shows the ratio between the cross polar cap potentials derived from the simulations using 56 cm −3 and 28 cm −3 as inner boundary density conditions as a function of run-time.

Fig. 10 .
Fig.10.The current into and out of the plane (i.e., Jy in µA/m 2 ) plotted as a color contour plotted under the grid and traces of magnetic field in the Y = 0 plane.The grid is once again shown as the block boundaries, where each square contains 4x4 cells.From upper left to bottom right, the highest resolution regions are 1 Re, 1/2 Re, 1/4 Re, 1/8 Re and 1/16 Re.

Fig. 11 .
Fig.11.The top figure shows the ionospheric cross polar cap potential as a function of run-time, with the base-line 2nd order scheme (dashed) and a simulation with similarly high resolution on the dayside (solid).The smallest grid cell size (in R E ) is indicated above the solid line.The bottom figure shows the ratio between the cross polar cap potentials derived from the simulations using a refined dayside grid and a fixed dayside grid.

Fig. 12 .
Fig. 12. From top to bottom: The solar wind density and temperature, the IMF B y and B z and the ionospheric cross polar cap potential from AMIE and the different simulations described in Fig. 1 on 4 May 1998.

Fig. 13 .
Fig. 13.Plots from Runs 1-6 showing the pressure and magnetic field topology at 04:40 UT on 4 May 1998.The grid resolution is indicated by the plots of the block boundaries, as in previous figures.

Fig. 14 .
Fig. 14.The measured (dashed black line) and simulated (colored lines) Dst index for the 4 May 1998 storm.The simulated results are described in Table1.

Fig. 15 .
Fig. 15.Comparison with the GOES-08 satellite.The top two plots show the location of the satellite in GSM coordinates, with the diamond, star and triangle indicated in all plots to show how the orbit corresponds to the data below.The bottom three plots show, from top to bottom, the Bx, By and Bz components of the magnetic field measured by the satellite (dashed black lines) and values from the simulations (colored lines).Each color represents a different simulation described in Table1.

Fig. 16 .
Fig. 16.Comparison with the GOES-09 satellite in the same format as Fig. 15.

Table 1 .
A description of the 21 different simulations that were conducted for this section of the study.The column "time" shows the final simulation time corresponding to 8 h of execution.The nRMS errors only include the first seven hours of the simulation (indicated by the vertical lines in Figs.14-16, to make sure that most of the simulations can be compared equally.The bold-faced simulation numbers are runs that crashed because of negative pressures (typically in the lobes).Further details are described within the text.
a Moved the ionospheric boundary condition to 10 • latitude instead of ∼ 50 • latitude.