Over 20-year global magnetohydrodynamic simulation of Earth’s magnetosphere

We use the fifth major version of Grand Unified Magnetosphere-Ionosphere Coupling Simulation (GUMICS-5), which is a parallelized version of GUMICS-4 described extensively in [16].

In the magnetosphere GUMICS-4 and 5 solve the ideal MHD equations using a first-order conservative finite volume method. The MHD equations are solved primarily using Roe’s approximate Riemann solver [29] and in rare cases when that solution fails a more robust and numerically diffusive solver is used for that particular calculation. Similarly to GUMICS-4 and [31], the magnetic field is split into a static background dipole ($B_{0}$) and a perturbed part ($B_{1}$) which is modified by the MHD solver(s). The ionosphere is modeled as a spherical surface with radially integrated current continuity, from which the electric potential is solved using a dipole-field-aligned electric current and a radially integrated conductivity tensor. Ionospheric conductivity is controlled by the magnetospheric temperature and mass density and by a parameter characterizing the loss cone filling rate. The ionosphere - thermosphere interaction processes impacting conductivities are based on the MSIS model [13].

On a high level the ionospheric and magnetospheric grids of GUMICS-5 are identical to GUMICS-4: The magnetosphere is discretized as a cell-by-cell run-time adaptive octogrid and in GUMICS-5 it is built on a custom version of DCCRG [15] with the solution calculated in parallel using the Message Passing Interface. The ionosphere is modeled as a spherical surface and is discretized as a static cell-by-cell adaptive triangular grid. In GUMICS-4 the ionospheric grid had a maximum resolution of approximately 100 km only in the auroral oval region, while in GUMICS-5 the maximum resolution has been extended to everywhere poleward of $\pm 58^{o}$ magnetic latitude.

The initial objective for developing GUMICS-5 was to obtain identical results compared to GUMICS-4 with significantly faster time to solution. During subsequent development several features of the model were improved, e.g. removal of divergence of perturbed magnetic field and mapping between ionosphere and magnetosphere along dipole field lines, and the code was updated to a more recent C++ standard used by current compilers. Perhaps most importantly, a rotating background dipole field was added that was not available in GUMICS-4. With these modifications, results between GUMICS-4 and -5 are very similar overall but not identical.

The interval 1998-2020 is divided into simulations of 1 UTC day with an additional 4 h overlap with the previous day’s simulation. A 4 h overlap was considered as a safe value, since [6] showed that, along the orbit(s) of Cluster spacecraft, switching from one simulation to another after 1 h of overlap produced on average a discontinuity smaller than the natural variation of the results. As shown in Section 3.1 a 4 h overlap is indeed sufficient for emulating one continuous simulation using several independent shorter simulations.

As input we use solar wind data measured by Advanced Composition Explorer (ACE) consisting of 1-min resolution plasma data and 16-s magnetic field data, available under cdaweb.gsfc.nasa.gov/pub/data/ace. Missing and invalid data are interpolated linearly from the nearest existing data regardless of gap length, except for gaps that spanned a change in calendar year. Because of this, approximately 20 days of simulations were not performed due to missing solar wind data, often during the last few days of December on several years. Figure 1 shows monthly 20, 50 and 80 percentiles of the solar wind velocity given to GUMICS-5 as input, with the median marked by a solid line and vertical bars showing the range of 20 and 80 percentiles.

The simulation extends to $\pm 248R_{E}$ in X dimension in geocentric solar ecliptic (GSE) coordinates, where $R_{E}=6371$ km is the radius of Earth. In GSE Y and Z dimensions the simulation extends to $\pm 64R_{E}$. We use a maximum magnetospheric resolution of 0.5 Earth radii ($R_{E}$) up to about 15 $R_{E}$ from Earth and a decreasing resolution further away down to a minimum resolution of 8 $R_{E}$. The boundary of the simulation on the Sun-ward side is set every minute to the position of ACE although it rarely changes from one layer of grid cells to another during one simulation. With the correct location for the source of solar wind data to within approximately $4R_{E}$, the simulation should give realistic timings for modeled phenomena within limitations of the solved MHD and electrostatic equations.

The background dipole field direction is updated every minute, the divergence of the perturbed magnetic field is removed every 20 s and the ionospheric potential solution is updated every 4 s. We save snapshots of the ionospheric state every simulated minute and of the magnetospheric state every 5 minutes, resulting in approximately 500 000 and 100 000 files per year respectively.

