Global MHD Simulations of the Earth’s Bow Shock Shape and Motion Under Variable Solar Wind Conditions

Empirical models of the Earth’s bow shock are often used to place in situ measurements in context and to understand the global behavior of the foreshock/bow shock system. They are derived statistically from spacecraft bow shock crossings and typically treat the shock surface as a conic section parameterized according to a uniform solar wind ram pressure, although more complex models exist. Here a global magnetohydrodynamic simulation is used to analyze the variability of the Earth’s bow shock under real solar wind conditions. The shape and location of the bow shock is found as a function of time, and this is used to calculate the shock velocity over the shock surface. The results are compared to existing empirical models. Good agreement is found in the variability of the subsolar shock location. However, empirical models fail to reproduce the two-dimensional shape of the shock in the simulation. This is because significant solar wind variability occurs on timescales less than the transit time of a single solar wind phase front over the curved shock surface. Empirical models must therefore be used with care when interpreting spacecraft data, especially when observations are made far from the Sun-Earth line. Further analysis reveals a bias to higher shock speeds when measured by virtual spacecraft. This is attributed to the fact that the spacecraft only observes the shock when it is in motion. This must be accounted for when studying bow shock motion and variability with spacecraft data.


Introduction
The interaction of the solar wind, streaming out radially from the Sun, with the magnetosphere leads to the formation of a bow shock upstream of the magnetopause, with a turbulent sheath region (the magnetosheath) in between (Burgess & Scholer, 2013;Fairfield, 1976;Krasnoselskikh et al., 2013;Lucek et al., 2005;Ness et al., 1964). While the main role of the shock is to decelerate and deflect the supersonic solar wind around the magnetosphere, reflection of particles at the shock leads to the formation of the foreshock, which leads to numerous wave-particle interactions and particle acceleration (Eastwood et al., 2005). In situ data provide a wealth of information about the local plasma properties and the various physical processes controlling thermalization, particle acceleration, etc. However, knowledge of the context of these in situ measurements is required to fully understand the observed behavior. Variability in the solar wind causes dynamic movement of the outer magnetospheric boundaries as they respond to changes in pressure and interplanetary magnetic field (IMF). Consequently, empirical models of the shock shape and location, derived statistically from spacecraft observations, are crucial in this regard. They are fundamental to many spacecraft investigations of shock and foreshock physics, as they are used to place local measurements in the global context, and to correctly interpret shock geometry. However, while such models agree well with the average location and shape of bow shock or magnetopause, they often need to be scaled in order to agree with the specific spacecraft crossing location for any particular event (Schwartz, 1998).
Early work on modeling the Earth's bow shock included quantifying the position and shape of the bow shock under different solar wind and IMF parameters, and the position and shape of the magnetopause. Following on from the Spreiter et al. (1966) "hydromagnetic" model, various empirical models were derived using a combination of spacecraft data for statistical fitting and gas dynamic theory for the basic shock shape (Farris & Russell, 1994;Filbert & Kellogg, 1979;Formisano et al., 1979;Peredo et al., 1995). Empirical models, constructed from databases of spacecraft shock crossings, usually assume several first-order MEJNERTSEN ET AL.
MHD SIMULATIONS OF EARTH'S BOW SHOCK 1 approximations. Most models employ axisymmetric conic sections for the bow shock, although more complex models have been developed (Schwartz, 1998). The magnetopause is also assumed to be a similar conic section: however, Lin et al. (2010) derives a full three-dimensional magnetopause shape, including cusps, from spacecraft observations.
These models have shown that the bow shock shape and position is largely determined by the solar wind dynamic pressure, and to a lesser extent, the solar wind Alfvén Mach number (Merka et al., 2005;Peredo et al., 1995). In using empirical models, it is effectively assumed that the solar wind pressure is uniform over the bow shock surface. It is also assumed that there is an instantaneous relationship between changes in the magnetopause and bow shock location, and they do not account for the relevant signal speeds in the magnetosheath, nor the motion of the boundaries. There are a number of empirical models in the literature (e.g., Cairns & Lyon, 1995;Jeřáb et al., 2005;Merka et al., 2005), all derived using different data sets and reference surfaces. As a result, these models differ from one another when applied under the same solar wind conditions. Comparisons of these models have been performed in numerous papers. For example, Turc et al. (2013) compared a number of models when the solar wind Alfvén Mach number was low, cross-referencing them with Cluster crossings of the bow shock. They found that in low Alfvén Mach numbers, the Jeřáb et al. (2005) model performed best.
The launch of Cluster led to routine calculation of shock normal speeds and orientations, which can be obtained using four spacecraft timing analysis (Schwartz, 1998). Subsequent work using the Cluster fourspacecraft timing analysis has shown that model prediction of position and shape typically perform reasonably well when compared with observations (Horbury et al., 2002;Meziane et al., 2014). Horbury et al. (2002) analyzed 48 crossings of the shock where the shock normal and velocity could be obtained. They found that typical shock velocities were~35 km/s and that the shock normal agrees with empirical models. Comparisons with the coplanarity method of finding shock normal, which uses single spacecraft, performed poorly. Meziane et al. (2014) extended Horbury et al.'s (2002 analysis by including more shock crossings, agreeing with Horbury et al.'s (2002) findings. It was found that slow shock speeds (<80 km/s) follow a Gaussian distribution of width 42 km/s. As found in the derivation of the empirical models, the spacecraft observations show that the bow shock velocity is largely driven by the changes in the solar wind dynamic pressure and has little dependence on Mach number (Meziane et al., 2014). Furthermore, the speed of shock, as calculated from the time derivative of empirical model shock position, agrees statistically with observation as well (Meziane et al., 2015).
Moving beyond typical solar wind variability, discrete structure (specifically, interplanetary shocks) causes rapid and significant boundary motion (Dryer et al., 1967;Grib et al., 1979;Ivanov, 1964;Shen & Dryer, 1972). Early work included Greenstadt et al. (1972), who used the Ogo 5, Heos 1, and Explorer 33 to calculate the shock velocities. Völk and Auer (1974) applied a theoretical analysis of shock motion by considering the shock under different solar wind Alfvén Mach numbers. The observed high shock velocities under low Mach numbers originate due to an increased sensitivity to Alfvén waves, whereas high Mach numbers and associated high shock velocities are due to the variability introduced by traveling interplanetary shocks (Völk & Auer, 1974). A number of studies have simulated the effect of an observed interplanetary shocks on the Earth magnetosphere, finding that the shock-bow shock interaction causes a new train of discontinuities (Pallocchia, 2013;Pallocchia et al., 2010;Přech et al., 2008;Šafránková et al., 2007). As they propagate through the shock and sheath, these have been shown to create subsequent shocks, which deflect and reflect off the magnetopause. Simulations by Samsonov et al. (2007) show Earthward bow shock and magnetopause motions as the interplanetary shock impacts, but outward motions due to a reflected shock. Similar results have been observed by spacecraft (Pallocchia, 2013). Such shocks have a profound effect on the magnetosphere, creating disturbances which propagate through the magnetosphere at speeds higher than the original disturbance in the solar wind (Andréeová et al., 2008;Wang, 2005;Xiao-Cheng et al., 2005), and create significant responses in the inner magnetosphere (Andreeova et al., 2011).
Empirical models are therefore limited in the extent to which they can capture shock motion that occurs in response to rapid variation in the upstream solar wind. In particular, it is not well understood as to when these statistical models become no longer appropriate, nor is it clear as to how these models differ from reality when taking into account the dynamic transition of boundaries to new equilibrium states. Global Journal of Geophysical Research: Space Physics 10.1002/2017JA024690 magnetohydrodynamic simulations are a useful tool in this regard, as they can simulate the whole systemself-consistently.
In this work, global MHD simulations of the Earth's magnetosphere under variable solar wind are performed using the Gorgon code (Chittenden et al., 2004;Ciardi et al., 2007;Mejnertsen et al., 2016). This work examines the motion of the bow shock. The documented case study from Cluster covering variable solar wind from Maksimovic et al. (2003) is considered. The layout of this paper is as follows. The Gorgon simulation is briefly introduced, along with the solar wind profile used and the methods for finding the boundary locations. The case of the motion of the bow shock in the subsolar region is then examined, followed by a qualitative analysis of the 2-D bow shock shape. Finally, the velocity of the shock is statistically examined and compared with virtual spacecraft observations.

