Recreating the Horizontal Magnetic Field at Colaba During the Carrington Event With Geospace Simulations

An intriguing aspect of the famous September 2, 1859 geomagnetic disturbance (or “Carrington” event) is the horizontal magnetic (BH) data set measured in Colaba, India (magnetic latitude approximately 20°N). The field exhibits a sharp decrease of over 1,600 nT and a quick recovery of about 1,300 nT, all within a few hours during the daytime. The mechanism behind this has previously been attributed to magnetospheric processes, ionospheric processes or a combination of both. In this study, we outline our efforts to replicate this low‐latitude magnetic field using the Space Weather Modeling Framework. By simulating an extremely high pressure solar wind scenario, we can emulate the low‐latitude surface magnetic signal at Colaba. In our simulation, magnetospheric currents adjacent to the near‐Earth magnetopause and strong Region 1 field‐aligned currents are the main contributors to the large Colaba BH. The rapid recovery of BH in our simulated scenario is due to the retreat of these magnetospheric currents as the magnetosphere expands, as opposed to ring current dynamics. In addition, we find that the scenario that best emulated the surface magnetic field observations during the Carrington event had a minimum calculated Dst value between −431 and −1,191 nT, indicating that Dst may not be a suitable estimate of storm intensity for this kind of event.

geomagnetic field had its limitations (Blake, Pulkkinen, Schuck, Glocer, & Tóth, 2020). One such limitation was that measurements often needed to be taken manually, leading to long gaps between measurements (Curto, 2020). The manpower needed for high-frequency measurements to be taken continuously was often not available at these early observatories. This can be seen in the Madras Observatory, for example, which took horizontal measurements once per hour, except on Sundays when no measurements were taken (Jacob, 1884). For large geomagnetic storms in particular, even hourly measurements could lead to rapid geomagnetic variations being missed (Viljanen et al., 2014). Even where this was not a problem (such as with the continuously recording magnetograms used in British observatories; Boteler, 2006), the devices often went beyond their operational range during geomagnetic disturbances. For example, the horizontal magnetometer operated in Rome during the Carrington event had an operational range of only around 300 nT (Blake, Pulkkinen, Schuck, Glocer, & Tóth, 2020). This lack of complete and useful measurement sets makes it difficult to study the evolution of the geomagnetic field during the Carrington event.
Presently, the most important set of horizontal magnetic field measurements that exists for the Carrington event is from Colaba (near modern-day Mumbai, 18.9°N, 72.8°E;Tsurutani et al., 2003;Kumar et al., 2015). Unlike other data sets, these measurements were taken with a relatively high cadence (beginning with one measurement every hour, increasing to once every 5 min as the geomagnetic field became more disturbed), and does not appear to go off-scale at any point. This horizontal magnetic field (B H ) data set is seen to decrease by ∼1,600 nT over 2 h, then rapidly increase by ∼1,250 nT over 20 min. The rapidity and magnitude of the magnetic deflection and recovery is unique among low-latitude B H measurements in the records. The mechanism which explains this signal is uncertain, and different hypotheses have been suggested in the literature. The simultaneity of disturbances at different latitudes on the morning of September 2, 1859 has led some to suggest that the signal is at least partly due to the encroachment of currents from the ionosphere (Cliver & Dietrich, 2013;Green & Boardsen, 2006). Siscoe et al. (2006) interpreted the Colaba B H time series as having been produced by magnetospheric instead of ionospheric currents. Keika et al. (2015) suggest that the rapid recovery is due to the relaxation of the equatorial ring current due to the flow out of energetic ions. Cid et al. (2015) suggest that field-aligned currents (FACs) played a major role in the magnitude of B H .
Colaba is a low-latitude location and so the B H measurements taken there are often used to characterize the Carrington event in terms of Dst (Keika et al., 2015;Li et al., 2006;Siscoe et al., 2006), which is normally calculated from four low-latitude observatory measurements (Love & Gannon, 2009). Dst for the September 2, 1859 has been estimated using the Colaba data set as a spot measurement (Tsurutani et al., 2003) or an hourly average (Gonzalez et al., 2011;Siscoe et al., 2006). Using these Dst estimates, the solar wind conditions that gave rise to the geomagnetic storm can then be estimated using empirical relations such as the Burton equation (Burton et al., 1975), or its modified version (O'Brien & McPherron, 2000).
Studies that attempt to reproduce the large decrease and rapid recovery in Dst using these models often used very high speed (>2,000 km s −1 ), very negatively B Z (<−80 nT) and high density solar wind. For example, Li et al. (2006) required an extremely high solar wind density of almost 2,000 particles cm −3 . By extending the Burton model to include ion flow out, Keika et al. (2015) were able to reproduce a Dst of ∼−1,600 nT and a rapid recovery using a solar wind density of 100 particles cm −3 .
In this study, we use the BATS-R-US as part of the Space Weather Modeling Framework (SWMF) to simulate the Carrington event, with particular focus on reproducing the intense dip and recovery at Colaba. As inputs for this storm simulation, we used solar wind conditions that resulted in an extremely high solar wind pressure exerted on the magnetosphere. These solar wind conditions were scaled to extreme but not unreasonable values for speed, and density. We find that as the simulated magnetopause is compressed, currents adjacent to the magnetopause are the likely cause of an extreme decrease in the dayside surface magnetic field at low latitudes. As the solar wind pressure decreases, these currents retreat, allowing the rapid recovery of the low-latitude surface magnetic field. Finally, we calculate the Dst for our simulated storm event. We find that due to the unusual nature of the simulated event, Dst may not be a suitable measure for its intensity.

Simulation Setup
For this study, we use the SWMF to simulate our Carrington-scale storm event. The SWMF is a software framework for physics-based simulations of the Sun-Earth system (Tóth et al., 2005(Tóth et al., , 2012. It allows for different physics models to be combined in order to simulate across a large range of spatial scales, from the solar corona to the Earth's upper atmosphere and ionosphere. Our study uses solar wind conditions at L1 as inputs and simulates the magnetosphere and ionosphere, as well as the FACs which connect them. The model used in this study consists of the BATS-R-US coupled to the Rice Convection Model (RCM) for the inner magnetosphere, as well as the Ridley Ionosphere Model (RIM).
BATS-R-US is an MHD model that calculates the plasma conditions in the general magnetosphere domain of the simulation and uses a block-adaptive mesh grid. Blocks of higher resolution cells can be specified for regions of interest such as near-Earth locations or the magnetotail (De Zeeuw et al., 2000;Powell et al., 1999). RCM captures ring current dynamics by receiving magnetic field values and plasma moments from BATS-R-US, and returns plasma density and pressure back to BATS-R-US (De Zeeuw et al., 2004;Toffoletto et al., 2003). RIM is a height-integrated ionospheric electrodynamics model, which calculates the currents on a two-dimensional ionospheric shell . It receives field-aligned current density from BATS-R-US and communicates electric potential to RCM and BATS-R-US. The communication between each of these separate models is facilitated by the SWMF, and is heavily dependent on the ionospheric conductance (Ngwira et al., 2014). This conductance is due to solar EUV ionization (a function of solar zenith angle) and particle precipitation. This latter component is dependent on FACs in RIM, with the relation between FACs and conductance having been derived from empirical data (Welling, 2019). For our simulation, the nightside and polar cap Pedersen conductances were assumed to be 1 and 0.25 mho, respectively.
This combination of BATS-R-US, RIM, and RCM is widely used for geomagnetic storm simulations. It has been shown to perform well when statistically replicating surface dB/dt  and Dst (Haiducek et al., 2017), and has been used to study extreme geomagnetic storm events (Blake, Pulkkinen, Schuck, Glocer, & Tóth, 2020;Ngwira et al., 2013Ngwira et al., , 2014Welling et al., 2020). In order to calculate the surface geomagnetic field variations during the simulation, the Biot-Savart law is applied to the currents in each of the simulated domains: The magnetosphere, the ionosphere (Hall and Pedersen), and the gap-region FACs which connect them (Yu & Ridley, 2008;Yu et al., 2010). The total magnetic field variation at a point on the Earth's surface is equal to the sum of these contributions. As the 1859 Colaba data set was sampled at best at 5 min, the surface magnetic field for our simulation was calculated on a 1° × 1° geographic grid every 60 s.
The magnetospheric grid used for this simulation had approximately 5.89 million cells, which ranged in size from 8 R E in the far tail to 1/16 R E near the Earth. The grid used for this simulation is shown in Figure 1. The inner magnetosphere boundary of the simulation was set at 1.5 R E and the magnetopause currents were mapped from a radius of 1.8 R E to the ionosphere along assumed dipole lines. The coupling time between the different modules was set to 2 s, in order to capture the rapidly varying ionospheric currents under the extreme solar wind driving conditions. The grid resolution for RIM was 2° in longitude and 1° in latitude.
By default, the Earth's rotational axis is determined by the date and time given by the time of the simulation. In order to better recreate the conditions of the Carrington event, three additional steps were taken in the setup of the simulation. The dipole strength of the Earth in the SWMF simulation was changed to 32,000 nT. This is the approximate value given by the GUFM1 model for 1859 (Jackson et al., 2000). In addition, the position of the magnetic axis at the start of the run was set such that the magnetic North Pole was at 69.174°N, 96.757°W, its approximate position in 1859 (calculated using the International Geomagnetic Reference Field (IGRF) model, www.ngdc.noaa.gov/geomag/GeomagneticPoles.shtml). This corresponds to colatitude Θ and longitude Φ angles of 16.52° and −130.32° in Geocentric Solar Ecliptic (GSE) coordinates. Finally, the solar radio flux F 10.7 (which affects ionospheric conductance) was set to 192 solar flux units. The solar F 10.7 flux is highly correlated to sunspot number (Okoh & Okoro, 2020), and the monthly mean sunspot for September 1859 of 200, taken from www.sidc.be/silso was used to estimate this F 10.7 value.

Solar Wind Conditions
For this study, we created two synthetic solar wind scenarios to use as inputs for our simulated Carrington events. The measured horizontal magnetic field time series taken at Colaba was used as a template for the shape of each of the solar wind components: Magnetic field, velocity, particle density and temperature. Our general aim for the solar wind inputs was to produce an extremely fast solar wind with a high ram pressure, which would highly compress the magnetopause. To achieve this, each of the solar wind components were scaled to particular values so that they peaked at ∼06:30 GMT, when the Colaba B H was at its most disturbed condition. The first solar wind scenario ("Scenario 1") gradually reduces in intensity after 06:30 GMT, following the Colaba B H time series. The second solar wind scenario ("Scenario 2") immediately reverts to quiet solar wind conditions after 06:30 GMT. Finally, Scenario 1 was run again without the RCM component ("Scenario 3"). These contrasting scenarios were simulated in order to better quantify the effect of the ring current, as shall be shown in Section 3. Figure 2 shows the measured B H at Colaba in 1859 (top panel), and each of the solar wind conditions for the two scenarios (bottom panels).
For simplicity, the V Y and V Z components of the solar wind were set to zero. The incoming solar wind structure was therefore purely frontal, which maximizes the effects of its symmetric impact (Oliveira & Samsonov, 2018;Rudd et al., 2019). The V X (Earthward) component of the solar wind peaks at 2,500 kms −1 . This represents an extremely fast Earthward velocity for a coronal mass ejection (CME) at 1 AU. This is slightly slower than the extreme scenario of ∼2,700 km s −1 outlined in Tsurutani and Lakhina (2014) (and simulated in Welling et al., 2020). Tsurutani and Lakhina (2014) suggests that for such a high speed CME, interplanetary space would likely need to be "cleaned out" by a previous CME, in order to reduce drag. Manchester et al. (2005) simulated a single very fast CME that traveled to Earth in 18 h (equivalent to the estimated 17 h 40 m for the Carrington event). In order to achieve this, they needed an initial speed of BLAKE ET AL.  4,000 km s −1 at the Sun. The speed reduced to ∼2,000 km s −1 by the time the CME had reached 1 AU. Such extreme speeds have been measured in solar wind at 1 AU. In situ measurements of the July 2012 CME taken by NASA's STEREO-A (Solar TErrestrial RElations Observatory) spacecraft recorded a speed of approximately 2,500 km s −1 (Ngwira et al., 2013), and the August 1972 event had an average speed of around 2,850 km s −1 (Hoffman et al., 1975;Zastenker, 1978).
The solar wind density is set to around 25 cm −3 at the start of the simulation. It then begins to gradually increase at 04:50. Following this, the density increases rapidly to a peak of 95 cm −3 at 06:30 GMT. This peak, coinciding with the time of maximum Earthward velocity, results in a maximum solar wind pressure of ∼1,000 nPa at 06:30 GMT. While this peak density is large, similar values have been measured in the solar wind before (Tsurutani et al., 2017).
For simplicity again, the solar wind B X and B Y components were also set to zero. A maximum intensity for the magnetic field of the solar wind inputs was set at 118 nT. This value is calculated from the empirical relationship between velocity and B z of magnetic clouds at 1 AU recorded by Gonzalez et al. (1998) (top panel), and the solar wind conditions used to drive the two simulation scenarios (bottom five panels). Solar wind conditions for Scenarios 1 and 3 are in bold red, and for Scenario 2 are in blue. These conditions are southward interplanetary magnetic field (B Z ), velocity (V X ), particle density, pressure and temperature. The shaded blue region is the time-period simulated. The Colaba data were digitized from Figure 3 of Tsurutani et al. (2003).
Although we note here that this relationship is derived from a limited data set with peak B intensities of <40 nT. With our maximum V X equal to 2,500 kms −1 , and our B X , B Y set to 0, this gives a minimum B Z equal to −118 nT. This is slightly smaller in amplitude than the "perfect" CME scenario given by Tsurutani and Lakhina (2014) (which had a B Z = −127 nT). Lastly, the temperature was scaled to a maximum of approximately 6 MK. This is similar to the maximum temperature of the July 2012 fast CME. Once all of the solar wind components were constructed, small levels of noise drawn from normal distributions with standard deviation proportional to the amplitude of each of the solar wind components were added, in order to mimic the small-scale fluctuations in the solar wind.

Simulation Results
The simulations were each run for approximately 12 h (from 04:10 GMT to 16:00 GMT). The top rows in Figure 3 show the normalized total current density (J) in the near-Earth GSE XZ-plane (top) and the pressure in the X − Y plane (middle) for three different times during the Scenario 1 simulation. These times correspond to pre-storm (column A), when the magnetopause was at its most compressed (column B) and during the main phase of the storm (column C). The black and white lines are open and closed field lines in the XZ-plane. The red lines are the closed field lines whose footprints are closest to the poles. The location of the magnetopause standoff distance is calculated from these lines. The third row shows the magnetopause standoff distance for Scenario 1 (red line) and Scenario 2 (blue line) for the duration of the simulations. The bottom row shows the SWMF-derived Dst for each storm scenario. Vertical dashed black lines correspond to the times shown in the top two rows.
At the start of the simulations, the magnetopause is situated at approximately 8 R E . The position of the magnetopause was forced to a minimum standoff distance of 2.26 R E at 06:30 GMT for both Scenarios 1 and 2 by the solar wind. This corresponds to the time when the input solar wind pressure was at its highest (at ∼1,000 nPa). This extremely close standoff position is well within the orbit of many satellites, and would expose them to the solar wind environment, causing damaging surface charging (Koons & Fennell, 2006). In Scenario 1, the subsolar magnetopause is expanded to a position >5 R E by 07:05 GMT, and it remains in the range 5-7.5 R E for the rest of the simulation. For Scenario 2, as the solar wind immediately reverts to pre-storm conditions, the magnetopause rapidly expands and returns to a position >10 R E by 07:40 GMT. After this, the magnetopause slowly returns to around 8 R E . For both scenarios, the minimum Dst calculated by the SWMF is −989 nT. This Dst is discussed further in Section 4.3.

Simulated B H at Colaba
As mentioned above, the surface magnetic field was calculated for every minute of the simulation on a 1° × 1° geographic grid. The magnetic field at each point on this grid is the sum of magnetic contributions from 1. All of the currents in the magnetospheric domain of the simulation 2. FACs in the region between the magnetospheric and ionospheric domains 3. The ionospheric domain of the simulation (Hall and Pedersen currents) Panels (a), (b), and (c) of Figure 4 show the historical (green) and SWMF simulated ( , blue). Also shown are the dynamic solar wind pressures (P SW , dotted green). For Scenario 1, at the peak of the B H depression (06:30 GMT), currents in the magnetosphere can be seen to contribute most to the B H at Colaba, with about −1,050 nT. Next in importance are the FACs, contributing around 720 nT. In comparison, the contribution from the ionosphere is largely insignificant, contributing around 50 nT.  Panels (e) and (f) show the northward contribution from FACs and magnetosphere respectively. The color bars for these plots were limited to ±1,000 nT to better highlight magnetic conditions at the dayside.
Colaba is marked with a yellow star in each panel, and the solar equator is marked with a white line. From Panel (a), the largest geomagnetic disturbances are generally in the regions around the geomagnetic poles (North of Canada and Antarctica at around 90°E). In addition, the low-latitude region around noon (at ∼90°E) sees smaller, although still significant geomagnetic disturbances. Although being too southerly to be affected by ionospheric currents (Panel (b)), Colaba is situated such that it sees significant negative B N field contributions from both the FACs (Panel (e)) and the magnetosphere (Panel (f)), with the magnetosphere contributing the most in the region just north of the equator. These will each be discussed in more detail in the next sections.

FAC Contribution to Surface B
FACs are mapped along assumed dipole field lines between the magnetosphere and ionosphere domains of the SWMF simulation. From Panels (d) and (e) of Figure 4, the amplitude of the contribution from FACs to the B H at Colaba follows the amplitude of the dynamic solar wind pressure acting upon the magnetosphere for both solar wind scenarios. As the solar wind pressure increases, more energy is transported from the magnetosphere to the ionosphere, resulting in magnetic field contributions to the surface magnetic field.
The top panels of Figure 6 show the ionospheric radial current density in the Northern and Southern hemispheres at 06:30 GMT for Scenario 1. These show the footprints of the FACs in the ionosphere. At this time, the Region 1 currents dominate, and can be seen to extend significantly equatorward on the dayside. Region 1 currents, which flow at higher latitudes than Region 2 currents, connect to the outer magnetosphere, and given the proximity of the extremely compressed magnetopause and the position of the last closed field lines at 06:30 GMT (Panel (b)   The middle and bottom panels of Figure 6 show the ionospheric Hall and Pedersen conductances for the same time. Both conductance distributions show complicated nightside structures, which are likely due to the Region 2 FACs distributions. Panels (c) and (e) of Figure 5 shows the magnetic effect of FACs on the ground. Colaba is situated nearly equidistant between the two sets of Region 1 currents in each hemisphere. The two Region 1 current systems contribute a very negative component in terms of B N (Panel (e), Figure 5). The 720 nT signal at Colaba from the FACs at 06:30 GMT appears to be due to a combination of FACs whose ionospheric footprints are in both hemispheres North and South of Colaba.

Magnetosphere Contributions to Surface B
The extreme solar wind pressure exerted on the magnetosphere forces the magnetopause to within 8,700 km of the Earth's surface at 06:30 GMT in our simulation. At this time, it is nearly directly above Colaba, which lies at approximately noon. The magnetopause flows east around the Earth, and so this setup should result BLAKE ET AL.
10.1029/2020SW002585 9 of 21   Immediately adjacent to the compressed magnetopause are regions of westward flowing current. This is true in the XZ-pane, and XY-plane. The current in these regions are less intense than that of the magnetopause, however, they cover a larger volume than the magnetopause, and can be seen to extend closer to the Earth in the XZ-plane. Panel (f) in Figure 5 shows north-south component of the magnetic field from the magnetospheric domain. The large negative B N region over Asia is not symmetric around the solar equator, but is instead shifted northward by approximately 20°. This corresponds to the position of the magnetopause, and the westward flowing current systems adjacent to the magnetopause.
Next, we look at the role of the ring current on the magnetic field at the dayside of the Earth. For our standard SWMF setup (BATS-R-US with RCM and RIM), the contribution from the magnetosphere to the ground magnetic field perturbation includes the effect of the ring current. While the ring current contribution itself is not separable from the total magnetospheric contribution in the model, we can use plots of the near-Earth pressure and the contrast between Scenario 1 and Scenario 2 to infer its importance. In addition, we can contrast Scenario 1 to Scenario 3, which was run without the RCM.
RCM solves for the bounce-averaged drift convection of ions and electrons in the closed field line region. As the magnetopause is heavily compressed on the dayside of the Earth in our scenarios, this reduces the volume where ring currents can be trapped. For this case, the higher energy portion of the ring current will drift further outward due to conserve the third adiabatic invariant. In addition, particles drifting from night to day are more likely to be lost. These effects should cause the depletion of the ring current on the dayside however during the most compressed part of the storm, there is a partial ring current concentrated in the pre-midnight sector, which should mitigate this effect.
On the other hand, strong storms have been shown to result in deeper penetration of energetic particles to lower L-shells due to the strong convection electric field at this time (Li et al., 1993;Zhao & Li, 2013   portion of the ring penetrating to very low L-shell could still exist on the dayside and would be closer to Earth than normal. Figure 3 shows a ring current concentrated at low L compared to normal conditions. Figure 3 shows the near-Earth pressure for three time-periods for Scenario 1. Up to 06:30, there is an intense build-up of pressure at the nightside of the Earth. After this time, as the magnetopause rapidly retreats then remains beyond 7.5 R E , the pressure wraps westwardly around the Earth until it envelops the Earth by around 08:00. This is the ring current that remains energized throughout the rest of the simulation as the solar wind continues to deposit energy into the magnetosphere. Panel (d) of Figure 4 shows this sustained B H contribution to Colaba from the magnetosphere hours after the massive dip and recovery.
In Scenario 2, the solar wind conditions immediately revert to solar-quiet conditions at 06:30 GMT, and the magnetopause rapidly expands (bottom row, Figure 3). This corresponds to a marked recovery in the for Scenario 2 between 06:30 GMT and 06:35 GMT is due to the expansion of the magnetopause, and the retreat of the previously identified westward current lobes adjacent to it. In Scenario 1, we see a similar recovery of around 500 nT between 06:30 GMT and 07:30, which corresponds to the slower expansion of the magnetopause (Figure 3).
Finally, we look at the B H at Colaba for the Scenario 1 simulation without RCM (Scenario 3, further images and movies of this simulation are given in the supporting information). As could be expected, for this simulation there was no significant build-up of charged particles in the near-Earth magnetosphere, and no ring current formation. Despite this, dip and recovery in B H similar to Scenario 1 can still be seen at Colaba (Panels (c) and (f) in Figure 4). In this case, the magnetospheric currents and FAC contributions to B H at Colaba are seen to increase in amplitude up to 06:30 GMT, then sharply decrease, mirroring the solar wind pressure on the magnetosphere and the motion of the magnetopause. The magnetospheric contribution can also be seen to reduce by 650 nT between 06:30 and 07:30 GMT, where it remains for the rest of the simulation.
Between the three different scenarios in this study, we can infer the behavior of the simulated ring current, and its role in the dip and recovery of the Colaba B H . The dip (up to 06:30 GMT) appears in all three simulations, indicating that it is not due to a build-up of the ring current. However, the magnetosphere clearly contributes more to B H at Colaba when RCM is included as part of the simulation. When the magnetosphere was at its most compressed, there is a difference of around 250 nT between the MAG H B of the Scenarios 1 and 3 simulations. This is possibly due to the nightside buildup of particles at this time (Panel (b) of Figure 3). Despite this, the magnetosphere of Scenario 3 is still able to produce around −800 nT by 06:30 GMT at Colaba without a ring current. The recovery also appears not be due to a ring current. Between 06:30 and 07:30, both Scenarios 1 and 3 show a recovery of around 500 nT in MAG H B . The difference between the two is that there is a continuous B H contribution from the magnetosphere in Scenario 1 (where a ring current is formed then kept energized by the solar wind driver), whereas this is absent in Scenario 3. in Scenario 2 shows how the simulated ring current would decay in the complete absence of solar wind drivers. This would appear to be too slow to account for the rapid recovery seen at Colaba (as it is simulated in our setup). While not modeled in the version of RCM used in our simulations, a critical loss mechanism in the ring current is charge exchange with geocoronal neutral hydrogen. In the present simulations, ring current populations are being injected deeply to low L-shells where the geocoronal density is high and, equivalently, the charge exchange loss rate will be high (Ilie et al., 2013) An inclusion of realistic charge exchange loss terms would drive additional losses not captured in the current model version, potentially improving the comparison to the Colaba magnetometer station.

Discussion
The extreme geomagnetic storm event simulated here offers possible insight for three questions regarding the Carrington event: 1. What mechanism could have caused the extreme dip and rapid recovery of B H at Colaba? 2. What were the solar wind conditions that could have driven such an event? 3. What was the Dst of the event?
We will discuss each point here, followed by a brief discussion of where the simulation failed to reproduce Carrington event measurements.

Colaba B H Dip and Recovery
For our simulation, the extreme dip and rapid recovery of B H at Colaba (and indeed the low-latitude region of Earth at solar noon) was a combination of near-Earth magnetospheric currents and FACs. As the solar wind pressure increased to a peak of around 1,000 nPa, the magnetosphere was compressed, and the magnetopause was forced extremely close to the Earth. Adjacent to the magnetopause were symmetric regions of westward current. These lobes of current led to a very negative B N signal at low-latitudes at the dayside of the Earth. Coupled with current contributions from Region 1 FACs in the simulation (which connect the magnetosphere and ionosphere domains), the B H signal at Colaba was negative enough to match the −1,750 nT dip measured in 1859.
The rapid recovery of B H in our simulation appears to be mostly due to the expansion of the magnetopause (as the solar wind conditions are relaxed), taking the magnetospheric lobes of westward current further away from the Earth. As the solar wind pressure reduces and the magnetopause retreats, the negative B N contribution from FACs similarly decreases as well. For Scenario 1, the combination of these factors means that the recovery of B H at Colaba between 06:30 GMT and 07:00 GMT (approximately + 1,100 nT) broadly matches the rapidity and amplitude of the recovery measured at Colaba in 1859.
As shown in the previous section, the role of the ring current in the rapid dip and recovery at the dayside of the Earth is minimal. Although there was a nightside build-up of near-Earth particles which did contribute to the negative B H at Colaba when RCM was included (and the ring current simulated), an approximation of the dip and recovery was reproduced when RCM was not included in the simulation. After around 07:30 GMT, the contribution from the magnetosphere in Scenario 1 is due almost entirely due to ring current, similar to normal large-scale geomagnetic storm events at low-latitudes. Without RCM, the signal at Colaba is underestimated after 07:30 GMT.

Solar Wind Driver
The solar wind conditions which drove our simulation were chosen to exert an extreme dynamic pressure upon the magnetosphere, in order to drive the scenario described above and emulate the B H at Colaba in 1859. We chose this measured B H time series as the template for our synthetic solar wind time series, so that the proximity of the magnetopause to Earth would mirror the dip and recovery of the historical time series. For simplicity, the other solar wind conditions (B, n, T) also used the same template.
The peak values used in this simulation for the solar wind velocity, density and B Z value are large, but not so that they are impossibly large. Maximum Earthward solar velocity equal to 2,500 km s −1 has been measured before in the July 2012 CME (Ngwira et al., 2013). The solar wind density peaks at 96 cm −3 , which is unusually high, but lower than extreme densities which have been measured before in CME filaments (Tsurutani et al., 2017). The minimum solar wind B Z was empirically scaled with the velocity to give −116 nT. This is an extreme value for a CME, only slightly smaller than the B Z of Tsurutani and Lakhina (2014)'s "perfect" CME scenario. The combination of solar wind velocity and B Z gives a maximum interplanetary electric field equal to 290 mVm −1 , approximately 50% larger than the estimation given by Tsurutani et al. (2017). The solar wind conditions used in this simulation could be considered an incremental improvement on previous simulations of Carrington-type storms, as aspects of the storm were successfully recreated without requiring unrealistically high solar wind densities (as in Ngwira et al. (2014)). However, there remains many aspects of the solar wind conditions used that can be improved upon. This is evident when looking at measured magnetic field data at sites other than Colaba (this is explored further in Section 4.4).
Each of the solar wind components used to drive our simulation increase and decrease in intensity in concert. This adds to the geoeffectiveness of the solar wind (e.g., solar wind pressure, which is a product of V X and n, or the magnetic reconnection rate which is a product of V X and B Z ). In reality, CMEs have a filament structure where the different components will peak at different times. Also effectively neglected in our simulation are the small-scale variations in the solar wind. Small-scale noise was added to the components, however, the amplitude of this noise was arbitrarily chosen. Future simulations will include realistic noise and CME dynamics based on statistical analyses of historical solar wind conditions. Finally, little attention was paid to the solar wind conditions at the arrival of our CME, where we instead opted for a gradual increase in the intensity.

Dst of the Event
The Dst index is an hourly global index calculated from four midlatitude geomagnetic stations. It is one of a number of geomagnetic indices in use (such as AE, Kp, Ap, aa, etc.), each of which are designed to measure different aspects of geomagnetic disturbances. The Dst index is often used to determine the occurrence and length of geomagnetic storms, and is used as a proxy for the strength of the symmetric ring current. However, Dst is affected by other storm-time phenomena, and reducing all of the global storm dynamics to a single hourly Dst value (or its minutely equivalent, SYM-H) is at best a simplification (Borovsky & Shprits, 2017).
Despite this, the minimum Dst value for a storm event is commonly used as a simple indicator of storm severity. Dst is also an often used metric in statistical studies to estimate the intensity of extreme but infrequent geomagnetic events (Chapman et al., 2020;Morina et al., 2019;Riley, 2012). Recently, multiple large pre-1957 storms have been studied and had quasi-Dst values calculated Hayakawa et al., 2020;Love et al., 2019aLove et al., , 2019b. Finally, as outlined in Section 1, the Dst of the Carrington event itself has been subject to much inquiry. For these reasons, we discuss the calculated Dst of our simulated event.
As shown above, the extremely impulsive scenario outlined in this study resulted in a very asymmetrical B field at low-latitudes particularly when the magnetopause was at its most compressed. In addition, the unusual signal at the dayside was not caused by the ring current, which the Dst series was devised to measure. This poses a problem when attempting to classify the strength of the storm using the Dst index. In the SWMF, the Dst value is calculated by summing the magnetic contribution from the magnetosphere domain and the ionosphere to a point in the center of the Earth using the Biot-Savart law (Liemohn et al., 2018). In reality, the Dst index is calculated from four low-latitude geomagnetic observatories with the following equation: where Dist i is the "local disturbance" at observatory i (the hourly B H variation minus baseline and solar-quiet signal), and the λ i term is the magnetic latitude of the ith observatory (Love & Gannon, 2009). The four sites currently used to calculate Dst are located at Hermanus, Kakioka, Honolulu, and San Juan. The top panel in Figure 8 shows the location of these four observatories as green diamonds, overplotted on the maximum calculated B H at all points for the duration of the simulation.
Substituting the SWMF B H time series at each of these sites into Equation 2 instead of Dist, a quasi-Dst value was calculated. The magnetic latitude of each observatory (λ) was calculated for 1900 (as this is the earliest date for which magnetic coordinates can be calculated using the IGRF model). Although the distance between the GUFM calculated 1859 and 1900 North Pole is approximately 140 km, this would minimally affect the Dst calculation for Equation 2. The resulting Dst for these four sites is shown in the bottom panel of Figure 8 in green, and for the first 7-h of simulation, the minimum Dst was −449 nT. Next, the longitudinal position of each of these sites was shifted in 1° increments, and this Dst calculation was repeated, giving 360 Dst time series. The positions of these shifted sites which resulted in the minimum calculated Dst are shown as red stars in the top panel of Figure 8. When the sites were located at these positions, the minimum calculated Dst was −1,303 nT. This time series is shown in red in the bottom panel of Figure 8.
As can be seen in this figure, the longitudinal position of the Dst recording sites has a significant effect on the magnitude of the calculated Dst, particularly during the first few hours of the disturbance. In our impulsive storm scenario, the high pressure solar wind forces the magnetopause and magnetospheric currents very close to the dayside of the Earth, leading to a localized low-latitude B H signal. As pointed out in Cid et al. (2015), for B H spikes that are large in amplitude but spatially localized, a calculated Dst index may not accurately reflect the intensity of a storm.
Using the process outlined above, we have a range of Dst-equivalent values with which we can describe our The SWMF calculated Dst for the simulation (minimum Dst = −989 nT) and the maximum calculated Dst with the recording sites shifted longitudinally (minimum Dst − 1,303 nT) are within or near the spread of Dst estimates for the Carrington event in the literature. These are between −850 and −1,160 nT for empirical estimates and between −625 and −1,160 nT for modeled estimates (see Cliver and Dietrich (2013) and references within). This is unsurprising given that the previously calculated Dst values use the 1859 Colaba data set, which our simulation emulates.
The Dst of a storm event can be used as an input for calculating the drag effects on satellites during geomagnetic storms (Emmert, 2015;Prölss, 2011;Zesta & Huang, 2016). Oliveira et al. (2020) used recreated Dst estimates for different large-scale historical superstorms, and calculated the magnitude of satellite drag effects. Given the unusual nature of our Carrington event simulation (where the calculated Dst is particularly not a good proxy for the ring current), we believe that our Carrington event simulation is not suitable for these calculations.

A Realistic Description of the Carrington Event?
The main focus of this study was to emulate the extraordinary B H data set at Colaba using a physics-based simulation. The solar wind drivers used in the simulation were tailored to force the magnetopause close to the Earth and then rapidly move away. This resulted in the intense dip and rapid recovery of B H at Colaba, and was achieved using extreme but not unreasonable solar wind maxima conditions.  Figure 4, we can interpolate B H at these locations and compare with the measured values. The comparison between measured and simulated B H at these three sites is shown in the top row of Figure 10.
It is clear that the measured B H time series at these three sites are not recreated as well as for Colaba in the simulation. In the case of Rome, there was a reported off-scale deviation of approximately 3,000 nT between 06:00 and 06:30 GMT. Allowing for timing errors, the amplitude of the simulated dip at Rome is less than 1,000 nT, a significant underestimate of the total ΔB H . Blake, Pulkkinen, Schuck, Glocer, and Tóth (2020) suggests that the extremely large deviation measured at Rome indicates that it was beneath the auroral oval on the morning of the 2 September 1859. From Panel (b) of Figure 5, the magnetic effects of the ionosphere can be seen to reach low latitudes equivalent to Rome (near Central America on the nightside for example), but it does not appear to reach Rome for our simulation. In the case of St. Petersburg and Barnaul, measurements were taken only once every hour. While this makes it difficult to directly compare to the 1-min simulated data (see Viljanen et al. (2014)), at both sites there are clearly large positive B H deviations in the time series toward the beginning of the storm (+500 nT at St. Petersburg and +720 nT at Barnaul). Neither of these positive deviations are recreated in our simulation. Indeed, the simulated B H at each of the three sites follow much the same shape, dipping and recovering much like the simulated B H at Colaba. However, as shown in the bottom panel of Figure 10, the ratios of the B H contributions from the magnetosphere, ionosphere and FACs are different for each site.
While our simulation can emulate the Colaba B H data set well, from these three examples, it is clear that our simulation is not a good global recreation of the Carrington event. This is to be expected, particularly given the simplifications in our chosen solar wind profile (see Section 4.2), as well as the physical limitations and caveats of the SWMF model, especially when simulating such an extreme event.
The particular configuration of SWMF that we have used for this study (BATS-R-US, RCM, and RIM) have been shown before to be capable of reproducing surface magnetic fields for other storm events (Welling, 2019). However, there are a number of aspects of the simulation that could be appraised for future simulations. The Biot-Savart approximation used in our simulations is a free-space implementation, which does not account for conducting bodies or boundaries. For a simulation with significant short-lived magnetospheric current structures near the Earth (such as the event in this study), the application of the Biot-Savart approximation could be investigated.
The resolution of the cells closest to the Earth can also affect the calculated surface magnetic fields. While our simulation used cells near the Earth sized at 1/16 R E , it is unclear if these were small enough to adequately resolve phenomena such as bursty bulk flows and other substorm dynamics. While our simulation was primarily concerned with the ground magnetic contributions from large-scale structures, this is an avenue that should be explored in the future for more physically accurate simulations, particularly when analyzing the ground magnetic field at auroral latitudes.
The RCM used for this simulation used a fixed composition of 0.9:0.1 for H+:O+. Simplifications such as these can be changed for future simulations.

Conclusion
Using the SWMF, the quick dip and rapid recovery of the measured B H at Colaba during the Carrington event were recreated during a large geomagnetic storm scenario. This was achieved by using an extremely impulsive (∼1,000 nPa) solar wind driver, which forced the magnetopause to within 2.3 R E of the Earth. The main contributions to the extremely negative B H at Colaba were due to the proximity of westward electric currents adjacent to the magnetopause, as well as from strong Region 1 FACs. As the solar wind pressure abated and the magnetopause moved away from the Earth, the quick recovery of B H at low-latitudes was primarily due to the motion of the westward electric currents in the near-Earth dayside magnetosphere. For our simulation, the ring current did not appear to play a significant role in the Colaba dip and recovery.
The solar wind amplitudes needed to drive such a scenario were extreme, but not unrealistic, with maximum V X = −2,500 kms −1 , B Z = −127 nT and maximum n = 95 cm −3 . These peak values for V X and n are smaller than has been measured before in the solar wind. The Dst calculated for this event by the SWMF had a minimum of −989 nT. However, given the localized nature of the low-latitude dayside B H during the storm, the longitudinal location of Dst recording sites would have a marked effect on the calculated Dst.
Depending on the timing of this storm, minimum calculated Dst could range from −457 nT to −1,087 nT.
While the storm scenario in this study was able to reproduce the B H at Colaba, it performed significantly less well when recreating the B H at other geomagnetic observatories. It would be unreasonable to declare the geomagnetic storm scenario in this study as being an accurate replication of the magnetospheric and ionospheric conditions during the Carrington event. Rather, it highlights a possible mechanism for the unusual dip and recovery of B H at Colaba using (somewhat) realistic solar wind conditions. That is to say, the simulated low-latitude dayside B H which emulated the Colaba 1859 data set is heavily dependent on the location of the magnetopause during this extreme storm scenario.
This study should be viewed as an initial step in simulating the Carrington event, where a tailor-made simulation was tweaked to reproduce the extraordinary Colaba B H data set. By using different combinations of SWMF modules parameters and by particularly varying the solar wind inputs, we may be able to better recreate and understand the global processes that occurred during the extraordinary 1859 Carrington event storm.

Data Availability Statement
The Colaba B H data were digitized from Tsurutani et al. (2003) using the WebPlotDigitizer software (https:// automeris.io/WebPlotDigitizer/). The Barnaul and St. Petersburg B H data were taken from Nevanlinna (2008). Rome B H data were taken from Blake, Pulkkinen, Schuck, Glocer, and Tóth (2020). Data sets used to produce the images in this study can be found at https://doi.org/10.5281/zenodo.4534660. The complete run folders for each SWMF simulation can be requested from the Community Coordinated Modeling Center (https://ccmc.gsfc.nasa.gov/index.php).The results presented in this study rely on data collected at magnetic observatories. The SWMF simulations were run on the Discover supercomputing cluster at NA-SA's Center for Climate Simulation (HEC-SMD-18-1749). Data were also obtained from the SuperMAG database (http://supermag.jhuapl.edu/info/?page=faq).