The impact on global magnetohydrodynamic simulations from varying initialisation methods : results from GUMICS-4

We investigate the effects of different initialisation methods of the GUMICS-4 global magnetohydrodynamic (MHD) simulation to the dynamics in different parts of the Earth’s magnetosphere and hence compare five 12 h simulation runs that were initiated by 3 h of synthetic data and followed by 9 h of solar wind measurements using the OMNI data as input. As a reference, we use a simulation run that includes nearly 60 h of OMNI data as input prior to the 9 h interval examined with different initialisations. The selected interval is a high-speed stream event during a 10-day interval (12–22 June 2007). The synthetic initialisations include stepwise, linear and sinusoidal functions of the interplanetary magnetic field with constant density and velocity values. The results show that the solutions converge within 1 h to give a good agreement in both the bow shock and the magnetopause position. However, the different initialisation methods lead to local differences which should be taken into consideration when comparing model results to satellite measurements.


Introduction
The solar wind plasma density, velocity and temperature as well as the orientation of its embedded interplanetary magnetic field (IMF) are fundamental in determining the efficiency of solar wind-magnetosphere-ionosphere coupling processes.The plasma mass and momentum transport into the magnetosphere are driven by local plasma processes, which have complex dependencies on the large-scale configuration.Magnetic reconnection (Dungey, 1961) close to the sub-solar point accounts for majority of plasma transport during southward IMF, while viscous processes (Axford and Hines, 1961) such as the Kelvin-Helmholtz instability (Nykyri and Otto, 2001) or diffusion (Johnson and Cheng, 1997) may play a strong role, particularly when the IMF is directed northward (Osmane et al., 2015).Quantification of the relative importance of these processes under varying solar wind conditions is a central element to developing capabilities to predict space weather in the geospace.
Several overlapping and multi-spacecraft missions, such as Cluster (Escoubet et al., 2001) and MMS (Burch et al., 2016), are operational today, offering statistically global coverage via their conjunction.However, it is not feasible to use these for real-time space weather prediction.Instead, global magnetohydrodynamic (MHD) models are often used to provide a large-scale view of the dynamical evolution in real-time.The best known models include the Lyon-Fedder-Mobarry (LFM) (Lyon et al., 2004), BATS-R-US (Powell et al., 1999), Open General Geospace Circulation Model (OpenGGCM) (Raeder et al., 2009), Ogino Model (Ogino et al., 1994) and the Grand Unified Magnetosphere-Ionosphere Coupling Simulation (GUMICS) (Janhunen et al., 2012).
Global MHD simulations typically start with a dipole in an almost empty space.In the initialisation process, the solar wind driver is turned on to create the upstream bow shock and magnetopause structure as well as the downstream magnetotail and plasma sheet configuration.Typically, con-Published by Copernicus Publications on behalf of the European Geosciences Union.
stant or slowly varying parameter values are used during the initialisation period (Bhattarai et al., 2012;Hoilijoki et al., 2014), and its duration is dictated by the time needed (of the order of 1 h) to form the magnetosphere (Janhunen, 1996;Raeder et al., 2009) and reset the magnetospheric memory (a few hours) of past dynamical processes (Büchner et al., 2003, pp. 212-246).
While the different MHD simulations are based on the same plasma theory, the exact form of the equations, the numerical solutions, and the initial and boundary conditions as well as the ionospheric coupling solutions vary.Several studies have focused on assessing the performance of the models through comparisons of the simulation results with in situ or remote observations of dynamic events or plasma processes (Birn et al., 2001;Pulkkinen et al., 2011;Honkonen et al., 2013).Recent comprehensive studies (Juusola et al., 2014;Gordeev et al., 2015) statistically quantify the strengths and weaknesses of each of the models, while they could not identify any of the codes as clearly superior to the others.
As a part of validation process the simulation results are often compared with point measurements made by satellites.As the magnetosphere includes regions with large spatial and temporal gradients, even small errors in the simulation configuration can cause large differences with respect to the observations locally at a single point.In such a case, it is difficult to interpret how well the simulation reproduces the largescale dynamic sequence.Moreover, local and large-scale effects of initialisation process are not known.In order to quantify variations that may arise from the initialisation procedure (history of the simulation), we investigate the effect of artificial initialisation period containing stepwise, linearly or sinusoidally changing IMF on a GUMICS-4 global MHD simulation run over a 9 h period following the initialisation.The quantitative differences are examined both locally considering difference in convergence of several variables at grid points located in different parts of the magnetosphere and globally studying differences in momentum transport over the magnetosphere.
This paper is constructed in the following way.Section 2 describes the features of GUMICS-4, and Sect. 3 characterises the chosen interval and discusses the ability of GUMICS-4 to reproduce the large-scale magnetospheric boundary structure.Section 4.1 illustrates the initialisation methods.The model solutions and their differences obtained from the multiple initialisation procedures are presented in Sect.4.2-4.4.The paper ends with discussion and conclusions.