The Gorgon MHD Code
Gorgon is a 3-D magnetohydrodynamic code, initially developed for studying high energy, collisional plasma interactions such as Z-pinches (Chittenden et al., 2004;Jennings, 2006;Jennings et al., 2010), laser-plasma interactions (Smith et al., 2007), and magnetic tower jets (Ciardi et al., 2007). It has recently been adapted to simulate planetary magnetospheres and their interaction with the solar wind (Mejnertsen et al., 2016).
Gorgon uses a fully explicit, Eulerian formulation of the resistive semiconservative MHD equations for a fully ionized hydrogen plasma, as given by equations (1) to (6): The equations describe conservation of mass (1), momentum (2), proton (3) and electron internal (fluid) energy (4), the magnetic induction (equation (5)), and Ohms law (6) where ρ is the mass density, v ! is the bulk plasma flow, P p, e is the proton/electron pressure, J ! is the current density, B ! is the magnetic field, A ! is the vector potential, ε p, e is the proton/electron internal energy density, η is the plasma resistivity, c is the speed of light, and ε 0 is the permeability of free space.
Due to its roots in high-energy plasma physics, the MHD formulation in Gorgon is atypical: it treats the electron and proton internal energy equations separately (equations (3) and (4)), allowing them to be out of thermodynamic equilibrium. These equations include terms for Ohmic heating η J ! 2 , optically thin radiation losses Λ, and electron-proton energy exchange Δ pe . In the parameter regime of magnetospheric space plasmas, these terms are negligible, and hence have been disabled in the code or are also negligibly small. The pressure is calculated using the ideal gas law, with γ = 5/3. The choice of solving for the internal energy, rather than the total energy, means that the total energy is not strictly conserved: however, this mitigates effects such as negative pressures occurring in regions of high magnetic field strength.
Gorgon makes use of second-order Van Leer advection to solve the advection terms in equations (1)-(6). It also employs a variable time step, which automatically satisfies the relevant Courant conditions. A Von Neumann artificial viscosity, as well as the deBar correction, is implemented in order to accurately capture the shock and improve energy conservation (Benson, 1992).