Most of the previously manual steps of using a GMHD model had to be automated because of the scale of the endeavor, as nearly 8000 simulations were prepared, executed, post-processed and transferred to long-term storage. In total over 2 000 000 magnetospheric output files and 10 000 000 ionospheric output files were produced with a total size of over 100 TB ($10^{14}$ bytes).

We use Finnish Meteorological Institute’s (FMI) partition on the Puhti supercomputer operated by the Finnish Center for Scientific Computing. The theoretical peak performance of Puhti is 1.8 petaflops using 682 nodes of which 238 belong to the partition dedicated for FMI. Puhti’s Lustre filesystem has 4.8 petabytes of space of which 0.6 PB is dedicated for FMI. The resources of FMI partition are shared between all users at FMI which resulted in some complications.

In practice disk space was always limited, despite a total size of almost 600 TB, and allowed the processing of only 2-4 years of simulation results simultaneously, with the rest having to be kept on a separate system designed for long-term storage. Also the computational load of the system varied widely, usually from almost full during business hours to almost empty during weekends and holidays.

In order to take full advantage of available computational resources without interfering with other users, a simulation scheduling program was developed. The scheduler was executed every hour to keep 20 simulations constantly running or queuing, using almost 10 % of total computational capacity, and to additionally start as many simulations as would fit into the FMI partition during weekends. Often the scheduling program was also run manually to start additional simulations in case significant computational capacity was available. We estimate that taking advantage of additional capacity available during weekends corresponded to constantly running 40-60 simulations representing 15-25 % of total capacity of FMI’s partition.

We solved limitations of disk space during post-processing by implementing a parallel wrapper program which downloaded results from long-term storage transparently and on-demand, and subsequently executed multiple copies of the specified serial program for post-processing results in parallel. By default this also included verifying checksums of downloaded files. A separate program was developed for initially calculating checksums, uploading the results and verifying that the upload succeeded by downloading the results into a temporary directory for verification. Deletion of results from local disk after sufficient post-processing was left for the user, as the manual steps are trivial.

While most tasks related to GUMICS-5 were automated, some were left to be done manually as they were either infrequent enough or too laborious to automate at the time. For example several times the local disk became full due to other users, meaning that most simulations and post-processing failed and had to be cleaned up and restarted. Latest version of the scheduling program is also able to restart such failed runs. On many occasions data transfers from permanent storage to local supercomputer disk failed silently or slowed down enough to require manual intervention. These cases might also be handled automatically in the future versions of data transfer and post-processing programs. Simulations also ”failed” on approximately 20 days when solar wind data was not available for the entire 28 h interval.

In order to test the convergence of GUMICS-5 results, Figure 2 shows the correlation coefficients ($CC$) between ionospheric potentials of overlapping simulations for years 2001 and 2009. For each point in time that was modeled by two simulations, the $CC$ is calculated between electric potentials of all ionospheric cells of both simulations. Thus for each minute between 20:00 and 23:59 we obtain up to 364 $CC$s each calendar year.

The 10th, 50th, and 90th percentiles of $CC$s are shown as functions of the number of minutes since beginning of the overlaps. The convergence results for other years are very similar to those shown in Figure 2.

The ionospheric results converge after about 3 h of overlap. Even after 2 h of overlap, using separate simulations instead of one continuous simulation most likely results in smaller errors than from other sources discussed in Section 4.2. As shown in [6], the dayside magnetosphere converges sufficiently even in 1 h and we assume that the magnetosphere converges faster than the ionosphere, hence the 4 h overlap used here should guarantee that the error from switching from one simulation to the next is insignificant compared to other sources. After about 1 h of overlap, convergence in the ionosphere appears to be very similar regardless of the phase of solar cycle.

In order to assess the differences to previous version of GUMICS, we compare our results to those presented in [17] which were obtained with the previous single-core version GUMICS-4, using a twice higher maximum spatial resolution of 0.25 R_E in the magnetosphere ([6]).

Figures 7 and 8 compare the magnetic field strength in magnetospheric lobe(s) to the empirical model of [7]. From GUMICS-5 we select the maximum magnetic field strength at three distances from Earth along the -X axis at any position in Z direction.

