Comparison of the observed dependence of large-scale Birkeland currents on solar wind parameters with that obtained from global simulations

Spatial distributions of the large-scale Birkeland currents derived from magnetic field data acquired by the constellation of Iridium Communications satellites have been compared with global-magnetosphere magnetohydrodynamic (MHD) simulations. The Iridium data, spanning the interval from February 1999 to December 2007, were first sorted into 45 -wide bins of the interplanetary magnetic field (IMF) clock angle, and the dependencies of the Birkeland currents on solar wind electric field magnitude,Eyz, ram pressure, psw, and Alfvén Mach number, MA , were then examined within each bin. The simulations have been conducted at the publicly-accessible Community Coordinated Modeling Center using the University of Michigan Space Weather modeling Framework, which features a global magnetosphere model coupled to the Rice Convection Model. In excess of 120 simulations with steady-state conditions were executed to yield the dependencies of the Birkeland currents on the solar wind and IMF parameters of the coupled model. Averaged over all IMF orientations, the simulation reproduces the Iridium statistical Birkeland current distributions with a two-dimensional correlation coefficient of about 0.8, and the total current agrees with the climatology averages to within 10 %. The total current for individual events regularly exceeds those computed from statistical distributions by factors of ≥2, resulting in larger disparities between observations and simulations. The simulation results also qualitatively reflect the observed increases in total current with increasingEyz andpsw, but the model underestimates the rate of increase by up to 50 %. The equatorward expansion and shift of the large-scale currents toward noon observed for increasing Eyz are also evident in the simulation Correspondence to: H. Korth (haje.korth@jhuapl.edu) current patterns. Consistent with the observations, the simulation does not show a significant dependence of the total current onMA .