10.1002/2017JA024690
The magnetic field solver (equations (5) and (6)) uses the vector potential representation of the magnetic field , spatially separated from the fluid equations on a staggered grid (Yee, 1966). This allows the magnetic divergence condition to be satisfied automatically to machine accuracy, without using divergence cleaning algorithms. Equation (5) is derived directly from the Ampere's law and contains a plasma ( J ! /ε 0 ) and The first describes the electromagnetic effects in plasma, using a standard Ohm's law (equation (6)) to account for resistive diffusion and advection. The second term describes electromagnetic wave propagation through a vacuum. This allows the electromagnetic fields to propagate through regions of very low mass density, defined by a set mass density threshold in the code, set at ρ = 10 À23 kg m À3 , below which plasma properties and forces are set to zero. Since the electromagnetic waves propagate at the speed of light, this part is sub cycled over the fluid time step, with a reduced speed of light.
Many global magnetosphere codes employ a split magnetic field setup (Hu et al., 2007;Lyon et al., 2004;Ogino et al., 1992;Powell et al., 1999;Tanaka, 1994), where the magnetic field is split into two components containing the static dipole field and the deviation from the dipole field: this removes discretization errors near the planet. In this paper, we do not use a split magnetic field as internal investigations have revealed that the split magnetic field version of Gorgon produce the same results as the version in this paper regarding the position of the magnetopause and bow shock.