Figure 9 shows the 20th, 50th and 80th percentiles of GUMICS-5 and empirical magnetospheric lobe field strength for each calendar year between 1998 and 2020 at -20 R_E GSE X in Earth’s magnetotail.

GUMICS-5 reproduces the solar cycle time variations of lobe field well when compared to [7]. In particular, both approaches show the exceptionally strong lobe fields in 2003 and 2012. GUMICS-5 underestimates the median field intensity by approximately 10 nT but the underestimation seems to decrease further from Earth.

The empirical model produces values in a much larger range for the solar maximum in 2001 than for the minimum in 2009. In 2001 the empirical model produces values as low as 0 nT at all studied positions, while GUMICS-5 yields rather fixed minimun values: approximately 15 nT at -15 R_E, 8 nT at -20 R_E and 6 nT at -25 R_E. This is not the case in 2009, when both the empirical model and GUMICS-5 yield values above approximately 10 nT.

In [9], lobe field strength was underestimated by GUMICS-4, BATS-R-US, LFM, and for large fields by OGGCM as well.

Figure 10 compares the magnetopause standoff distance on the Sun-Earth line to the empirical model of [22]. From GUMICS-5 we select the position on the GSE $X$-axis between 5 and 30 R_E from Earth with maximum $J_{y}=\nabla\times B_{1}|_{y}$, where $B_{1}$ is the perturbed magnetic field solved by GUMICS-5. This magnetopause current is a simple and reliable way of detecting the location of the magnetopause along the Sun-Earth line in GUMICS.

Figure 11 shows the 20th, 50th and 80th percentiles of GUMICS-5 and empirical magnetopause standoff distance for each calendar year between 1998 and 2020.

GUMICS-5 reproduces the solar cycle time variations in magnetopause standoff distance quite well when compared to [22]. For example the years 2003 and 2012 associated with strong lobe fields (Figure 9) are local minima in stand-off distance in both GUMICS-5 and the empirical model. GUMICS-5 underestimates the standoff distance on average by about 0.5 R_E, i.e. by one grid cell, but the two median curves are within their error bars as characterized by the 20th and 80th percentiles.

In 2001 and to some extent in 2009 in Figure 10, the results consist of two distinct populations mostly separated by the sign and magnitude of the solar wind electric field (E_y, not shown). In the larger population IMF $B_{z}$ is mostly negative and the standoff distance is underestimated by about 0.5 R_E. In the smaller population IMF $B_{z}$ is positive and, when E_y is larger than about 3 mV/m, GUMICS-5 overestimates the standoff distance by several R_E.

The larger population above is close to the results of GUMICS-4 in the lower panel of Figure 4 in [16], while in [9] the standoff distance was also underestimated by OGGCM and LFM and slightly underestimated by GUMICS-4.

Figure 12 compares plasma sheet pressure to the empirical model of [32]. From GUMICS-5 we select the maximum thermal pressure within $\pm 10$ R_E in $Y$ and $Z$ dimensions, at three distances from Earth along the GSE $X$ axis. We only show the results at -20 R_E GSE X as at other distances the results are essentially the same except for a smaller range of values further from Earth and larger range closer to Earth.

Figure 13 shows the 20th, 50th and 80th percentiles of GUMICS-5 and empirical plasma sheet pressure at -20 R_E each calendar year between 1998 and 2020.

The modeled and empirically estimated plasma sheet pressures do not appear to be correlated. This is visible also in the solar cycle time variations, which are less consistent with each other than in the cases of lobe magnetic field and magnetopause standoff distance. The correlation also does not seem to improve by grouping results based on time of year, $Kp$, hourly average of solar wind $B_{z}$ nor $V_{x}$ (not shown).

The behavior of GUMICS-5 is similar to that of GUMICS-4 in [9], except that here the spread of input and output values is much wider and correlation is significantly smaller. This is likely due to the stark contrast in coverage of simulation inputs between approximately 200 hours of synthetic versus over 20 years of real solar wind.

Figure 14 shows the 20th, 50th and 80th percentiles of cross-polar cap potential from GUMICS-5 and the empirical models of [1, 27] on each calendar year between 1998 and 2020. Using the Alfvén Mach number and magnetopause size corrections of [27] seems to minimize the difference between GUMICS-5 and the empirical model, although their effect is small. Similarly to plasma sheet pressure, there is not much correlation in 5 minute values, nor longer averages, of CPCP between GUMICS-5 and [27] (not shown).