The GUMICS-4 global MHD model
The fourth edition of the Grand Unified Magnetosphere-Ionosphere Coupling Simulation (GUMICS-4) is a global 3-D MHD simulation code for the magnetosphere and the solar wind coupled to an electrostatic spherical ionosphere region (Janhunen et al., 2012).The MHD solver utilises the finitevolume method.The simulation box has dimensions of +32 to −224 R E in the x GSE direction, and −64 to +64 R E in the y GSE and z GSE directions.The inner boundary is spherical with a radius of 3.7 R E , with a gap separating the magnetospheric and ionospheric simulation regions.The magnetosphere provides the ionosphere with the field-aligned current pattern and the electron precipitation, while the ionosphere gives the electric potential to the magnetosphere.This feedback loop is updated every 4 s.
GUMICS-4 solves the ideal MHD equations with the separation of the magnetic field to a curl-free (dipole) component and perturbed component created by currents external to the Earth (B = B 0 + B 1 (t); Tanaka, 1994).Two important features of GUMICS-4 are its adaptive grid and temporal subcycling.Both make the computations feasible on a single PC.An adaptive grid enables the coupling of a high-resolution ionospheric grid with coarse magnetospheric grid to execute computations in a reasonable time.To limit the number of grid cells in the MHD region, the grid cell size is varied depending on the local gradients (the plasma mass density, the total energy density without the contribution of the dipole field, the momentum density and the perturbed component of the magnetic field).Temporal subcycling reduces number of MHD solutions an order of magnitude while guaranteeing that the local Courant-Friedrichs-Lewy (CFL) constraint (Lions and Ciarlet, 2000, pp. 121-151) is satisfied.
The input to the GUMICS-4 are solar wind plasma number density, temperature, velocity and the IMF conditions at the upstream boundary X = 32 R E .The smallest grid cell size in this study was 0.5 R E .While it would be desirable to use the maximum resolution of 0.25 R E , due to the long simulation physical time (10 days), this was not feasible, nor was it necessary for the purpose of this study.Geocentric solar ecliptic (GSE) coordinates are used throughout this paper.