Simulation Setup and Solar Wind Input Conditions
In this work, the region of interest contains the outer magnetospheric boundaries. Hence, a simulation area spanning x GSE = (À50, 30) R E , y GSE = (À20, 20) R E , and z GSE = (À30, 30) R E is used, with a uniform grid size of 0.2 R E . The simulation coordinate system is consistent with geocentric solar ecliptic (GSE): the dipole is located at the origin, and to avoid complications associated with dipole tilt or rotation in this analysis, it is aligned to the z GSE direction. The inner boundary is set at 3 R E and is centered on the origin. It consists of a dipole source at the center of the Earth, as well as an artificial cutoff, which numerically removes plasma crossing the inner boundary. The vacuum solution of the code is solved for the magnetic field in this area, which allows the magnetic field to vary on the inner boundary: this provides a resistive inner boundary. The simulation is initialized as empty (vacuum), which allows the terrestrial dipole to expand into the box.
The outer boundaries are set to be free-flow (Dirichlet conditions), except for the sunward (+x GSE ) edge, where solar wind and IMF parameters as observed by ACE were applied uniformly across this face of the simulation domain. One-minute resolution data observed on 31 March 2001 were used, corresponding to the Cluster spacecraft case-study performed by Maksimovic et al. (2003). The data were first ballistically propagated from the location of ACE to the sunward edge of the simulation box (x = 30 R E ). Since the simulation time step is smaller than 1 min, the parameters were linearly interpolated between observations. The interval of interest is 17:14:00 UTC-19:34:00 UTC. The simulation was initialized with the preceding 2 h of ACE observations. Consequently, a simulation time of t = 0 min corresponds to 15:14:00 UTC, and the analysis interval begins at t = 120 min.
The simulation input data corresponding to the analysis interval is shown in Figure 1. The proton number density ( Figure 1a) ranges between 10 cm À3 and 30 cm À3 . The average solar wind speed is approximately 550 km/s (Figure 1b). The strength of the IMF (Figure 1c) was larger than typical IMF strength, ranging between 20 nT and 35 nT. The IMF orientation is strongly southward, with B z ranging between À32 nT and À20 nT. Since the solar wind velocity is predominantly in the Àx GSE direction, B x cannot be propagated into the box. A constant B x is initialised throughout the simulation domain, whose magnitude corresponds to the average B x measured by ACE. We note that with these solar wind and IMF parameters, the solar wind Alfvén Mach number is low, averaging M A ≈ 5. Empirical models typically do not perform as well under low Mach numbers as there are fewer crossings in this regime: the best model for low Alfvén Mach numbers has been established to be the one proposed by Jeřáb et al. (2005) (Turc et al., 2013) and is used here to compare with the simulation output.
The state of the magnetosphere at the beginning of the analysis period (t = 120 min) is shown in Figure 2. It shows the number density and the dynamic pressure in slices along the ecliptic and noon-midnight planes. At this time step, the simulation has formed a full magnetosphere with its characteristic features: a low density solar wind, arriving from the left; the bow shock and the turbulent magnetosheath; and a low density