Introduction
Global simulations are becoming the research tool of choice for exploring and understanding state and dynamics of planetary magnetospheres.The physical processes in the coupled solar wind-magnetosphere-ionosphere system can be simulated using a variety of modeling techniques, which are coarsely divided into three categories: full particle simulations (e.g., Pritchett, 2003), which treat both electrons and ions kinetically; hybrid simulations (e.g., Winske et al., 2003), in which ion motion is kinetic while approximating electrons as a conducting fluid; and magneto-hydrodynamic (MHD) simulation (e.g., Raeder, 2003), which assume fluid motion of both ions and electrons.The practicability of a simulation tool for a particular environment is predominantly determined by the ability to meet the computational requirements, which scale, among other parameters, with the size of the simulation volume.At Earth, the dimensions of the magnetotail are so large that fully kinetic and hybrid simulations of the entire magnetosphere are computationally unfeasible at this time.However, as a result of recent advances in computing technology and parallelization of the simulation codes, MHD global simulations of the terrestrial magnetosphere are now well within reach, and community access to MHD simulations, such as through the Community Coordinate Modeling Center (CCMC), has become commonplace.
H. Korth et al.: Global Birkeland currents comparison Because of the increased popularity these simulations enjoy in attempting to understand the natural system and because the simulations can only approximate the physics of the coupling between solar wind, magnetosphere, and ionosphere, it is important that the model output is validated against observations on a global scale and for a wide range of boundary conditions.
A number of studies comparing products of MHD models with observations have been carried out in the past.Early studies compared magnetic field and plasma parameters along single spacecraft trajectories (e.g., Frank et al., 1995;Raeder et al., 1997;Palmroth et al., 2001).However, interpreting the results of these comparisons to identify specific areas where the models should be improved proved to be difficult because the observations sampled only a small fraction of the simulation volume.As a result, there remained an ambiguity in distinguishing between missing or incompletely implemented dynamics in the simulation and displacement in time or location of the relevant process in the simulation relative to the point of observation.
Comparisons of simulation results with two-dimensional distributions of auroral emissions offered the advantage of providing coverage that maps to a large portion of simulation space (Goodrich et al., 1998;Rae et al., 2010).However, significant uncertainties in the analysis may result from the indirect relationship between the observed and simulated quantities, namely auroral emissions and particle precipitation, because the physical processes in the auroral acceleration region linking these quantities are not well enough understood to be captured in a predictive model (Carlson et al., 1998;Ergun et al., 2002).Hence it can be difficult to determine whether discrepancies between observed and predicted emissions imply errors in the technique for estimating aurora, missing plasma populations in the simulation, or important differences in magnetospheric electrodynamics.
Two-dimensional distributions of simulated and observed quantities have also been compared by use of integrated quantities.Palmroth et al. (2005) compared the total Joule heating in the Northern Hemisphere from the Grand Unified Magnetosphere Ionosphere Coupling Simulation (GUMICS-4) with estimates computed from electric fields measured by the Super Dual Auroral Radar Network (SuperDARN) and Pedersen conductances inferred from POLAR ultra-violet and X-ray images.These authors found that the temporal variability in that simulation resembles the observations quite well, while the Joule heating magnitude in the simulation is about a factor of 10 lower than observed.Similar underestimates were obtained when the simulation results were compared to Joule heating magnitudes computed from an AEindex-based proxy (Ahn et al., 1983) and using the Assimilative Mapping of Ionospheric Electrodynamics (AMIE) technique (Richmond and Kamide, 1988).
Ground-based observations of magnetic field perturbations induced by toroidal currents flowing in the ionosphere were used to assess model performance (Raeder et al., 2001b;Ridley et al., 2001;Shao et al., 2002;Yu and Ridley, 2008;Yu et al., 2010).From the simulations, the ground magnetic perturbations are computed from two-dimensional Biot-Savart integration over the ionospheric toroidal currents (Fukushima, 1976;Raeder et al., 2001a).The technique depends critically on estimates of ionospheric conductance and is therefore most suitable for use on the dayside, where the conductance is produced predominantly by solar extreme-ultraviolet radiation.Nevertheless, the nearcontinuous availability of magnetic field observations from over 200 ground stations (Gjerloev, 2009), probing the simulation space on a scale similar to that of auroral observations, makes them attractive for evaluating the global models.
More recently, climatology averages of in-situ plasma and magnetic field observations in the magnetosphere from multi-year satellite missions have been used for comparative studies with MHD simulations (e.g., Guild et al., 2008;Daum et al., 2008).These data sets extend the spatial coverage of single satellite trajectories, and the averaged plasma moments and vector magnetic field can be compared directly to the simulation variables.Validation studies based on these data sets have shown that state-of-the-art MHD models are capable of reproducing the large-scale features of the magnetosphere.However, the above climatological data sets can only provide insight into the average conditions within the magnetosphere and do not allow individual states of the solar wind-magnetosphere-ionosphere interaction to be distinguished.
Direct comparisons of the instantaneous magnetospheric state that are (1) global in nature, (2) do not rely on supplemental assumptions to relate simulation results and observations, and (3) allow different solar wind driving conditions to be distinguished, can be made with estimates of field-aligned Birkeland currents.Birkeland currents play a fundamental role in conveying stress between magnetosphere and ionosphere and their distribution reflects the dynamic state of the coupled system.They are therefore an ideal diagnostic quantity to benchmark model performance.The magnetic perturbations associated with the Birkeland currents are sampled globally by the constellation of Iridium Communications satellites (Anderson et al., 2000), which consists of 66 satellites in six 780-km altitude, circular polar orbit planes.Current density distributions are routinely derived from the Iridium engineering magnetometer data using the approach detailed by Waters et al. (2001).The uncertainty of the derived Birkeland current densities was established by Korth et al. (2004).Statistical distributions of the Birkeland currents sorted by IMF orientation were compiled from "stable" current patterns identified between February 1999 and December 2005 with no additional assumptions (Anderson et al., 2008).The results not only confirmed but also improved upon previous statistical models derived from Dynamics Explorer 2 (Weimer, 2001) and Ørsted data (Papitashvili et al., 2002), both of which demonstrated that the IMF orientation fundamentally Ann.Geophys., 29,[1809][1810][1811][1812][1813][1814][1815][1816][1817][1818][1819][1820][1821][1822][1823][1824][1825][1826]2011 www.ann-geophys.net/29/1809/2011/determines the magnetosphere-ionosphere (M-I) coupling geometry.The stable current database was extended through December 2007 by Korth et al. (2010), hereinafter P1, who examined the control of solar wind electric field, dynamic pressure, and Alfvén Mach number on the distribution and intensity of the currents.It was found that the solar wind electric field is the dominant factor determining the Birkeland current intensity and that the large-scale Birkeland currents shift toward the dayside and expand equatorward with increasing solar wind electric field.The current intensity shows only a modest dependence on solar wind dynamic pressure after correcting for the solar dynamo dependencies.Finally, normalizing the Birkeland current densities to both the median solar wind electric field and dynamic pressure effects, there was no significant dependence of the Birkeland currents on solar wind Alfvén Mach number.
Global distributions of the Birkeland currents observed by the constellation of Iridium satellites were previously compared to simulation results from the Lyon-Fedder-Mobarry (LFM) MHD model (Fedder and Lyon, 1995;Fedder et al., 1995;Mobarry et al., 1996) both for steady-state intervals with opposite sign of the IMF B y (Korth et al., 2004) and for magnetic cloud events with slowly rotating IMF (Korth et al., 2008a).The simulations were driven by solar wind (McComas et al., 1998) and IMF (Smith et al., 1998) data from the Advanced Composition Explorer (ACE), located at the first Lagrangian point, L1, starting at least twelve hours prior to the beginning of the respective event period.Comparison of the observed and simulated Birkeland current distributions, which are intimately related to the plasma drifts at the ionosphere, showed that the ionospheric convection pattern in the MHD model and its dependence on the IMF orientation is essentially correct.For southward IMF, the Birkeland total currents in the simulation were about a factor of 2 larger than observed, which was attributed to a combined effect of the simulation overestimating the ionospheric electric field and the Iridium fits underestimating the magnetic perturbations.Finally, in all events examined, the Region-2 current system was found to be largely under-represented because the magnetospheric drift physics was not included in the model.This fundamental difference in the current distribution implies that the path of current closure in the simulation deviates from that of the natural system.Therefore, faithful simulation of the magnetosphere requires the drift physics to be included in magnetospheric global models.
The above comparisons of large-scale Birkeland current observations with LFM simulation results demonstrated the qualification of the Iridium data set to benchmark the capabilities of state-of-the-art magnetospheric models and pinpoint their strengths and shortcomings.However, the studies carried out thus far are not all inclusive and considered only a limited range of the solar wind and IMF conditions.The purpose of this work is to close this gap and assess comprehensively the current state of our modeling capabilities using an MHD model with wide community access.This is accomplished by comparing the observed dependence of the large-scale Birkeland currents on IMF and solar wind parameters to global simulations.The basis of the observations for the analysis below are the statistical results of Anderson et al. (2008) and of P1, which define the M-I coupling geometry and intensity by means of a small set of variables representing the solar wind conditions, namely the IMF clock angle and the solar wind electric field and ram pressure.The simulations were performed using the University of Michigan Space Weather Modeling Framework implemented at the Community Coordinated Modeling Center.The Iridium constellation and database are described in Sect.2, and the global magnetosphere model and simulation runs are presented in Sect.3. The comparison of the observed solar wind control of the Birkeland currents with simulation results are presented in Sect. 4. The results are discussed in Sect. 5 and summarized in Sect.6.