Model inputs
We use solar wind parameters from the OMNI database at 1 min resolution as inputs to GUMICS-4.OMNI data are a composite of solar wind measurements typically measured at the Lagrangian point L1 and then subsequently propagated (King and Papitashvili, 2005) to the bow shock nose using the Farris and Russell (1994) bow shock model.OMNI data were obtained through the NASA OMNIWeb service (http: //omniweb.gsfc.nasa.gov).
3 Simulation of a high-speed stream event

Event description
Continuous solar wind measurements are a necessary precondition for obtaining a reliable simulation result of the magnetospheric dynamics.While such periods are surprisingly uncommon, the period from 12 June, 04:00 UT, to 22 June, 04:00 UT, 2007 contains a mostly uninterrupted solar wind record (from the Wind satellite during 12 June, 04:00 UT-13 June, 20:59 UT, 2007, and from the ACE spacecraft thereafter) in the OMNI database with only short data gaps.These gaps make up 14 % of the period.We have simulated the full 10-day period that comprises a high-speed solar wind stream interval, and we use the simulation results as a reference for further numerical experiments.
Figure 1 shows the IMF B Y and B Z components, solar wind speed, and density.A high-speed stream was initiated during 14 June 2007, as indicated by the enhanced speed (∼ 650 km s −1 ) around 00:00 UT on 15 June 2007.The leading edge was associated with a B Y reversal from slightly negative to positive and rapidly fluctuating IMF B Z .The speed gradually decreased toward the end of the interval from the peak value to about 400 km s −1 .During the high-speed stream, the density was low at a few particles per cubic centimetre.The IMF B Y was positive at a level of a few nT, and the IMF B Z fluctuated around zero; 21 June features rapid density increase with IMF B Y and IMF B Z fluctuating between −10 and 10 nT.We focus our attention particularly to the period when the solar wind speed increases during 14 June, 12:00-24:00 UT, to cover both negative and positive IMF B Z and a range of solar wind speed values.
In the magnetosphere, the high-speed stream did not cause major activity as the IMF B Z did not include periods of intense southward fields after the initial sheath region on 13 June (not shown).The Dst index remained above −20 nT throughout the period.The sheath region leading to the highspeed stream included a period of AE activity that exceeded 500 nT, but during the high-speed stream the AE activity was quite low, with only moderate substorms until the end of the stream, which again included a period of higher AE activity during 21 June 2007.Thus, from the magnetospheric activity point of view, the period can be characterised as driven by a fast solar wind stream, rather fluctuating IMF driver, which resulted in only moderate magnetospheric and ionospheric responses.
The 1 min OMNI data were interpolated across short data gaps to obtain a continuous simulation input.One hour of fixed solar wind values was included at the start of the interval (12 June 2007 03:00-04:00 UT).In order to preserve ∇ • B = 0, the IMF B X was set to zero.

Bow shock and magnetopause locations
The location of the bow shock upstream of the magnetosphere is examined utilising the Rankine-Hugoniot conditions (Fitzpatrick, 2014).Assuming that the plasma upstream and downstream of the bow shock is spatially uniform, we find that the shock is time-stationary and planar, and its front is aligned with the Y Z plane; the shock can be located by examining the solar wind velocity and magnetic field in the XY plane.
Rankine-Hugoniot relations make the assumption that the shock is a thin layer, as it is assumed that dissipation takes place at collisional kinetic scales below the fluid representation of the plasma.The plasma fluxes (mass, momentum, energy) are conserved, resulting in jumps of the moments across the shock transition layer.However, the layer can only be located up to the smallest scale in the GUMICS-4 simulation grid, which in this study is 0.5 R E in the vicinity of the shock.Thus, the location of the shock layer can be determined to an accuracy of ±0.5 R E .To trace the changes, we cover sufficiently large distance away from the shock region upstream and to beyond the shock front in the downstream direction ending before reaching the magnetopause.
We concentrate on locating the shock distance at the Sun-Earth line from the Rankine-Hugoniot relations.By letting both the upstream and downstream plasma flows to be perpendicular to the magnetic field the MHD equations in the rest frame of the shock yield a perpendicular shock: where B is the magnetic field, ρ is the plasma density, and V is the plasma speed.Subscripts u and d refer to regions upstream and downstream of the shock, respectively, and the ratio r is given by where p is the plasma pressure, γ = 5/3 is the ratio of the specific heats, and M u the upstream sonic Mach number.Equation (1) provides three parameters (magnetic field, density and flow velocity) for identifying the bow shock standoff distance.
To compute the shock position from the simulation sufficiently far from the Earth averaged over several grid points, we first evaluate the sonic Mach number (Eq. 2) by using values of V , P and ρ averaged over a 4 R E distance upstream of the shock between X = 24 and 28 R E , which gives us the ratio r relating the upstream and downstream values in Eq. (1).For example, plasma flow speed of 493 km s −1 results in r ≈ 3.92, thus suggesting The values of B, ρ and V are evaluated along the Sun-Earth line towards the Earth until they reach the predicted downstream value; the bow shock position is defined as the first grid point that satisfies Eq. ( 1) for each parameter individually providing three estimates of the shock position.Note that this method identifies the inner edge of the shock, as we trace the first point inside of the shock layer.Due to 0.5 R E maximum resolution the accuracy of each shock position estimate is less than or equal to 0.5 R E .This procedure is repeated at 1 min temporal resolution throughout the simulation run, but for slight smoothing we use 10 min sliding averages both for the simulation results and the OMNI data.The standard deviation of each 10 min ensemble is used as an uncertainty estimate.
To provide a reference to the three shock position estimates and to be consistent with the used OMNI solar wind data, the estimates are compared with empirical values from the OMNI database derived from the Farris and Russell (1994) bow shock model.A comparison of the shock location between OMNI and GUMICS-4 is shown in panel (a) of Fig. 2 together with the relative difference between the shock position obtained using the three variables B, ρ and V (panel b).The standard deviations over 10 min intervals used as a proxy for uncertainty are shown by vertical bars.The relative difference in the bow shock standoff distance between GUMICS-4 and OMNI is mostly positive with 88 % of the values in the range of 0-20 % for B and over 99 % for both ρ and V , thus suggesting that using ρ and V results as lower difference compared to OMNI.Panel (b) in Fig. 2 shows that the blue plot (V ) is consistently closer to zero line than green (ρ) and especially magenta (B), which is confirmed by calculated average relative differences over the 10-day period that are 10.8, 10.7 and 6.7 % for B, ρ and V , respectively.Regardless of the chosen variable, the applied method consistently predicts the bow shock position closer to the Earth by up to several R E than that obtained from OMNI.The general trend is expected given the 0.5 R E spatial resolution, which leads to a wider making MHD shock in the simulation than in reality.Furthermore, global MHD models generally underestimate the ring current in the inner magnetosphere, which positions the magnetopause and hence the shock closer to the Earth (Janhunen et al., 2012), especially during times of strong geomagnetic activity.However, GUMICS-4 predicts the temporal evolution of the bow shock location quite accurately over timescales of tens of minutes.This holds even when abrupt changes in solar wind conditions take place (see number density values on 13, 14 and 21 June 2007 in Fig. 1) as such effects are transient in nature.
The magnetopause subsolar point position given by GUMICS-4 is compared to the empirical Shue et al. (1997) model, which along the Sun-Earth line is reduced to the form where r 0 is the magnetopause standoff distance at the subsolar point, p dyn is the solar wind dynamic pressure in nPa, a = 11.4,b = 0.013 for B Z ≥ 0 and b = 0.14 for B Z < 0. The position of the subsolar magnetopause in GUMICS-4 is computed using the current density J Y component.Starting from the bow shock position obtained from Eq. ( 1) and ending at X = 5 R E , we identify the maximum value of J Y as the magnetopause position.Ten-minute averages similar to the bow shock location determination are used both in evaluating the Shue model and the GUMICS-4 magnetopause position.
Figure 2 shows the comparison of the magnetopause subsolar point location (panel c) and the relative difference between the position obtained using GUMICS-4 J Y and the Shue model (panel d).The general trend is similar to the bow shock position: GUMICS-4 magnetopause position is closer to the Earth than the predicted by the Shue model, in good agreement with previous studies (Palmroth et al., 2001), and results mostly from the ring current underestimation.However, the relative difference is smaller and fluctuates between −10 and 20 % with 98.5 % of the data in the range of 0-10 % and 77.4 % of the data in the range of 0-5 % during the 10-day period.Average relative difference over the 10-day period is 3.4 %.The relative difference is more enhanced when the number density and the solar wind speed (Fig. 1) increases and thus dynamic pressure increases and IMF B Z fluctuates with an amplitude up to 10 nT (14 and 21 June 2007).

Initialisation sequences
In order to study the effects of simulation initialisation, we ran a set of five simulations during the 10-day period, starting at 12:00 UT on 14 June 2007.Input solar wind data in each simulation consists of 9 h identical subset of the 10-day run input data section (15:00-24:00 UT) preceded by varying 3 h synthetic initialisation period (12:00-15:00 UT).The 3 h length was chosen to distinctly exceed 1 h time needed for the magnetosphere to form (Janhunen, 1996).To quantify the effect of the different initialisation procedures in GUMICS-4, we compare the model solutions obtained from each initialisation procedure by examining time series of several variables in different points across the magnetosphere (Sect.4.2) and surface plots of the magnitude of the relative difference in ρV X variable between the 10-day reference run and the runs using different initialisations (Sect.4.3).
We performed five simulations with constant solar wind parameters (n, T , V X , V Y , V Z ) but variable IMF during the initialisation periods.The fixed values for density and velocity are chosen to be typical values in the solar wind, while the value for temperature is an average over the 10-day period.Two simulations employ a stepwise-varying B Z with constant B Y , two linearly varying B Z with constant B Y , and one simulation uses a sinusoidally varying B Z together with a sinusoidally varying but anti-phase (90 • phase difference) B Y .The chosen initialisation methods both utilise simple functions to determine the time evolution of IMF B Z and are used in previous studies (see e.g.Bhattarai et al., 2012;Palmroth et al., 2001).As the IMF components do not change spatially at the solar wind inflow boundary (XZ plane, X = 32 R E ), the magnetic field B X component is fixed (zero in the present study) to preserve the ∇ • B = 0 condition at the boundary.The IMF B Y component is fixed to 2 nT for the four first sim- The initialisations are shown with different colours (st1 in blue, st2 in blue with symbols, lin1 in magenta, lin2 in magenta with symbols and sin in green, 10-day (10d) run in black).The 10-day run data are from OMNI.The yellow background highlights the 1 h interval during which the magnetosphere should be formed (Janhunen, 1996;Raeder et al., 2009).
ulations.The IMF B Z component is varied between −5 and 5 nT, values typically observed in the solar wind (Dimmock et al., 2014).The first stepwise-varying initialisation starts from the B Z maximum, while the second one uses the minimum as a starting point.The same applies to the two linearly varying initialisations.Sinusoidally varying B Z is initially at 5 nT, combined with sinusoidally varying B Y starting from 0 nT. Figure 3 illustrates the initialisation conditions covering the 3 h initialisation period and the first hour with the observed solar wind data.From this point, we refer to the different initialisations as st1 (stepwise varying starting from +5 nT), st2 (stepwise varying starting from −5 nT), lin1 (linearly varying starting from +5 nT), lin2 (linearly varying starting from −5 nT), and sin (sinusoidally varying starting from +5 nT).Tables 1 and 2 summarise the numerical values of the solar wind and IMF parameters.
The first 3 h of the GUMICS-4 simulation were run using the initialisations described above.Following the 3 h period, the simulation was run using the observed solar wind parameters similarly to the 10-day run.This produces five 9 h simulation periods that we compare with the 10-day simulation that starts from 04:00 UT on 12 June 2007.The period was Table 1.Fixed solar wind parameters with the associated numerical values.Note that B Y was not fixed in sinusoidally varying B Z method.
Fixed parameters Numerical value chosen to comprise abrupt changes in the solar wind parameters that result in dynamic changes in the magnetosphere (Fig. 1); not only do the polarities of B Y and B Z change multiple times but also the plasma speed undergoes a rapid change from 300 to 600 km s −1 .

Differences in lobe and plasma sheet dynamics
In order to quantify the simulation differences, we have selected points in the plasma sheet and in the lobe at four locations along the magnetotail and one point at the dayside magnetosheath.By using these points we cover three important regions in the magnetosphere.Choosing identical X coordinate values for both the plasma sheet and the lobe points allows us to compare the conditions in the two regions approximately equal distance away from the Earth in the X direction.On the other hand, selecting a point in the magnetosheath several Earth radii away from the Sun-Earth line can actually lead to studying a point that is constantly moving between the solar wind and the magnetosphere and thus contribute to revealing how GUMICS-4 reproduces the bow shock and the magnetopause positions using a different approach than what was described in Sect.3.2.At these locations, we produce absolute difference time series of several parameters.For the plasma sheet and the lobe, these parameters are V X , n and B Y , while in the magnetosheath we examine time series V , B, plasma beta β, and E Y .The selected points are shown in Fig. 4.
Figure 5 shows the absolute difference time series in the tail at X = −8 and −15 R E (panels a2, a3, a4, b2, b3, b4, c2, c3, c4, d2, d3, d4).Panels (a1, b1, c1, d1) illustrate plasma flow speed at the selected points.Time series plots at X = −10 and −20 R E are found in Fig. A1 in Appendix A. The general trend is that, by 16:00 UT, the convergence has taken place almost everywhere for the three variables, thus being consistent with earlier studies (see e.g.Janhunen, 1996) which state that the formation of the magnetosphere takes approximately 1 h.Some interesting features remain, though: while generally V X shows similar time evolution for each simulation, there are significant differences (several hundreds of km s −1 ) in the lobe close to the Earth at −8 R E lasting several hours (Fig. 5).The largest differences are produced by the sin initialisation.The differences decay gradually and eventually diminish by 17:00 UT (23:00 UT) at X = −8 R E .Further away at X = −15 R E the differences in V X get as high as 50-100 km s −1 and vanish by 23:00 UT.Remarkably low number density values and the fact that similar effects are not observed in the plasma sheet, where the number density is higher suggest that differences on the large-scale mass transport in the magnetotail are small.Z sheet refers to varying Z coordinate in the plasma sheet.In each panel, subpanels from top to bottom: plasma flow speed in km s −1 , absolute difference in plasma velocity V X in km s −1 , density in cm −3 and B Y in nT, all in GSE coordinates.Absolute difference is calculated as 10d − x, where x = st1, st2, lin1, lin2, sin.
To further test the origin of the differences, we performed yet another simulation with sinusoidally varying B Z but constant B Y .This led to substantially smaller differences in the lobe plasma transport values (not shown).Thus, the differences arise from the configurational changes associated with IMF B Y due to the significant role of B Y in tail rotation.In GUMICS-4, the magnetotail B Y is usually approximately 50 % of the solar wind value (Janhunen et al., 2012), which roughly corresponds to observations (Sergeev, 1987).This result highlights the importance of using a B Y that is close to the observed value in the initialisation process, as the effects are clearly visible in the simulation results for several hours after the initialisation period.
Note that since the magnetopause and the bow shock are in motion as depicted in panel (a), this point is not necessarily always located within the magnetosheath but fluctuates relative to the magnetopause position crossing to the upstream solar wind when the magnetopause moves inward.This can be seen in the time series as small values of the magnetic field magnitude in panel (c) that often coincide with the bow shock and the magnetopause moving closer to the Earth in panel (a).The time series measured at X, Y, Z = [9, 6, 0] show only minor differences between the runs following the first hour of similar solar wind input, in-   GSE coordinates).Yellow background highlights 1 h interval during which the magnetosphere should be formed (Janhunen, 1996;Raeder et al., 2009).
dicating no locally introduced dynamic evolution that would carry the memory of the previous driving conditions.