Locating the Bow Shock
To characterize the shape and location of the shock, we examine these properties in the x GSE y GSE (ecliptic) plane and the x GSE z GSE (noon-midnight) plane. In this section, the shock finding algorithm is explained by referring to the shock in the ecliptic plane. The same method is used for the polar plane, except all instances of y GSE are substituted by z GSE .
First the bow shock is located along the x axis, by seeking the compression of plasma due to the shock. Technically, this is done by calculating the mass density gradient along the x axis, and identifying negative peaks, as illustrated in Figure 3a. In the simulation, there are often many such peaks, only one of which corresponds to the bow shock. To determine which is the shock, "significant" peaks are found: ones whose scaled magnitude is larger than a set value. This method is found to work well as the bow shock is a strong shock and is usually one of the, if not the largest, such peak. If there are multiple significant peaks, the peak most sunward is chosen, as there tend to be no such peaks in the solar wind.
In global MHD simulations, the position of the shock is always located on the grid cell interface, since Eulerian simulations are a   (Figures 2b and 2d). The magnetosphere is largely denoted by the regions of low number density and dynamic pressure, shown in dark blue. The white dashed line denotes the inner boundary of the simulation, at 3 R E . In Figure 2d, there is an enhanced dynamic pressure along the magnetopause, which is due to reconnection jets on the dayside. In Figures 2a  and 2b, there is a kink in the magnetopause position at x = 4, y = À7: This is due to a flux rope propagating downstream along the magnetosphere.
Journal of Geophysical Research: Space Physics 10.1002/2017JA024690 set of calculations solving the MHD equations on finite sized cells. To overcome this, a subgrid method was used, as well as smoothing of the bow shock standoff position curves. Once the location of the grid position of the bow shock has been found, the bow shock peak position and the two adjacent grid points in the x direction are taken and a parabolic function was fitted dρ dx ¼ ax 2 þ bx þ c , passing through these three points. This is shown as the purple curve in Figure 3a. The subgrid shock position is given as the maximum of this parabola, x BS ¼ À a 2b , and corresponds to the black dashed line. This process was repeated for all y GSE grid values, and for each time-step. Hence, the 2-D shock position x BS (y GSE , t) is obtained as a function of y GSE and t. The result of this analysis is shown in Figure 2b (red curve).
The object of this analysis is to explore the motion of the bow shock. The corresponding shock speed along the x direction is shown in Figure 2c (red curve). Although the large-scale motion of the shock is visible, the subgrid parabolic fitting method does not completely eliminate the impact of the grid as the shock passes from one cell to the next. Therefore, the shock position is smoothed as a function of time using a Savitsk-Golay filter: this algorithm fits a polynomial to the data in a rolling window (blue curve in Figure 2b). Once smoothed, a second-order centered difference gradient is used to obtain the shock velocity, as shown in blue in Figure 2c. This combined filtered subgrid method removes grid fluctuations in the shock position and velocity, leaving the underlying behavior visible. This method is simpler and more effective than using a Rankine-Hugoniot approach since it assumes a steady state solution, which is not generally the case in this simulation, and mitigates issues, which arise in how the upstream and downstream definitions are determined.
The slight overshock present in the mass density curve (red) in Figure 3a is a numerical effect of the method used in solving the MHD equations on a grid. Many MHD codes use an artificial viscosity to prevent such overshoot, though if the viscosity is too large, the shock is smeared over more grid cells, leading to a wide shock. The parabolic fitting method and temporal smoothing used here overcome this numerical effect and minimize the effect of the overshoot.
Finally, the subsolar magnetopause is located by identifying the magnetopause current sheet, and the corresponding peak in the magnitude of the current density. The subgrid magnetopause position and speed can then be located using the same technique as for the bow shock.
The empirical models are calculated as described in Jeřáb et al. (2005) for the bow shock and Shue et al. (1998) for the magnetopause, using the solar wind conditions as input into the simulation. The velocity of the model position is calculated using the same smoothing techniques as with the simulation positions.   The MHD model performs well in capturing the overall motion of the bow shock and magnetopause, compared to the scaled empirical models. For example, Figure 4a shows that large variations in the solar wind dynamic pressure occur starting at t = 130 min and t = 190 min. The sharp increase in dynamic pressure at t = 130 min causes the bow shock and magnetopause standoff positions to move earthward, as expected. The pressure then reduces slightly, causing the shock and magnetopause to move sunward. Both the empirical models and Gorgon show the boundaries responding on similar overall timescales. However, they do share common features, such as the alternating magnetopause velocity peak at t = 130 min, corresponding to an earthward movement, followed by a sunward movement.

Comparison Along the Sun-Earth Line
The drop in pressure at t = 190 min has more complex features. Leading up to the peak, from t = 170 min both the dynamic pressure and the Mach number oscillate as they rise. This causes the bow shock and magnetopause to rapidly reposition, moving back and forth along the x GSE direction. The bow shock velocities agree well, with similar peaks overlapping between Gorgon and model. During this interval, the empirical model bow shock velocities tend to be larger than the simulation velocities, which is a general trend throughout. However, the simulated magnetopause is much more variable than predicted by the Shue et al. (1998) empirical model. This is echoed in Figure 4c, which shows the ratio between the bow shock and magnetopause distances, and where there is more small-scale variability in the simulation result (red). Closer analysis of the simulated magnetopause shows that the magnetopause current sheet is often disrupted by the presence of reconnection jets and associated structure. However, this does not appear to affect the overall location of the shock, with the magnetosheath acting in some sense as a low-pass filter. We note that in more recent Vlasov simulations, magnetopause reconnection and flux rope formation can affect the bow shock locally, but this does not impact the goal of the present work to understand the shape and motion of the shock on global scales (Pfau-Kempf et al., 2016).