There does not seem to be much correlation between yearly GUMICS-5 and empirical CPCP at any percentile, although both show a decreasing linear trend from 1998 to 2020. In GUMICS-5 the behaviors of minimums and maximums of both hemispheres’ CPCP resemble each other closely.

We convert horizontal ionospheric currents from GUMICS-5 to variations of external ground magnetic field ($dB$) using the method of [33], after which we convert $dB$ to $Kp$ using the same procedure used for real $Kp$: We convert $dB$ to a local $K$ index using the FMI method described e.g. in [30, 24] and software available at space.fmi.fi/image/software. We then transform this $K$ index to the standardized index $Ks$ and finally obtain $Kp$ as described in [23], using the same stations that were used for the real $Kp$.

The left panel of Figure 15 compares the $Kp$ index obtained from GUMICS-5 to the measured one. Due to the significantly smaller amplitude of GUMICS-5-derived ground $dB$ compared to measurements, GUMICS-5 severely underestimates $Kp$. Hence, we derive GUMICS-5-specific $K$ indices for all stations used in $Kp$ by adjusting the station-specific K9 threshold such that the frequency of $K\geq 5$ is approximately 5 % ([23] and references therein). After deriving the GUMICS-5-specific $K$ index, the same procedure as above is used to compute the GUMICS-5-specific $Kp$ index. This brings the dynamic range of the GUMICS-5-specific $Kp$ to the same level as the observed one and is compared to real $Kp$ in middle and right panels of Figure 15.

In Figure 15 the minimum value of GUMICS-5-specific $Kp$ seems to be mostly between 1 and 2, contrary to measurements, especially during solar minimum in 2009 when real $Kp<1$ is seen most often. During solar maximum in 2001, real $Kp$ values are divided somewhat evenly in the range from 0 to 2, while a value of 2 is most often seen in GUMICS-5-specific $Kp$.

Figure 16 shows the yearly averages of real $Kp$, real $Kp$ calculated from GUMICS-5 and GUMICS-5-specific $Kp$ while Figure 17 shows the 20, 50 and 80th percentiles of real and GUMICS-5-specific $Kp$ each calendar year between 1998 and 2020.

In [14], horizontal ionospheric currents accounted for approximately 90% of modeled ground $dB$. However the effects of ring current and radiation belts on ground $dB$ were not studied, hence their contribution to $dB$ produced by GMHD models is less well known. As the ring current and radiation belts are neglected also here, our $Kp$ estimates are more dependent on high-latitude activity than the observed $Kp$ is.

Figure 18 compares the auroral electrojet index ($A_{E}$) obtained from GUMICS-5 to the real one and Figure 19 shows their 20th, 50th and 80th percentiles each calendar year between 1998 and 2018 when observed $A_{E}$ is available. The method for deriving $A_{E}$ is described by [5]. Determining $A_{E}$ involves removal of the quiet-day baseline, which is defined as the average of the horizontal magnetic field component over the 5 international quietest days of every month. Quiet-day determination is based on the $Kp$ index and we find relatively small differences in $A_{L}$ and $A_{U}$ indices produced by GUMICS-5 (not shown) whether we use observed quiet days or quiet days determined from GUMICS-5-specific $Kp$. The $A_{E}=A_{U}-A_{L}$ index is not affected by the determined baseline.

There are some similarities in the yearly percentiles of $A_{E}$ between measurements and GUMICS-5 but their correlation is not very good. The worse correlation of GUMICS-5 versus $A_{E}$ compared to $Kp$ is likely due to the fact that $A_{E}$ reflects more the effects of e.g. substorms, whose formation is less directly driven by solar wind and hence is more difficult to capture using GMHD models. The higher time resolution of $A_{E}$ likely also contributes to worse correlation with GUMICS-5 compared to $Kp$.

The strength of magnetospheric lobe field is modeled similarly by both GUMICS-5, GUMICS-4 and other GMHD models. In [9] lobe field was underestimated by GUMICS-4, BATS-R-US and LFM, and for large fields by OGGCM as well.

The behavior of GUMICS-5 w.r.t plasma sheet pressure is also similar to GUMICS-4 in [9], except that here the spread of input solar wind and output pressure values is much wider and correlation is significantly smaller. This is likely due to the fact that we have simulated an approximately three orders or magnitude longer time interval than was simulated in [9].