Differences in the noon-midnight meridian plane
To gain a large-scale view of the simulation differences, we focus on mass transport X component ρV X in the XZ plane.
Figure 7 shows time instance 18:30 UT on 14 June 2007, which represents the time interval when the highest differences occurred over the course of the 12 h simulation period.The top left panel (a) in Fig. 7 shows the ρV X scaled by the upstream value at X = 14.5 R E for the 10-day run and reveals the large-scale characteristics of the solar wind, magnetosheath, and magnetosphere in the form of colour-coded rel-ative plasma transport (scaling not shown).The black sphere in the middle of each panel (3.7 R E radius) masks the region not covered by the MHD simulation.The white sphere inside the black sphere depicts the Earth.
The following panels (b-f) show the magnitude of the relative difference in ρV X between the 10-day run and st1, st2, lin1, lin2 and sin, respectively, with the colour scale shown indicating the differences as percentages.The pattern of the large relative differences approaching and in excess of 100 % in the inner part of the magnetosphere (shown with red colours) is qualitatively similar in each of the panels (bf) in Fig. 7 (18:30 UT).These differences are mainly caused by differences in V X .The localised nature of the differences indicates that these arise at the boundaries of plasma regions, where a slightly different orientation of the overall magnetospheric configuration can lead to large differences in the local values of the plasma parameters.While such effects may not affect the large-scale evolution of the magnetotail dynamics (as demonstrated by the convergence of the time series above), they can cause large differences at local point measurements.Such differences are important when simulation results are compared with spacecraft observations -just a slight change in the satellite orbit might lead to significantly different time series along the spacecraft trajectory.
Appendix Figs.A2-A4 show time instances 16:03, 17:12 and 19:30 UT on 14 June 2007 in a format similar to Fig. 7.The overall pattern of the large relative differences in the inner part of the magnetosphere is the same in all four figures.However, panel (e) in Fig. A2 shows a bright region in the nightside tail lobe extending from −15 to −25 R E in the X direction, indicating an absolute difference in the range of 30 to 40 %.It appears to be an extension to one of the bright regions in the inner tail, and is mainly caused by the decrease in the tailward flow speed V X in the lin2 run with respect to the 10-day run.Such differences imply differences in the overall magnetotail configuration and dynamic evolution, for example, generation of flow bursts in the magnetotail.As the differences move around the noon midnight meridian plane and do not grow as a function of time, the results indicate that the simulations converge towards similar long-term dynamical evolution.