Overview of Two-Dimensional Boundary Dynamics
We now examine in more detail the interval t = 131 min to t = 135 min when the bow shock moved earthward and then sunward over a period of 4 min in response to a sharp increase and small decrease in solar wind dynamic pressure (denoted by the grey region in Figure 4). Figures 5a and 5b each show eight snapshots of the dynamic pressure in the ecliptic and noon-midnight planes, with time increasing from left to right. Figures 5c and 5d show the positions of the shock in the ecliptic and noon-midnight planes at those times, calculated according to the procedure described in section 2. The color of the line indicates the time. Figure 5e shows the corresponding position of the empirical model shock, which is cylindrically symmetric. Figure 4d immediately shows that the simulation bow shock shape has cusp like features and is not a simple conic section in the noon-midnight plane, as assumed by many empirical models. Furthermore, as the bow shock moves, the position of the cusps does not stay in the same place. When comparing the shock normal between the simulation shock and the model shock, the simulation shock normal deviates on average by 8°from the model shock normal, with the majority of the deviation being within 20°: this is similar to the statistical results obtained in (Meziane et al., 2014).
As mentioned in the previous section, during this interval the solar wind dynamic pressure increased and then reduced. Figures 4c and 4d show that the subsolar point is the first to move earthward: this is expected as the solar wind front reaches this position on the shock first. However, the flanks of the shock during the first three snapshots (131.0 to 132.2 min) barely move from their original position. At intermediate times, 132.8 min to 134.00 min, the flanks are now moving antisunward whereas the shock subsolar point is approximately stationary and appears to reach equilibrium. The time from when the subsolar point first moves and when the flank at 90°latitude moves is approximately 2 min, which corresponds to a speed along Àx GSE of approximately 530 km/s, comparable to the solar wind speed. In the final two snapshots the subsolar point moves sunward. This is due to the decrease in solar wind dynamic pressure near the nose as can be seen in the change in color from yellow to green in panels 5 to 8 of Figures 5a and 5b. However, this change has not yet reached the flanks, which continue to move antisunward. This behavior is not seen in the model bow shock position (Figure 4e), whose overall profile moves instantaneously and uniformly with the upstream solar wind conditions.

Journal of Geophysical Research: Space Physics
10.1002/2017JA024690

Shock Motion on the Dayside Flanks
To better understand the global motion of the shock, its variability and speed away from the Sun-Earth line are now examined. In particular, we focus on the radial position and velocity of the shock as a first step. Using the method described earlier, the bow shock position is found in the Cartesian GSE coordinate system. To calculate the radial motion, the GSE position is first converted to a two-dimensional polar coordinate system centered on Earth, where r is the radial distance from the center of the Earth and ϕ is the azimuth from the positive x GSE axis (defined such that tanϕ ¼ y GSE xGSE ). The location of the shock is then interpolated along lines extending radially from the origin. This gives the position of the shock as a function of ϕ. By doing this for each time-step and taking the gradient between each time step, the radial velocity is obtained. Figure 6 shows the radial bow shock position along three radial lines: at ϕ = 0°(the Sun-Earth line, shown in red), at ϕ = 40°(blue) and at ϕ = 80°(green).
The standoff positions (Figure 6a) along each flank direction show similar behavior to the Sun-Earth line. The flank shock positions have larger standoff positions, which is expected due to the overall shape of the shock (Figure 4). While the general pattern of the flank shock motion follows that of the subsolar shock (with some delay, as noted above), it is more exaggerated, creating larger bow shock velocities. Hence, it is found that the shock velocities on the flanks are, in general, much larger than at the Sun-Earth point. This can also be seen in Figure 5. Figure 6a also highlights the effects of the solar wind propagation delay on the flank shock location. This can best be seen in the velocity graph in Figure 6b, comparing the 0°line (red) and 80°line (green): the 80°l ine trails the 0°line. One example of this is at t ≈ 130 min, denoted as the grey shaded area in Figure 6: this region corresponds to the same period discussed in section 3.2. The subsolar shock initially moves earthward due to the increase in solar wind dynamic pressure. The flanks also begin to move earthward, a couple of minutes later than the subsolar shock. However, by this time a new phase front with lower solar wind dynamic pressure has arrived at the subsolar shock, which begins to move sunward. Consequently, the subsolar region is moving oppositely to the flanks toward the end of the shaded interval. Figure 5 therefore shows that the shock does not move coherently in and out at all locations in response to variations in solar wind pressure that occur on timescales less than the transit time of a phase front over the shock surface.
In principle, to improve the fidelity of an empirical model, one could account for the difference in the propagation time of the SW phase fronts between subsolar and flank, apply the local solar wind conditions to each location on the model shock, and thus trace out the distorted shock surface. However, we see here that the velocity of the shock complicates this, since the shock reacts differently at the flanks compared to the subsolar region.