The magnetopause standoff distance is modeled similarly by GUMICS-5 and GUMICS-4, when comparing to the lower panel of Figure 4 in [16]. More specifically the behavior of both models seems similar at least when IMF $B_{z}$ is mostly southward, although this could due to the much shorter amount of simulations performed with GUMICS-4. In [9] magnetopause standoff distance was also underestimated by OGGCM and LFM and slightly underestimated by GUMICS-4.

The poor performance of GUMICS-5 in some of the comparisons can at least partly be assigned to the following causes.

Due to limitations in computational and storage capacity we changed the highest magnetospheric resolution from 0.25 R_E, used mostly with GUMICS-4, to 0.5 R_E which had been used in early GUMICS-4 simulations. In future studies, with larger computational and storage resources available and with further optimization of GUMICS-5, we hope to use the more standard value of 0.25 R_E. Lower magnetospheric resolution increases numerical diffusion which might cause GUMICS-5 to e.g. underestimate the maximum magnetospheric lobe field strength. Lower resolution also decreases the magnitude of field-aligned currents, causing e.g. lower horizontal ionospheric currents.

GUMICS-5 lacks an inner magnetosphere model which is essential for reproducing many low-latitude phenomena [21]. Incorporating an inner magnetosphere into GUMICS would require a large effort and likely should be done in the context of overall modernization of GUMICS’ solvers.

The poor agreement of plasma sheet pressure with the empirical model might be explained by the lack of an ionospheric outflow model in GUMICS-5. Adding solar activity-dependent outflow might improve plasma sheet pressure results considerably over solar cycle time scales, while adding ionospheric activity-dependent outflow could improve the results over hourly and/or daily time scales.

We have set $B_{x}=0$ in the solar wind, in order to have a divergence-free magnetic field at the solar wind boundary. In the future, a straightforward solution could be to introduce a slow, e.g. 1 mHz at most, global variation in the background $B_{x}$ within the simulated volume, in order to better take into account a changing IMF $B_{x}$ while maintaining a divergence-free field. This could perhaps be calculated automatically from the solar wind magnetic field already given as input to GUMICS-5.

Similarly to most simulations carried out with GUMICS-4 [16], we use a constant solar EUV flux that is parametrized by solar radio F10.7 cm flux. We expect the effect of changing F10.7 flux to be small based on [8].

We have not omitted simulations with gaps in solar wind data even though sufficiently long gaps could cause GUMICS-5 to miss certain short-term effects in magnetosphere or ionosphere. Most likely this does not impact results presented here but it should be kept in mind when analysing e.g. individual storms.

The relatively low magnetospheric resolution in GUMICS-5 used here, as well as the first-order MHD solvers, likely also contribute to the low correlation between GUMICS-5 and $A_{E}$ compared to the correlation with $Kp$. We also neglect the effect of ground conductivity when deriving geomagnetic indices. Its effect on the largest horizontal ground magnetic fields can be a few tens of percent [18] but such a relatively small error will likely not change our results significantly.

One or more of the following optimizations will likely be required in order to e.g. increase the highest magnetospheric resolution from 0.5 to 0.25 R_E of the simulations presented here.

We have already reduced the size of magnetospheric output files by approximately 50 % in post-processing by omitting the face-average background magnetic field stored for all faces of magnetospheric cells. This is used by the MHD solvers but is not necessary for e.g. obtaining the results presented here. The ionospheric files are about 1/100 of the size of the original magnetospheric files and most likely will not require optimization.

A further 50% reduction in size can be obtained by converting results from 64 bit to 32 bit floating point representation, either in GUMICS-5 or as an additional post-processing step. Significant further savings would require e.g. leaving out some parts of the magnetosphere most of the time, which would affect future analysis and the ability to restart the simulation. In that case, existing analysis methods could be applied directly after the simulation, before discarding parts of magnetospheric results, thereby saving space in the long run without affecting comparisons to previous results.

Optimizing CPU performance of GUMICS-5 is currently not as important as lowering the storage requirements of results, but there is at least one significant potential optimization. Currently in GUMICS-5 the magnetosphere is parallelized only using MPI, but taking advantage of OpenMP threading could further increase processing speed and reduce memory requirements.