Quantification of the differences
Above, we have shown that there are large, localised differences and smaller, wider-spread differences in the plasma transport properties.In order to quantify the overall convergence of the different initialisations with the 10-day reference run, we compute kernel density functions of the differences integrated over all grid cells and over 4 h covering the time period 16:00-20:00 UT.This time period includes all four instances studied in Figs.7 and A2-A4.The grid used is evenly spaced, interpolated from the original non-uniform GUMICS-4 grid, including 13 015 grid points at 0.5 R E resolution, covering the spatial volume of −25 < X < 29 R E , −30 < Z < 30 R E .
Figure 8 shows the kernel density functions for the differences for each of the initialisations.It is clear that for most grid points the differences are very close to zero, below 5 %.Such points show as dark-blue regions in Figs.7 and A2-A4.Furthermore, it is also evident that only a small minority of grid points show relative differences over 5 % which cover the large localised errors (extreme values shown in red) and the region in between (error between 5 and 30 %), which includes the smaller errors over larger regions (shown in light blue) in Figs.7 and A2-A4.
The very thin tails of the distribution functions with over 5 % errors demonstrates that the large and medium-sized errors are statistically insignificant and do not lead to largescale differences in the temporal evolution of the simulations.
Comparing the kernel density functions in Fig. 8 shows that the performance using the different initialisations shows some differences, but overall all initialisations work reasonably well.However, two simulations starting with positive IMF B Z and utilising constant IMF B Y (st1, lin1) seem to produce fewer errors in the below 5 % range (lower peak values) and thus produce larger errors in the 5-30 and > 30 % ranges, indicating some differences in the magnetotail dynamics.The third simulation starting with positive IMF B Z (sin) differs from the two in terms of peak value in Fig. 8.However, as was shown in Sect.4.2, it produced largest localised differences partly because of the usage of fluctuating IMF B Y .This would indicate that it is useful to start the initialisation with strongly reconnecting magnetosphere driven by negative IMF B Z and constant (close to observed value) or near-zero IMF B Y .
We classify the simulations by computing probabilities for having 0-5, 5-30 and > 30 % differences during two time intervals, 16:00-20:00 and 16:00-24:00 UT, on 14 June 2007.Table 3 shows the probabilities for having the relative difference in ρV X in the aforementioned percentage ranges during the indicated time intervals.The first set of probabilities (16:00-20:00 UT) contains the same data as Fig. 8.As the integration times are much longer than the solar wind transit time through the simulation box, the differences in the numbers between the two integration times give an indication of the statistical errors (a few percent in each direction).
About 5-7 % of the time the relative errors fall between 5 and 30 %, indicating that there are instances when the largescale patterns are slightly different, but these do not last for extended periods of time.The small localised errors are of the order of less than 1 %, and do not generate larger errors as the simulation continues.Based on these statistics, it is not easy to pick one initialisation method over another.However, we would slightly prefer st2 and lin2 initialisations Kernel density . Kernel density functions of the magnitude of the relative difference in ρV x for st1, st2, lin1, lin2 and sin runs integrated over 4 h from 16:00 to 20:00, 14 June.The kernel density values are normalised by the total grid volume.