Iridium database
This study uses a database of "stable" current distributions identified in the Iridium data from February 1999 to December 2007.The selection of these intervals from a total of ∼85 000 two-dimensional spatial distributions of the Birkeland currents is described in detail by Anderson et al. (2008) and in P1.Briefly, for each one-hour interval the Iridium magnetometer samples were reduced to magnetic perturbation data as described by Anderson et al. (2000), and Birkeland current patterns with 4 • latitude resolution were derived from these perturbation data using the approach of Waters et al. (2001).Since the magnetosphere can exhibit largescale dynamics on time scales that are shorter than the onehour data accumulation time, intervals suitable for the analysis were carefully selected from the overlap of the dominant large-scale current regions between consecutive onehour distributions (Anderson et al., 2008).A fractional overlap of at least 45 % was required to ensure that the automatic procedure selected only those pairs of distributions that one would visually identify as consistent.In addition, plasma (McComas et al., 1998) and magnetic field (Smith et al., 1998) observations by the Advanced Composition Explorer (ACE) satellite specifying the solar wind and IMF conditions at the first Lagrangian point, L1, were required to be available.The selection criteria were met by a total of 1536 events, corresponding to 4 % of the data.This set of events is the basis for the comparisons described below.
The Birkeland currents are modulated by systematic variations in the ionospheric conductivity, which must be normalized to distinguish the solar wind control of the currents from conductance effects.The two main sources of ionization are solar extreme ultraviolet (EUV) radiation and energetic particle precipitation.Because of the uncertainties in accounting for the contribution from particle precipitation, we account only for conductance variations due to solar illumination (see H. Korth et al.: Global Birkeland currents comparison P1).Solar illumination controls the Birkeland currents primarily by altering the Pedersen conductance, P , which can be expressed in terms of the solar zenith angle as (Rasmussen et al., 1988): where ν = χ /90, u = F 10.7 /90, B is the magnetic field strength, F 10.7 is the 10.7-cm solar radio flux, and χ is the solar zenith angle.We normalized the dayside current densities to equinox conditions by multiplying each Birkeland current distribution, j i (θ,φ), with the ratio of the conductance distributions at equinox, P,eq (θ,φ), and at the time of the observations, P,i (θ,φ).The normalized current density, j * i (θ,φ), is given by: where θ and φ are the magnetic colatitude and local time, respectively.Equation ( 1) is applicable within the solar zenith angle range 0 • ≤ χ ≤ 85 • .For χ > 85 • , the conductance variations due to solar EUV are not expected.The normalization factor is therefore set to smoothly transition to unity as described in P1 so that the current densities on the nightside remain unaffected by the normalization.
The identified Iridium Birkeland current distributions span a wide range of solar wind conditions and can be used to test the steady-state configuration in global magnetosphere simulations.As stated above, the Birkeland currents were previously found to be well correlated with the solar wind conditions as described by the IMF clock angle, α, the magnitude of the solar wind electric field in the plane normal to the Earth-Sun line, E yz , and the solar wind dynamic pressure, p sw .From the Iridium observations, the dependencies of Birkeland currents on solar wind observations could be determined for low to moderate solar wind electric fields, E yz ≤ 4.5 mV m −1 , for which the number of Birkeland current distributions in the database was statistically significant (P1).The corresponding range for the solar wind dynamic pressure is 0 < p sw < 4 nPa.The Birkeland currents did not show a significant dependence on solar wind Alfvén Mach number, M A , when normalized for both electric field and pressure effects.Nevertheless, we retain M A as a simulation parameter to test the null result in the range 0 < M A < 11 spanned by the Iridium events.