Histograms of Bow Shock Speed
By taking the radial shock velocity for all ϕ and t, one can produce a histogram showing the distribution of radial shock speeds. Figure 7 shows the distribution of the bow shock velocity in the ecliptic plane. The time series data from the previous analysis from all radial lines is taken and binned in 10 km/s bins to create a histogram of the shock velocity distribution. A Gaussian profile is then fitted to the data. In general, Gaussian profiles fit the data well, especially for low velocities. At high velocities, there is a slight deviation with an enhancement in the tail of the distribution. Figure 6a shows that the radial shock velocity has a mean near zero, and a standard deviation of 21.4 km/s. The data were then divided into subsolar and flank regions (Figures 6b and 6c). There is a difference between velocities measured in the subsolar region ( Figure 7b) to those in the flank regions ( Figure 7c). The flank regions have a larger standard deviation, which reflects the fact that the shock tends to move faster on the flanks than it does at the subsolar point. These distributions are similar to those in Meziane et al. (2014), who derived the bow shock velocity distribution using Cluster observations. We also note that Meziane et al.'s (2014) analysis of shock velocities observed by Cluster The red is the Sun-Earth line, the blue is 40°duskward, and the green is 80°d uskward. The grey shaded region shows an example where the flank shock velocity peak is much larger in magnitude than the subsolar shock, and is delayed. indicated that the shock motion is radial. The mean of the shock speed in all three panels of Figure 7 is not exactly zero, but approximately 1-2 km/s radially outward. This is likely a consequence of the non-steady solar wind input used, where the shock does not necessarily end where it starts in the analysis period, implying a mean displacement further away radially from the Earth. However, this result is not significant, as this mean is comparable to the Gaussian error on it, 1.5 km/s.

Virtual Spacecraft
Part of the variability between different empirical models can be ascribed to the fact that they are based on different spacecraft data sets acquired in different regions of space (Schwartz, 1998). To better understand this, here we examine the simulation data using virtual spacecraft, and in particular focus on the distribution of shock speeds. This is done with the following method.
A prototypical satellite orbit is first defined. Here we use the apogee and perigee of Cluster's initial orbit (Escoubet et al., 2001). The orbit is confined to the ecliptic plane, and its apse line orientation is chosen at random between ϕ = [À110°, 110°]. The starting position of the virtual spacecraft at t = 120 min of simulation time is set randomly. In order to find shock crossings, the spacecraft and bow shock positions as a function of time are compared, locating points of intersection. The shock velocity as observed by the virtual spacecraft is then found by sampling the shock velocity time series at that position and time. This process was repeated for 201 virtual spacecraft and resulted in a similar number of crossings to that of Meziane et al. (2014). Unlike the four-spacecraft analysis, a timing analysis is not performed in this work, since we have the complete spatial information of the shock and thus fewer assumptions need to be made. Figure 8a shows all the virtual spacecraft orbits as black dashed lines. Shock crossings are shown as red dots, with the corresponding shock surface also shown.
At the time of each shock crossing, the radial speed of the shock is also known from the analysis described in the previous sections. The histogram of radial shock speeds associated with the virtual spacecraft crossings is shown in Figure 8b. The green curve is a Gaussian fit to the data in blue, and the black curve shows a similar fit for the velocity over the whole shock (the same Gaussian fit as in Figure 7a). The virtual spacecraft radial shock speed distribution has a standard deviation of 43 km/s, much larger than that of the full simulation distribution. If the number of virtual spacecraft is increased (hence increasing the number of crossings), the distribution in Figure 8 becomes smoother but does not tend toward the global distribution (Figure 7a). This suggests a fundamental bias to larger shock velocities in the virtual spacecraft analysis.