Discussion
In this paper we study the performance of the GUMICS-4 global MHD simulation.First, we compare the bow shock standoff distance and magnetopause subsolar point position with empirical models using the OMNI database containing the empirical shock model by Farris and Russell (1994) and the Shue et al. (1997) empirical model for the magnetopause location.The performance was evaluated using a simulation covering a 10-day period during 12-22 June 2007, when nearly continuous solar wind observations were available from the OMNI database.We conclude that depending on the method used (B, ρ, V ) the GUMICS-4 shock position stays within 20 % of the empirical value 88 or 99 % of the inspected time interval, while the magnetopause position is closer to the empirical estimate, within 10 % for 98.5 % of the 10-day period.The relative differences averaged over the 10-day period are 10.8-6.7 % (depending on used method) for the bow shock position and 3.4 % for the magnetopause position.Using our methodology and the GUMICS-4 simulation position, we find the shock and the magnetopause consistently closer the Earth than the empirical references, consistent with earlier findings (Janhunen et al., 2012).As relative differences of the order of 10 % suggest absolute differences of the order of 1 R E , the differences almost fall within the errors introduced by the 0.5 R E spatial resolution in the simulation run.The differences between empirical and simulation estimates did not increase even when abrupt changes in the solar wind parameters take place, indicating that GUMICS-4 is correctly capturing the large-scale dynamics.
Next, we took a 9 h subset of the 10-day data interval, combined it with five different 3 h synthetic initialisation phases, simulated the resulting five 3 + 9 h intervals, and compared the results to the original 10-day run in order to identify the effect of using different initialisation methods.Magne-tospheric response to synthetic data in global MHD simulation has been studied previously using GUMICS-4 by Palmroth et al. (2001).Our results demonstrate that, regardless of the initialisation parameters, the GUMICS-4 simulations converge toward a common evolution within 1 h.
However, comparing difference images produced on the noon-midnight meridian plane (XZ plane) using ρV x variable suggests that some differences in the magnetospheric state remain for extended periods following the initialisation phase.These differences depend on the initialisation parameters used and can affect the results several hours afterwards, as is shown in Table 3.The table shows that 5.1-7.2 % of the grid points show differences above 5 % even 9 h after the initialisation stage has ended.On the other hand, it can also be seen that the relative differences in ρV x greater than 30 % are statistically insignificant on timescales of several hours.This implies that the small localised deep red regions visible in Figs.7 and A2-A4 have only minor effects on differences between the simulations.
While the results do not clearly point to any of the initialisations as being superior to the others, there are couple of indications that may be useful for future studies.First, the dayside magnetosheath does not show any differences between the different initialisations, indicating that the magnetosheath has no internal processes that would carry the memory of the past driving conditions.Second, introducing a sinusoidally varying IMF B Y during the initialisation caused strong twisting of the tail, whose effects remained in the magnetotail for several hours.Thus, we would recommend using constant and/or near-zero IMF B Y preferably close to the observed value in the simulation initialisation.Third, the initialisations starting with strongly reconnecting magnetosphere driven by negative IMF B Z seem to produce more reliable results in the magnetotail.As shown for the B Y case, the magnetotail has a long memory that can retain effects of past configurations for quite some time.Strong reconnection has the effect of removing these by enhancing the convection cycle.Thus, we would recommend starting the initialisation with negative IMF B Z .
The results suggest that choosing the switching moment (14 June 15:00 UT) from artificial initialisation to nonartificial input data differently would have no impact on our recommendation by creating qualitative differences between the runs.Such change would naturally alter the behaviour of IMF B Z at 15:00 UT for each simulation.However, the current study is already covering the switching of IMF B Z from +5 to −7 nT (st1, lin2) and from −5 to −7 nT (st2, lin1, sin) as depicted in Fig. 3.If the change in the switching time would cause qualitative differences the simulations would be divided in the aforementioned two groups by the behaviour of IMF B Z at 15:00 UT.Such a division is not apparent in Table 3, however.
We identified two types of errors that arise in the simulations.Large but localised errors are observed in regions of strong gradients, at the plasma sheet boundary, when the pa-rameters change significantly over a short distance.Such errors were shown to be statistically insignificant, but can cause major discrepancies when the simulation results are compared with in situ measurements along the spacecraft orbit.Thus, caution is necessary when interpreting such comparisons.
Medium-size relative errors in the range 5-30 % occur over larger portions of space and imply dynamic events that do not occur simultaneously in the two simulations.Such errors vary in location and do not have a continuous temporal evolution, which again points out that, despite such differences, the simulations on the large scale converge toward a common solution.
In order to examine the cause of the differences, we plot in Fig. 9 time series of the relative errors in the three error ranges for the time interval 16:00-20:00 UT on 14 June 2007.The top panel shows the IMF clock angle for reference.The second panel shows the energy flux incident at the magnetopause evaluated from (Palmroth et al., 2003) where u is the total energy density, p pressure, B magnetic field, V flow velocity and E ×B the Poynting flux.The three bottom panels show the errors in the ranges 0-5, 5-30, and > 30 % for each of the initialisation runs.The grey bars note periods when the accuracy for at least one of the initialisation runs falls below the typical 93-95 % (panel c), indicating periods when errors are observed between the 10-day run and the different initialisation runs.Dynamics in the magnetosphere is controlled by the energy entering through the magnetopause in the form of Poynting flux (Palmroth et al., 2003).The total energy flux of the solar wind expressed by Eq. ( 4) is a strong function of B Z and thus moderate difference in B Z value results in moderate difference in transferred energy amount from the solar wind to the magnetosphere.The first two difference peaks are associated with a peak in the energy flux, featuring an about 70-80 % increase in the incident energy flux.Interestingly enough, the associated simulations lin2 and st1 are the only ones with northward IMF B Z at the end of the initialisation period.It is possible that the two first peaks are associated with the rapidly changing energy entry to the magnetosphere with a profile different from the 10-day simulation.
The last two maxima are preceded by a rapid change in sign of IMF B Y from positive to negative and back to positive.The actual peaks coincide with rapid rotations between northward and southward orientations.However, the start of the increase in the errors is not marked by distinct changes in the IMF or the energy flux variations.
The present study has utilised the relative differences in ρV X to quantify the global differences in the simulation runs at the maximum grid resolution of 0.5 R E .This configuration has been used extensively in previous studies (Palmroth et al., 2001) and therefore any differences we observe are likely induced by the initialisation as opposed to the chosen model setup.For completeness, we should mention that an alternative approach would be to examine global variables such as the cross-polar cap potential in the ionosphere.However, this would require increasing the adaptation level to 5 (0.25 max grid resolution) rather than 4, which we have used here.Even so, even with a finer grid resolution, GUMICS-4 produces relatively low polar-cap values (Gordeev et al., 2015) due to all the currents closing through the polar cap; this is a result from too weak region 2 currents.Therefore, the validity of this approach with regards to the initialisation phase is unclear, and more investigation is required.For that reason, we have focused exclusively on the magnetospheric region which has been validated extensively, and reported in the existing literature.Having said that, the investigation of the polar-cap potential will be the focus of future studies, but this is beyond the scope of the present paper.
It should be noted that a global MHD simulation does not necessarily converge towards the same solution even if resolution is improved (Ridley et al., 2010).In fact, if resolution is improved, numerical diffusion decreases, which means that the role of MHD physics as an error source grows.Our results imply that the differences in the simulation runs arise from the different initial conditions rather than numer-ical accuracy, as the solutions converge towards the same large-scale state in the long run.Thus, the results might be generalisable to all the MHD codes applying similar physics (ideal MHD) to that of GUMICS-4.
Analysis of the temporal behaviour of differences appearing in both large and small spatial scales indicates that differences in the detailed locations of current and flow systems in the (inner part) of the magnetotail maintain a memory of about 3 h.Thus, while for most purposes the different initialisations lead to rapid convergence toward nearly the same solution, there are local-and large-scale differences that can be seen hours after the initialisation period.Future work should extend the investigation of the large-scale dynamics causing the observed differences in the ρV x variable, including computation of the energy budget of the magnetosphere.