Global simulation
The simulations were performed using the University of Michigan Space Weather Modeling Framework (SWMF) implemented at the Community Coordinated Modeling Center (CCMC, http://ccmc.gsfc.nasa.gov).The SWMF integrates various numerical models from the solar corona to the upper atmosphere into a coupled model (Toth et al., 2005).In this study we use a subset of this framework, which includes only the global magnetosphere, inner magnetosphere, and ionosphere modules.The global magnetosphere is simulated by the Block-Adaptive Tree Solar wind Roe-type Upwind Scheme (BATS-R-US) global MHD code (Powell et al., 1999).The BATS-R-US model solves the ideal MHD equations in finite volume form using numerical methods on an adaptive grid composed of rectangular blocks arranged in varying degrees of spatial refinement levels.The conditions at the simulation inner boundary, located at 3 R E geocentric distance, are determined by an ionospheric electric potential solver, which computes the two-dimensional distribution of the ionospheric electric potential from the fieldaligned currents provided by the MHD simulation using a semi-empirical model for the ionospheric conductance (Ridley and Liemohn, 2002;Ridley et al., 2004).The electric potential is mapped back to the boundary where it is used to compute the E × B flow of the magnetospheric plasma.The MHD model does not account for gradient or curvature drifts, which dominate the plasma flow of 1-200 keV ringcurrent ions, and, consequently, does not adequately simulate the pressure distribution in the inner magnetosphere.The resulting underestimate of the Region-2 Birkeland currents was evidenced in previous comparisons of observations with simulation output from the LFM MHD model (Korth et al., 2004(Korth et al., , 2008a)).In the SWMF v7.73 configuration used here, a more realistic representation of the plasma population in the inner magnetosphere is achieved through bi-directional coupling of the MHD simulation to the Rice Convection Model (RCM) (Harel et al., 1981;Toffoletto et al., 2003) as described by De Zeeuw et al. (2004).Briefly, the MHD simulation provides estimates of the magnetic and electric fields and of the plasma density, velocity, and pressure moments to the RCM.The RCM then uses a multi-fluid formalism to describe adiabatically drifting isotropic particle distributions in the self-consistently computed electric field and a magnetic field specified by the MHD model.Finally, the density and pressure spatial distributions computed from the RCM particle distribution functions are used to correct the plasma parameters in the MHD model, resulting in improved estimates of the Region-2 Birkeland currents.In the coupling framework, execution of the MHD and RCM models is alternated continuously on a simulation time scale of 10 s.The SWMF was selected for this study because initially it was the only model available at the CCMC that offered coupling of global and inner-magnetospheric models, which was found to be fundamentally important for accurate simulation of the natural system.The CCMC was chosen as the host for the simulations because of its public accessibility and because of its infrastructure to support large numbers of simulations on demand.
The parametrization of the global simulation is selected to match the range of moderate solar wind conditions corresponding to the Iridium intervals selected in Sect. 2. To keep the number of simulations required for comparison within reason, the simulations were executed for the IMF clock angles α = 0 and  315 • modifying one of the parameters E yz , p sw , or M A at a time emanating from an average magnetosphere configuration.The state of origin for the simulations is thereby given by the nominal solar wind conditions E yz = 2.8 mV m −1 , p sw = 1.7 nPa, and M A = 6.5, which correspond approximately to the median values computed from the observations in P1.The solar wind electric field was simulated for the magnitudes E yz = 1, 2, 3, and 4 mV m −1 , yielding 32 individual model runs.The solar wind dynamic pressure and Alfvén Mach number were simulated for p sw = 1, 2, 3, 4,  and 5 nPa and M A = 2, 4, 6, 8, 10, 12, and 14, respectively, resulting in 40 and 56 model runs, respectively.The total number of simulations carried out for this study is thus 128.The detailed parametrization of the model runs is contained in the supporting online materials.
Prior to simulation, the physical quantities defining the Birkeland currents must be translated into the native simulation variables, namely, the solar wind density, n sw , velocity, v sw , and the IMF components, B y and B z .Given the definitions where B yz is the magnitude of the magnetic field in the plane perpendicular to the Earth-Sun line, m p is the proton mass, and µ 0 is the vacuum permeability, the simulation variables are obtained as Two-hour simulations with fixed upstream boundary conditions were carried out for each set of modeled conditions computed from Eqs. ( 7)-( 11) on a simulation grid with ∼2 million cells, the highest resolution made available for common public use at the time of this study.The simulation requires specification of additional parameters, which were chosen as follows.First, the dipole tilt angle was set to 0 • , consistent with the normalization of the Iridium Birkeland currents to equinox conditions.Second, the temperature of the solar wind was approximated as 10 eV, or 1.2 × 10 5 K. Third, the 10.7-cm solar radio flux intensity, which is used by the semi-empirical conductance model to estimate the contribution of solar EUV ionization to the Pedersen and Hall con-ductances, was set to F 10.7 = 150 × 10 −22 Wm 2 Hz −1 .Finally, the solar wind velocity components v y,z and the IMF B x were set to 0 km s −1 and 0 nT, respectively.
To ensure a high-fidelity representation of the magnetospheric state, the simulation must be allowed to settle for a given set of boundary conditions.Since reaching the equilibrium state may require a large number of time steps, the simulation time was chosen as a compromise between the final state conditions and the computation time that can be dedicated to the simulation of an individual set of boundary conditions.To establish the simulation time, the evolution of the Birkeland total current was used as a proxy for assessing the approach of the equilibrium state.In Eq. ( 12), δA = R 2 i δφ δθ sinθ is the area element computed using R i = 6481 km for the radius of the ionosphere at an assumed altitude of 110 km, δφ = 1 h, and δθ = 2 • .Figure 1 shows the total current, I , from BATS-R-US simulations of fixed nominal solar wind conditions with IMF B z = −5 nT, v sw = 400 km s −1 , n sw = 5 cm −3 , and T = 2.3 × 10 5 K with and without coupling to the RCM as a function of simulation real time.The figure shows that the pure MHD model reaches the final state after as little as 1 h, while the time to approach equilibrium is substantially prolonged with the bi-directional coupling between the global and inner-magnetospheric models.The total current in the coupled simulation after 1, 2, 3, and 4 h amounts to 1.31, 1.43, 1.50, and 1.53 MA, respectively.At the end of the 4-h settling period, the total current of the steady-state simulation no longer increases monotonically but fluctuates by less than 1 % of the total magnitude.This uncertainty in the equilibrium state is negligible for the purpose of this study.Assuming the state at the end of the 4-h simulation interval is approximately the equilibrium state, we find that the simulation resembles the final state after 2 h to within <10 %.The error associated with limiting the simulations to 2 h of real time is typically lower than the uncertainties in the analysis of the Iridium observations discussed below.The distributions of the Birkeland current density after 1-4 h of simulated real time are shown in Fig. 2. Visual comparison of the current patterns in the figure suggests that the geometry of the Birkeland currents is also well established after 2 h.Based on the analysis above, we therefore limited all simulation runs to 2 h of real time, and used only the final state, i.e., the last simulation time step, for comparative analyses.

IMF clock angle dependence
We first compare the observed Birkeland current patterns with those obtained from the simulations.For this we proceed as described by Anderson et al. (2008)  To quantify the similarity of observed and simulated distributions, we calculated both the average fractional overlap and the two-dimensional cross-correlation coefficient for the observed and simulated current distributions within each IMF clock angle bin.The fractional overlap, Q OS , is defined as where A O and A S are the areas associated with the largescale currents in the observations and simulations, respectively, and A OS is the area of the overlap.The fractional overlap is computed separately for the upward and downward current region to yield an average magnitude, Q OS .This technique was successfully used by Anderson et al. (2008) to identify the similarity of consecutive one-hour distributions of the Birkeland currents.In the computation of Q OS , the areas in Eq. ( 13) were restricted to include only current densities larger than the mean of the 3-σ standard deviations obtained for each grid location in the Iridium statistics (∼0.02 µA m −2 ).The values of Q OS , shown for each IMF clock angle bin in Table 1, demonstrate a significant overlap, ∼60 %, for southward IMF orientations.For increasingly northward-turning IMF, the fractional overlap diminishes and reaches a minimum of 13 % for the northward direction.
The two-dimensional discrete cross correlation of the observed current density distribution j O (θ,ϕ) ≡ g k,l and the simulation result j S (θ,ϕ) ≡ h k,l is defined by where i and j are the lags applied to g k,l , G k,l and H k,l are the discrete Fourier transforms of g i,j and h i,j , the asterisk denotes complex conjugation, and F signifies the Fourier transform (Press et al., 1994).The correlation coefficient between j O (θ,ϕ) and j S (θ,ϕ) is then given by r OS = Corr(g,h) 0,0 , and the equatorward shift and azimuthal rotation of the observed distribution relative to the simulation is provided by the lags associated with the maximum correlation, r OS,max = Max Corr(g,h) i,j .The correlation coefficients, relative equatorward displacements, and azimuthal rotations are shown for all IMF clock angle bins in Table 1.
Consistent with the analysis of the fractional overlap, the correlation between observed and simulated current distributions is largest, ∼0.8, for southward IMF and smallest, ∼0.5, for northward IMF.That both the fractional overlap and the cross-correlation coefficients are lower for northward IMF conditions is predominantly due to smaller areas occupied by the large-scale currents, for which slight relative displacements yield a large penalty in these performance metrics.
From analysis of the lag in the cross correlation it is evident that the location of the large-scale current regions agree to within 2 • in latitude and that, compared to the simulation, the axis of symmetry observed in the Birkeland current distributions for southward IMF is rotated by as much as 1 h in local time toward pre-noon.Table 1.Dependence of the similarity of observed and simulated Birkeland current patterns on IMF clock angle.The similarity is quantitatively defined by the average factional overlap of the large-scale current regions, Q OS , and the two-dimensional correlation coefficient computed without rotation or shift, r OS .Azimuth rotation and equatorward shift shown are the lags maximizing the correlation coefficient.Finally, the Birkeland total current flowing between magnetosphere and ionosphere in the simulation was compared with the Iridium average distributions.The Birkeland total current was calculated from Eq. ( 12).For the observations, only current densities with magnitudes above the 3σ confidence level were included in the total current calculation.The total current for each IMF clock angle bin is shown in Table 2, the indices "O" and "S" represent the observations and simulations and the indices "down" and "up" denote the total current in downward and upward direction, respectively.Also shown in Table 2 is the ratio between the simulated and observed total current, I S /I O .

Solar wind electric field dependence
The P1 statistical analysis identified the solar wind electric field as the primary quantity determining the Birkeland current intensity.Analysis of the dependence of I on E yz is therefore the natural next step in comparing the simulation www.ann-geophys.net/29/1809/2011/Ann.Geophys., 29, 1809-1826, 2011 with the observations.Figure 5 shows scatter plots of the observed total current, I O , versus E yz (blue dots) for the stable intervals identified in the Iridium data set within each IMF clock angle bin (cf.P1).The E yz was thereby calculated from ACE observations using Eq. ( 4).The statistical dependence was estimated by linear least-squares fits of these data (blue line).Shown in red are the total currents of the simulations, I S , and the corresponding linear least-squares fits.The parameters and uncertainties for the fits are provided at the top of each panel in the respective color.The total current computed for individual 2-h intervals for nominal E yz magnitudes of 2-3 mV m −1 are noticeably larger than those obtained from the averaged Birkeland current patterns noted in Table 2, which is indicative for considerable variability in the current distribution among the events in each IMF clock angle bin.The linear fits of the simulated total currents typically reflect the lower limit to the event-base current estimates, which in some cases can exceed the simulated value by a factor of ∼5.Since the areas of the simulated and observed Birkeland currents are similar (cf.Figs. 3 and 4), this discrepancy is consistent with previous comparisons of field-aligned current densities from SWMF simulations with CHAMP observations (Wang et al., 2008).The linear fits of the event-based current estimates show disparities with respect to the simulations, which extend over the entire E yz range examined here and which consist of two components.First, the average observed total current in the magnetospheric ground state, i.e., in absence of a solar wind electric, is 2-4 times larger than obtained from the simulations.Second, for southward IMF the differences between observations and simulation widen further with increasing E yz due to the stronger dependence of the observed total current on this quantity.The observations show an intensification of the total current with increasing E yz , which is strongest for southward IMF (0.8 MA (mV m −1 ) −1 ) and gradually weakens as the IMF turns northward.The simulation reproduces the observed dependence on E yz quali-tatively, although the maximum slope in the simulation for southward IMF (0.6 MA (mV m −1 ) −1 ) is about 25 % smaller than observed.
In addition to modulating the current intensity, the analysis in P1 showed that the solar wind electric field modifies the topology of the Birkeland currents.An increase in E yz thereby results both in an equatorward expansion of the large-scale currents and a shift toward the dayside.The technique to evaluate these parameters from the distribution of the horizontal magnetic perturbations, δB, associated with the large-scale Birkeland currents was described in P1.Briefly, the grid positions of the magnetic perturbations were parameterized with a colatitude expansion factor and an offset along the noon-midnight meridian.Choosing the distribution for E yz = 2 mV m −1 as the reference δB pattern, the best-fit latitudinal expansion factor and noonward shift were then determined by minimizing the root-meansquare (rms) of the vector magnetic field residuals relative to the reference distribution.The above analysis is well-suited for southward IMF conditions, where the large-scale Birkeland current are distributed over a wide range in latitude, but typically yields unreliable results for northward IMF orientation for which the currents are concentrated at high latitudes (see P1).We therefore restrict the comparison of the simulation results with the observations to the IMF clock angle range 90 • < α < 270 • .Figure 6 shows the colatitude expansion as a function of E yz for all IMF clock angle bins with B z ≤ 0 nT in the same format as used in Fig. 5. Comparison of the observations (blue dots) with the simulation result (red dots) demonstrates that the simulation captures the equatorward expansion of the Birkeland current qualitatively.The linear fits to these data yield observed and simulated colatitude expansion factors averaging 3.0 % (mV m −1 ) −1 and 2.4 % (mV m −1 ) −1 , respectively.The observed average is elevated by the expansion factor determined for α = 180 • and E yz = 4 mV m −1 , and disregarding this data point leads to a somewhat lower slope in this clock angle bin.Overall, the equatorward expansion in the natural system is modestly (∼10-20 %) larger.The shift of the large-scale currents along the noon-midnight meridian is presented in Fig. 7.In contrast to the colatitude expansion, the P1 statistical analysis did not reveal a definite dependence of the noonward shift on IMF clock angle.Therefore, a single linear function fit was estimated from all noon-shifts evaluated in the clock angle range 90 • < α < 270 • with a given E yz magnitude.As illustrated in Fig. 7, the simulation reproduces the noonward shift of the Birkeland currents, but is, however, less pronounced.While the observations show the currents shifting duskward with increasing E yz by about 0.4 • (mV m −1 ) −1 , the shift in the simulation is only about 0.2 • (mV m −1 ) −1 and thus 50 % smaller.

Solar wind ram pressure dependence
The solar wind dynamic pressure, p sw , was identified in P1 to play a secondary role in controlling the Birkeland current intensity.To determine the dependence of I on p sw , the observed current densities must first be normalized with respect to the solar wind electric field using the relationships given in Sect. 2. For this study, we used E yz = 2.8 mV m −1 as the reference electric field, while in P1, the current densities were normalized using E yz = 2.3 mV m −1 .Consequently, the normalized total currents of the observations presented here are modestly larger than in the previous study.The solar wind dynamic pressure was computed from ACE observations using Eq.(5).To simulate different p sw conditions, the model input variables were recomputed for each p sw step using Eqs.( 7)-(11).Maintaining average conditions for E yz and M A thereby necessarily implies that the quantities n sw , v sw , and B yz vary in the simulation of different p sw values.Figure 8 shows the scatter plot of the observed and simulated Birkeland total current versus p sw sorted by IMF clock angle in the same format as used in Fig. 5.The Birkeland total current is observed to increase with increasing p sw regardless of IMF clock angle.The blue lines represent linear fits of the observed pressure dependence, showing a maximum increase for southward IMF (0.4 MA nPa −1 ) and local minima in the factors of proportionality for dawnward and duskward orientation of the IMF.The Birkeland total current increase with rising p sw is reproduced by the simulation.Similar to the observations, the simulation shows the maximum rate of increase for southward IMF.However, the rate increase is only half that observed (0.2 MA nPa −1 ).Furthermore, in contrast to the observations, the minimum increase of I with p sw was obtained for northward IMF rather than for dawnward and duskward IMF directions.

Solar wind Alfvén Mach number dependence
In P1, it was initially anticipated that the Birkeland total current also depends on the solar wind Alfvén Mach number, M A .The reason for this assumption was that M A regulates the magnetosheath plasma beta, which in turn controls the reconnection efficiency (Anderson et al., 1997).The reconnection efficiency is known to control the magnetospheric electric potential (Siscoe et al., 2002b,a), which, during nominal solar wind driving conditions, equals the ionospheric electric potential that is proportional to the Birkeland current density.However, no such dependence was detected in the range of modest solar wind conditions examined for the possible reasons discussed in P1.Here we test the simulation against the observed null result.For comparison with the simulation, the Birkeland current densities were normalized to solar wind electric field and ram pressure of E yz = 2.8 mV m −1 and p sw = 1.7 nPa, respectively.Subsequently, M A was computed for each event from ACE observations using Eq. ( 6). Figure 9 shows the scatter plot of the observed Birkeland total current (blue dots) versus M A sorted by IMF clock angle in the same format as used in Fig. 5.The linear fits of the observed dependence of I on M A exhibit slopes that are for most clock angles near or below the 2σ confidence limit of the fit.As discussed in P1, a correlation between these quantities does either not exist or is too small to be resolved by the Iridium data set.Within each panel, the total current from the simulation for E yz = 2.8 mV m −1 , p sw = 1.7 nPa, and M A = 2, 6, 10, and 14 are shown as red dots.In agreement with the observations, the dependence of I on M A in the simulation is also weak and comparable to the observations.

Discussion
The qualitative and quantitative comparisons of observed and simulated Birkeland current distributions conducted within this study have illustrated a number of similarities and differences, which are discussed below.Visual comparison of Figs. 3 and 4 shows that the simulated Birkeland current patterns are overall in good agreement with the Iridium climatology averages.Owing to consideration of drift physics in the inner magnetosphere, the simulation captures both Region-1 and Region-2 Birkeland currents for southward IMF.The simulation also reproduces the observed characteristic trans-formations in response to the IMF B y .Finally, the NBZ current system observed during northward IMF at high latitudes is reflected in the simulation output.A feature persistently observed in the Birkeland current distributions during southward IMF is a clockwise rotation of the distribution symmetry axis by about 1-h of local time with respect to the noonmidnight meridian.This bias to the reconnection geometry is not reproduced by the simulation, implying a modest difference in the convection pattern compared to the natural system.A possible reason for the discrepancy is the assumption of v y = 0 km s −1 in the simulation.Quantitative analysis of the cross-correlation between the observed and simulated Birkeland current patterns yielded an average 2-D correlation coefficient of ∼0.8, with higher correlation found for southward IMF conditions.Thus, small differences in symmetry aside, the simulation reproduces the Iridium climatology, reflecting average states of M-I coupling geometry, and its dependence on IMF orientation reasonably well.
The total currents in the simulation match those computed from the Iridium climatology generally well.The total current for northward IMF is significantly smaller than for southward IMF so that the absolute differences between observations and simulation, which are of similar magnitudes as encountered for southward IMF, carry a much larger relative weight.As a result, the observed total current for northward IMF is about a factor of 2 larger than the simulations show.Averaged over all IMF clock angle bins, the simulation underestimates the total current derived from the Iridium climatology averages by ∼10 %.The simulation thus mimics the average intensities of the M-I interaction with only modest underestimation.
The Birkeland total current computed for individual events frequently deviates from the Iridium climatology in magnitude (cf.Table 2) and furthermore exhibits significant variability with respect to the statistical average (cf.Fig. 5).The event-based and climatology-derived current estimates typically differ by a factor of ∼ 2, although the former can exceed the latter by factors up to ∼5 for some events.Because the global model was found to be consistent with the Iridium climatology, this disparity also extends to the simulation results.It should be noted that the assessment for the disparity between the Iridium observations and the model output is likely conservative because previous comparisons of fits of the Iridium magnetic perturbations with in-situ data from Ørsted and DMSP showed that the Iridium fitted peak magnetic perturbations are typically 30-50 % too low (Korth et al., 2004(Korth et al., , 2005(Korth et al., , 2008a,b),b).Similar uncertainties apply to the total current derived by integration over the Northern Hemisphere.Thus, while the M-I coupling geometry and intensity of the average magnetospheric configuration is favorably reproduced by the simulation, the static global model often falls short in its ability to predict the degree of coupling observed for a particular state.It is possible that the variability of the IMF, which is not included in the static simulations but, although reduced to the extent possible, is inherently present during the Iridium intervals identified as stable in this study, yields stronger and more variable Birkeland currents.Testing this hypothesis requires a statistical comparison of the observations with event-based simulations, which is beyond the scope of this study.
The simulation captures the trend in the dependence of the Birkeland total current on the solar wind electric field.Both observations and simulation results show the total current to increase with rising E yz (cf.Fig. 5).However, the rate of the change differs and, compared to the observations, is typically 25 % lower in the simulation.As argued in P1, the Birkeland current density, j r , is related to the solar wind drivers by during moderate solar wind conditions (see also Siscoe et al., 2002a,b).In Eq. ( 15), E sw and p sw are the magnitudes of the solar wind electric field and dynamic pressure, respectively, f r is the reconnection efficiency, and  the clock-angle dependence of magnetic reconnection, where F (0 • ) = 0 and F (180 • ) = 1.Since the Birkeland total current is the areal integral of j r , the reconnection efficiency represents the factor of proportionality between this quantity and the solar wind and IMF properties.In the natural system, the efficiency of magnetopause reconnection is controlled by the magnetic field and plasma conditions near the reconnection site (Siscoe et al., 2002b).In the simulation, reconnection is not driven by physical processes per se but instead by numerical diffusion (Raeder, 1999).Since the mechanism for magnetic reconnection in the natural system differs from that in the simulation, a disparity in the total current dependence on E yz is not unexpected.In the future, the Iridium observations can serve a useful reference for parametric studies, such as the numerical considerations examined by Ridley et al. (2010), aiming at simulating the interaction between solar wind and magnetosphere more accurately and thus improve the performance metrics of the model.
The observed total current intensification with increasing solar wind dynamic pressure is also reproduced by the simulation.Although the rate of increase in the total current differs between the simulation and the observations by about a factor of 2, both data sets contrast the decrease in the total current with increasing p sw predicted by Eq. ( 15).It was suggested in P1 that a portion of the Birkeland currents is generated by viscous interaction and this contribution is not included in Eq. ( 15).This interpretation was supported by the fact that the total current dependence on p sw is weakest for dawnward and duskward IMF orientations, for which the magnetosphere is least susceptible to the Kelvin-Helmholtz instability, which presumably drives the viscous dynamo (Axford and Hines, 1961).However, this interpretation does not apply to the simulation.Twodimensional MHD simulations of the Kelvin-Helmholtz instability by Nykyri and Otto (2001) have shown that the spatial scales of the resulting plasma vortices are of the order of hundreds of kilometers.Since these structures are much smaller than the simulation grid resolution at the flanks of the magnetopause, which ranges from 0.25 R E at the dawn-dusk terminator to 2 R E at locations 20 R E downtail, the Kelvin-Helmholtz instability cannot be resolved by the global MHD model.Furthermore, minima in the total current dependence on p sw are obtained for northward and southward IMF, inconsistent with Kelvin-Helmholtz waves, which maximize under these IMF orientations.Thus, the total current intensification with increasing p sw seen in the simulation results cannot be attributed to viscous interaction contributing to the generation of the Birkeland currents.It is possible that pressure enhancements lead to higher stresses in the coupling between magnetosphere and ionosphere, which are mediated by an increased flow of Birkeland currents.In specific, as argued by Palmroth et al. (2004), the plasma sheet pressure is in balance with the solar wind pressure (see also Borovsky et al., 1998), so that variations in solar wind pressure indirectly modify the Region-2 Birkeland current densities as j ,R2 ∝ ∇V × ∇P , where V and P are the flux tube volume and inner-magnetospheric pressure (Vasyliunas, 1970), respectively.Such process is not considered in the Hill model of the transpolar potential represented by Eq. ( 15).Therefore, deviations of the observed and simulated pressure dependencies of the total current from predictions of the simplistic Hill model are conceivable due to incomplete implementation of the physical processes of the natural system.
The model's ability to reproduce the natural system depends critically on the implementation details of individual model components and of the coupling between these modules.While field-aligned currents are computed in both the inner-magnetosphere and global MHD modules, only currents from the MHD module are used as input parameters to the ionospheric electrodynamics module.In the MHD model, field-aligned currents are generated locally where curls exist in the magnetic field and flow along magnetic field lines to the simulation inner boundary.If the grid cells in the simulation are not small enough, the simulated currents can diffuse and close prematurely in the magnetosphere (Ridley et al., 2010).Furthermore, Region-2 currents in the MHD model are driven by pressure gradients generated by the inner-magnetosphere model.Assuming that the numerical dissipation in the inner-magnetosphere and MHD codes is similar, the field-aligned currents in the inner magnetospheric and MHD models are equivalent.However, this coupling scheme breaks down when pressure gradients in the inner magnetosphere extend to radial distances below the simulation inner boundary, because the magnetic field distortions associated with the diamagnetic currents cannot be adequately modeled.Above, only two of many ways in which the model implementation affects the simulation results are mentioned.In this study, we have identified both similarities and differences between model and observations, and it is hoped that our results prompt investigations into the origin of the latter.

Summary and conclusions
We have compared statistical dependencies of the Birkeland currents on solar wind parameters, which were derived from Iridium distributed magnetic perturbation data and discussed by Korth et al. (2010), with a global MHD simulation results.The simulations were conducted at the publicly-accessible Community Coordinated Modeling Center using the University of Michigan Space Weather Modeling Framework with coupling to the Rice Convection model.The key findings from the comparison are: 1.The global model reproduces the Iridium statistical Birkeland current patterns well, showing a twodimensional correlation coefficient of ∼0.8 averaged of all IMF clock angles bins.The simulation modestly underestimates the total current intensity of the Iridium climatology by <10 % on average.
2. The current intensity for individual events may deviate significantly from the climatology averages.Consequently, it is not uncommon for the observed total current to exceed this quantity in the simulation by factors of ≥2.Considering that the Iridium-based estimates are up to 50 % lower than the true total current, our assessment of the disparity is likely conservative.
3. The Birkeland total current increase with rising solar wind electric field magnitude and ram pressure is qualitatively reproduced by the model, but the rate of increase is observed to be twice that reported by the simulation.
4. The observed equatorward expansion and noonward shift of the Birkeland currents with increasing solar wind electric field are also evident in the simulation results but are not as pronounced.
The purpose of this study is to present an example for a comprehensive and objective assessment of global magnetosphere simulation tools with focus on a particular set of coupled models.Because the models and their performance evolve in time, this analysis can only represent a snapshot in the assessment of the state of our modeling capabilities.

Fig. 1 .
Fig. 1.Evolution of the Birkeland total current during a 4-h simulation of nominal southward IMF conditions using the SWMF model with (solid line) and without (dotted line) coupling to the RCM.
and sort the 1536 two-hour distributions observed by Iridium into 45 •wide bins of the IMF clock angle centered at 0 • , 45• , 90• , 135 • , 180• , 225• , 270• , and 315 • .The averaged current density distributions within each clock angle bin are shown in Fig.3, where upward and downward currents are shown in red and blue, respectively, and the IMF direction for each distribution is indicated by the arrows in the center of the panel.Only current densities above the 2σ confidence level are shown, and black contours are placed in 0.1 µA m −2 increments to emphasize the structure in the distributions.For comparison, the current distributions obtained from the simulations are shown in Fig.4using the same format as used in Fig.3.The Iridium average distributions include observations from the full range of solar wind electric fields included in the database, i.e., 0 ≤ E yz ≤ 4.5 mV m −1 .For Fig.4, we therefore averaged the simulation results for E yz = 1, 2, 3, and 4 mV m −1 within each clock angle bin.To facilitate the comparison between the distributions of Figs.3 and 4, regions where upward and downward currents are in common between the observations and the simulations are outlined by solid red and blue contours, respectively.From visual comparison of these contours we find that the simulated Birkeland current patterns are overall in good agreement with the Iridium average distributions.

AnnFig. 2 .Fig. 3 .Fig. 4 .
Fig. 2. Distributions of the upward (red) and downward (blue) Birkeland current density from the SWMF model with RCM coupling (a) one, (b) two, (c) three, and (d) four hours of real time into the simulation of nominal southward IMF conditions.

Fig. 5 .
Fig. 5. Dependence of total current on solar wind electric field, E yz , and IMF clock angle, α, for the Iridium observations (blue dots) and the BATS-R-US plus RCM simulations (red dots).The lines of the respective color represent the linear least-squares fit functions noted at the top of each graph.

Fig. 6 .
Fig. 6.Dependence of the colatitude expansion of the large-scale Birkeland current distribution on solar wind electric field, E yz , and IMF clock angle, α, for the Iridium observations (blue dots) and the BATS-R-US plus RCM simulations (red dots).The lines of the respective color represent the linear least-squares fit functions noted at the top of each graph.

Fig. 7 .
Fig. 7. Dependence of the noon-midnight shift of the large-scale Birkeland current distribution on solar wind electric field, E yz , for the IMF clock angle range 90 • < α < 270 • .The Iridium observations and the BATS-R-US plus RCM simulations are represented by blue and red dots, respectively.The lines of the respective color represent the linear least-squares fit functions noted at the top of each graph.

Fig. 8 .
Fig. 8. Dependence of total current on solar wind ram pressure, p sw , and IMF clock angle, α, for the Iridium observations (blue dots) and the BATS-R-US plus RCM simulations (red dots).The lines of the respective color represent the linear least-squares fit functions noted at the top of each graph.

Fig. 9 .
Fig. 9. Dependence of total current on solar wind Alfvén Mach number, M A , and IMF clock angle, α, for the Iridium observations (blue dots) and the BATS-R-US plus RCM simulations (red dots).The lines of the respective color represent the linear least-squares fit functions noted at the top of each graph.

Table 2 .
Dependence of the Birkeland total current I on IMF clock angle computed from the observed (index "O") and simulated (index "S") Birkeland current distributions.The average total currents, I O and I S , represent the averages of the magnitudes of the respective total downward (index "down") and total upward (index "up") currents.

material related to this article is available online at: http://www.ann-geophys.net/29/1809/2011/ angeo-29-1809-2011-supplement.pdf.
Global Birkeland current distributions of much higher fidelity than used for the present study are now publicly available through the Active Magnetosphere and Planetary Electrodynamics Response Experiment (AMPERE).These data should find wide-spread use within the research community in the future, among others, to facilitate continued critical assessment of simulations tools.