Discussion and Conclusions
A global MHD simulation has been performed on real solar wind parameters during an event analyzed by Maksimovic et al. (2003). The variability in the solar wind dynamic pressure causes the magnetosphere's outer boundaries to vary in response. Here the shape, size, and velocity of the bow shock is analyzed as it readjusts to the arrival of new solar wind phase fronts. The bow shock position is found to be Figure 8. The bow shock velocity distribution as measured by virtual spacecraft: (a) shock crossings calculated using the virtual spacecraft algorithm. The black dashed lines show the randomly generated elliptical orbits, and the red dots show where the virtual spacecraft crossed the shock. The corresponding shock surface is also shown. (b) The virtual spacecraft observed shock velocity distribution, together with a Gaussian fit (shown in green) and the Gaussian fit over the whole shock as in Figure 7a. The distribution observed by the virtual spacecraft is much wider than that observed in the simulation as a whole.
Journal of Geophysical Research: Space Physics 10.1002/2017JA024690 highly variable during this time: it expands and contracts largely due to changes in the solar wind dynamic pressure.
We note that the empirical models here are scaled to that of the simulation bow shock and magnetopause: this is done in order to compare the variability of the outer boundary motions, rather than the size and location. We note that empirical models almost always require scaling when used with satellite shock observations; furthermore, the discrepancy in the size of the outer boundaries between empirical models and global MHD simulations is common and could be attributed to a number of factors (Gordeev et al., 2015). However, it is generally thought that the relative dynamics are still applicable; more realistic inner magnetosphere models do not much improve the motion of the outer boundaries (Němeček et al., 2011).
In analyzing the motion of the subsolar region, it is found that the variability of the bow shock is comparable to that of the Jeřáb et al. (2005) empirical model. There is more variability in the magnetopause location in the simulation, but this does not appear to affect the location and shape of the bow shock on a global level, suggesting a low-pass filter effect.
Analysis of the shock surface as a whole shows that for typical variability on time scales of minutes, the nose and flank can be moving differently at the same time. This is essentially due to the solar wind transit time over the shock in the antisunward direction, with the subsolar region affected by the varying solar wind several minutes before the flanks. This can lead to situations where the subsolar region of the shock moves in a different direction to the flanks, due to the solar wind front arriving at the flanks a finite time later. Empirical shock models often assume a characteristic shape scaled by a single value of the solar wind ram pressure. The results presented here indicate that such models must be used with care when the solar wind dynamic pressure is varying on minute timescales and when spacecraft observations are far downstream of the subsolar point. In particular, there has been a long-standing discussion of whether the shock moves coherently in the x direction or self-similarly (i.e., radially). However, from Figure 5, this picture is too simplified. Since the shock is essentially a series of local phenomena, affected by local plasma parameters and waves, it can cause different parts of the shock to behave differently, due to the timing between solar wind fronts.
Regarding the velocity of the bow shock, the location of the flank shock is much more variable than the subsolar region, exhibiting a wider distribution in the radial velocity. Measurements of shock speed must therefore account for this in determining the dynamic properties of the shock itself. This was examined in more detail by flying a virtual spacecraft with a Cluster-like orbit through the simulation. A much wider distribution in shock speed is found, implying a bias in shock velocities observed by spacecraft. Spacecraft observations will observe a wider distribution of shock velocities due to two factors. The majority of orbits take the spacecraft to the flanks of the shock (Merka et al., 2005), which has a wider distribution than at the subsolar point (Figure 7). Virtual spacecraft move relatively slowly compared to the shock: this causes observations of the shock to be due to the shock moving over the spacecraft. When the shock moves quickly, it moves between a large range of positions and is observed more often. Hence, the spacecraft will observe slow-moving shocks fewer times than fast-moving shocks, leading to a wider distribution (seen in Figure 8).
In conclusion, while empirical models can be used to derive information about the shape and motion of the bow shock, they must be used with care when used in spacecraft data analysis. Variability on the timescale of minutes is sufficiently fast to invalidate assumptions of uniform ram pressure, which apply to many empirical models. More work is required to determine how (and if) empirical models can be modified to overcome these limitations. For example, one could dynamically scale the empirical model to the local upstream solar wind conditions for position on the shock, which could account for the solar wind propagation time. This could then be used to piece together a more accurate representation of the shock surface. However, this will not fully account for the motion, as the bow shock moves differently at the flank than at the subsolar region. Finally, spacecraft observations are likely to be biased in estimates of shock speed. The subsolar and flank shock can move in different directions simultaneously, and their speed distribution is also different. Although observations indicate that the shock is moving fast (tens of km/s) when it is observed, spacecraft can only observe the shock when it is in motion and when it passes over the spacecraft. Comparing the shock motion with simulations will provide more complete details. It is therefore necessary to revisit experimental Journal of Geophysical Research: Space Physics 10.1002/2017JA024690 observations of shock motion and their statistics to better understand the true underlying nature of the shock motion in response to the variable solar wind.