Conclusions
The main results can be summarised as follows: 1. Based on rather steady time evolution of the differences between empirical and simulation estimates to the bow shock and the magnetopause locations, GUMICS-4 is correctly capturing the large-scale dynamics during the simulated 10-day period.
2. Initialising the simulation with different upstream conditions affects the magnetospheric dynamics for approximately 3 h after initialisation in GUMICS-4.The differences in the magnetospheric states at the end of each initialisation phase are seen in both large-scale dynamics of the tail and small-scale localised differences of the flows and currents for several hours after the initialisation.
3. While the differences are small most of the time, they may grow as a response to, for example, fast IMF clock angle rotations after initialisation period, creating differences in momentum transport in the magnetospheric tail region.
4. A general conclusion can be given that, to quite a high degree of accuracy, the GUMICS-4 simulation converges within 1 h when considering the outer parts of the magnetosphere.Optimal initialisations use small and constant IMF B Y , start with negative IMF B Z , and leave preferably longer than 1 h period with actual observations before the simulation results are used for comparisons with observations.

Figure 1 .
Figure 1.Solar wind and IMF conditions during a 10-day period from 12 to 22 June 2007.Panels from top to bottom: (a) IMF B Y and (b) B Z in nT, solar wind speed (c) V in km s −1 , and density (d) n in cm −3 .Vertical red lines indicate 9 h data interval which was combined with initialisation periods.Non-labelled panels at the top and bottom indicate data gaps longer than or equal to 10 min.

Figure 2 .
Figure 2. Bow shock and magnetopause positions: (a) bow shock standoff distance in R E from the OMNI data (black) and GUMICS-4 (coloured lines), (b) relative difference between the shock position obtained using GUMICS-4 and OMNI, (c) magnetopause standoff distance in R E from the Shue model (black) and GUMICS-4 (magenta), and (d) relative difference between the GUMICS-4 magnetopause position and the Shue model.The relative differences are computed from 100•(reference−GUMICS)/reference.Vertical red lines indicate the 9 h data interval which was combined with initialisation periods.

Figure 3 .
Figure 3. Five different initialisations for the simulation.Panels from top to bottom: IMF (a) B Y and (b) B Z in nT, and solar wind velocity(c) V X in km s −1 and (d) density in cm −3 at x = 32 R E .The initialisations are shown with different colours (st1 in blue, st2 in blue with symbols, lin1 in magenta, lin2 in magenta with symbols and sin in green, 10-day (10d) run in black).The 10-day run data are from OMNI.The yellow background highlights the 1 h interval during which the magnetosphere should be formed(Janhunen, 1996;Raeder et al., 2009).

Figure 4 .
Figure 4. Magnetospheric configuration and the Earth's dipole field in 14 June 2007 at 16:30 (a) and 19:40 (b).Both the Earth's magnetic field and the IMF field as well as ρV x flow lines are showing.White dots show the time series points in the XZ plane.Labels (a-i) refer to Figs. 5, 6, and A1.Note that the figures are cut planes of 3-D plots.

Figure 5 .
Figure 5. Magnetotail plasma parameters in the lobe (a, c) and plasma sheet (b, d) at X = −8 R E Y = 0 R E Z = Z sheet (a, b) and X = −15 R E (c, d).Z sheet refers to varying Z coordinate in the plasma sheet.In each panel, subpanels from top to bottom: plasma flow speed in km s −1 , absolute difference in plasma velocity V X in km s −1 , density in cm −3 and B Y in nT, all in GSE coordinates.Absolute difference is calculated as 10d − x, where x = st1, st2, lin1, lin2, sin.

Figure 6
Figure 6 shows the bow shock and the magnetopause position (panel a), the V (panel b), B (panel c), plasma beta (β = 2µ 0 p/B 2 , panel d), andE Y (panel e) at X = 9 R E , Y = 6 R E , Z = 0 R E .Note that since the magnetopause and the bow shock are in motion as depicted in panel (a), this point is not necessarily always located within the magnetosheath but fluctuates relative to the magnetopause position crossing to the upstream solar wind when the magnetopause moves inward.This can be seen in the time series as small values of the magnetic field magnitude in panel (c) that often coincide with the bow shock and the magnetopause moving closer to the Earth in panel (a).The time series measured at X, Y, Z = [9, 6, 0] show only minor differences between the runs following the first hour of similar solar wind input, in-

Figure 6 .
Figure 6.Panels from top to bottom: (a) bow shock and magnetopause nose positions in R E , (b) plasma speed V in km s −1 , (c) magnetic field magnitude B in nT, (d) plasma β = 2µ 0 p/B 2 and (e) electric field E Y in mV m −1 .Parameter time series in panels (b-e) are produced at X= 9 R E Y = 6 R E Z = 0 R E (inGSE coordinates).Yellow background highlights 1 h interval during which the magnetosphere should be formed(Janhunen, 1996;Raeder et al., 2009).

Figure 7 .
Figure 7. (a) Scaled (by solar wind value) ρV x in the 10-day run at 18:30 on 14 June and (b-f) the magnitude of the relative difference (defined as |(ρV i x − ρV 10d x )/ρV 10d x | • 100) in ρV x in the noon-midnight meridian plane.Note that the scaling of (a) is not shown.

Figure
Figure A3.(a) Scaled (by solar wind value) ρV x in the 10-day run at 17:12 on 14 June and (b-f) the magnitude of the relative difference (defined as |(ρV i x − ρV 10d x )/ρV 10d x | • 100) in ρV x in the noon-midnight meridian plane.Note that the scaling of (a) is not shown.

Figure
Figure A4.(a) Scaled (by solar wind value) ρV x in the 10-day run at 19:13 on 14 June and (b-f) the magnitude of the relative difference (defined as |(ρV i x − ρV 10d x )/ρV 10d x | • 100) in ρV x in the noon-midnight meridian plane.Note that the scaling of (a) is not shown.