Employing the Coupled EUHFORIA‐OpenGGCM Model to Predict CME Geoeffectiveness

EUropean Heliospheric FORecasting Information Asset (EUHFORIA) is a physics‐based data‐driven solar wind and coronal mass ejections (CMEs) propagation model designed for space weather forecasting and event analysis investigations. Although EUHFORIA can predict the solar wind plasma and magnetic field properties at Earth, it is not equipped to quantify the geo‐effectiveness of the solar transients in terms of geomagnetic indices like the disturbance storm time (Dst) index and the auroral indices, that quantify the impact of the magnetized plasma encounters on Earth's magnetosphere. Therefore, we couple EUHFORIA with the Open Geospace General Circulation Model (OpenGGCM), a magnetohydrodynamic model of the response of Earth's magnetosphere, ionosphere, and thermosphere to transient solar wind characteristics. In this coupling, OpenGGCM is driven by the solar wind and interplanetary magnetic field obtained from EUHFORIA simulations to produce the magnetospheric and ionospheric response to the CMEs. This coupling is validated with two observed geo‐effective CME events driven with the spheromak flux‐rope CME model. We compare these simulation results with the indices obtained from OpenGGCM simulations driven by the measured solar wind data from spacecraft. We further employ the dynamic time warping (DTW) technique to assess the model performance in predicting Dst. The main highlight of this study is to use EUHFORIA simulated time series to predict the Dst and auroral indices 1–2 days in advance, as compared to using the observed solar wind data at L1, which only provides predictions 1–2 hr before the actual impact.


Introduction
The diverse socio-economic and health impacts of solar-terrestrial relations have drawn attention towards geomagnetic storms and space weather forecasting (FREng, 2013;Oughton et al., 2019;Pilipenko, 2021;Zenchenko & Breus, 2021;Buzulukova & Tsurutani, 2022).The solar wind causes spatiotemporal variability in the geospace environment, especially the magnetosphere and the ionosphere (Chapman & Ferraro, 1931;Akasofu & Chapman, 1963;Akasofu, 2021, and references therein).Earth's magnetosphere provides natural protection from the solar phenomena.However, a southward interplanetary magnetic field (IMF) caused by solar transients like coronal mass ejections (CMEs) can result in a rapid magnetic reconnection with Earth's northward pointing magnetic field in the day side of the magnetosphere causing geomagnetic storms (Chapman, 1918;Dungey, 1961;Bothmer & Daglis, 2007).A geomagnetic storm is characterized by three phases.First is the sudden commencement, when Earth's magnetic field is suddenly compressed, second is the main phase when the magnetic field rapidly decreases, and finally in the recovery phase the magnetic field strength returns back to normal (Chapter 4; Bothmer & Daglis, 2007).Moreover, when the two reconnection sites (the one on the day side and the other one on the night side) are sufficiently separated, different parts of the magnetosphere respond to the storm at different times, leading to intermittent storage and release of energy.This results in a phenomenon called magnetospheric substorm (Akasofu, 1968).These auroral/magnetospheric substorms are also linked to the ring current development during the main phase of the storm (Akasofu, 2020;Alberti et al., 2022).The substorms also interfere with GPS communication and impact electric power systems (Skone & de Jong, 2000;Boteler, 2001, and references therin).
The process of using the ground magnetic field observations directly in real-time, to evaluate geomagnetic indices (like the Disturbance storm time (Dst) index, auroral indices or A-indices -AU, AL, AE, and Kp index) is called 'nowcasting' or 'short-term forecasting'.Nowcasting is important for situational awareness, i.e., for verifying whether the occurrence of power outages or satellite malfunctions is a consequence of a geomagnetic storm or if they are caused by terrestrial or man-made accidents or activities of other socio-political agencies.However, the geomagnetic indices must be forecasted for mitigation purposes.The standard forecasting methodology involves predicting the indices with the near-Earth near-real-time (NRT) solar wind data, i.e. using the solar wind conditions at L1 to compute the geomagnetic indices at Earth (Keesee et al., 2020;Smith et al., 2021Smith et al., , 2022)).With NRT data, the earliest warning of impact at Earth that we can obtain at present times, is from the satellites at L1, unless there is a serendipitous multi-spacecraft lineup before 1 au.This gives us less than two hours (20-90 minutes; Wintoft et al., 2017) to take measures against these storms which is substantially insufficient for prompt and effective mitigation for the spaceand ground-based technological infrastructure.The ideal lead time for the best mitigation for the power grid operators is 2-3 days and for the aviation industry is 24 hours (Joint Research Centre (JRC) Science for Policy Report, https://publications.jrc.ec.europa .eu/repository/handle/JRC104231,Krausmann et al., 2016).This major limitation of a short forecast time window associated with using NRT observations, motivates us to employ modeled solar wind data (predictions) in computing the geomagnetic indices to increase the forecast time window.
Various state-of-the-art empirical (e.g., HUXt, Barnard and Owens (2022); DBM, Vršnak et al. (2013), etc), and MHD models (e.g., ENLIL, Odstrcil (2003); MAS, ?(?); SUSANOO, Shiota andKataoka (2016), AWSoM, van der Holst et al. (2014), EUHFORIA, Pomoell and Poedts (2018); etc) have been created to study CME propagation and predict their arrival time.Heliospheric models like EUHFORIA provide an additional advantage of employing magnetized CME models (Verbeke et al., 2019;Maharana et al., 2022) for improving the prediction of magnetic field components (especially the z component of IMF) in addition to plasma properties.The first step in the realistic prediction of CME arrival time at Earth is modeling as accurately as possible the ambient solar wind conditions through which the CMEs traverse.The second step involves constraining the initial CME parameters (geometric and magnetic field) from pre-and post-eruption observations (Scolini et al., 2019).The final step is driving the 3D heliospheric simulations with the modeled ambient solar wind conditions and inserting the CMEs as time-dependent boundary conditions.Ideally, with the best and fullest efforts involving the maximum human and machine resources, the arrival time prediction of a CME event must be accomplished within 1-2 days from the launch of the CME from the Sun.If achieved, the magnetosphere models can be employed to make the geomagnetic index forecast 1-2 days in advance of the actual recording of the geomagnetic storm at Earth.This calculation assumes an average actual CME arrival time of 3-4 days (Iwai et al., 2021).The increased lead time of 1-2 days is much better in comparison with the above-mentioned 20-90 minutes warning time, and can potentially help in monitoring and planning mitigation strategies more efficiently.For the prediction of even faster CMEs, further improvements are required in individual models in the chain.
The quantification of the comprehensive impact of the solar wind on the geospace, the connected layers of magnetosphere, ionosphere, and thermosphere (MIT), requires the consideration of interactions with the incoming plasma and magnetic field.The complexity of the geospace, because of its multiple coupled components, added to the non-linearity and time dependence of the physical processes, demands a numerical approach for its modeling.Three-dimensional magnetosphere models like Open Geospace General Circulation Model (OpenGGCM, Raeder et al., 1998), GUMICS (Janhunen et al., 1996), BATSRUS (Powell et al., 1999) -implemented as the Space Weather Modeling Framework, (SWMF, ?, ?), and PPMLR-MHD (Y.Q. Hu et al., 2007) etc. (see the review by Wang et al. (2013)) complement single point-based observations by satellites or ground-based instruments to understand the complex physics.In a study assessing different types of magnetospheric models, Rastätter et al. (2013) reported that the empirical models (analytic or iterative formula or neural network-based algorithm without involving calculations of the intrinsic energy flow in the geospace) perform well in Dst prediction, in general.The magnetosphere models that could match the accuracy of the empirical results are the global MHD models of the magnetosphere that include the inner magnetosphere dynamics, i.e., the ring current model or the standalone ring current models with well-defined boundary conditions.OpenGGCM performed average or poorly relative to models like BATSRUS and CMIT in the analysis by Rastätter et al. (2013) because of the absence of a ring current model.However, the latest version of OpenGGCM (v5.0.ccmc) is coupled to a ring current model, Rice Convection model (RCM; Toffoletto et al., 2003) capable of computing the kinetic ring current (Cramer et al., 2017).In addition, OpenGGCM is coupled to an ionosphere and thermosphere model, which makes it a global model of the geospace.
The physics of CME eruption, propagation, and its interaction with the magnetosphere happens at different spatial and temporal scales.Coupling of models from the Sun to Earth (Luhmann et al., 2004;Tóth et al., 2007) has become a popular 'trend' in space weather modeling as it enables faster and more efficient space weather predictions at Community Coordinated Modeling Center (CCMC,https://ccmc.gsfc.nasa.gov/)and Virtual Space Weather Modelling Centre (VSWMC,https://esa-vswmc.eu/) (Poedts, 2019).Combining models working at different spatiotemporal scales conserves the physics at each stage and also consumes less time than running a single model from Sun to Earth.The objectives of this study are two-fold: (1) To demonstrate the coupling of 3D MHD models of the heliosphere (EUHFO-RIA) and magnetosphere (OpenGGCM), and (2) to emphasize the possibility of forecasting geomagnetic indices with a chain of physics-based models instead of nowcasting or forecasting with L1 observations that provide a short forecast time window of around an hour.The paper is organized as follows.In Section 2, we introduce the heliospheric model of EUHFORIA and the MIT model of OpenGGCM and demonstrate the coupling between them.Section 3 lists and explains the geomagnetic indices computed using this coupling.The real events used for validating the coupling are detailed in Section 4. In Section 5, we analyze the capability of this coupling to forecast the geomagnetic indices for the chosen events 1-2 days in advance.We also define the metrics to gauge the performance of the coupling quantitatively.Finally, we summarize and discuss the results in Section 6.
2 Models and coupling methodology 2.1 EUropean Heliospheric FORecasting Information Asset (EUHFORIA) The EUHFORIA architecture comprises a coronal part, which is a 3D semi-empirical model based on the Wang-Sheeley-Arge (WSA, Arge et al., 2004) model, and a heliospheric part.The coronal part uses the photospheric magnetic field from the synoptic magnetogram maps and generates the solar wind plasma conditions using empirical relations at 0.1 au, i.e., the inner boundary of the heliospheric part (Pomoell & Poedts, 2018;Asvestari et al., 2019).The heliospheric part, a 3D time-dependent model of the inner heliosphere, numerically solves the ideal MHD equations (including gravity) employing a cell-average finite volume method in Heliocentric Earth EQuatorial (HEEQ) coordinate system.The divergence-free condition is ensured by implementing the Constrained Transport method (Evans & Hawley, 1988).CMEs are inserted into the heliospheric part as time-dependent boundary conditions, and then self-consistently advanced by the MHD equations in the heliosphere.In addition to the cone model (simplified non-magnetized spherical blob of plasma, Pomoell & Poedts, 2018;Scolini et al., 2018), magnetized CME models (Linear Force Free (LFF) spheromak model, Verbeke et al. (2019); Flux Rope in 3D (FRi3D) model, Maharana et al. (2022)) are incorporated in EUHFORIA.This makes it possible to predict the magnetic field components in the heliosphere, which mainly determine the severity of the geomagnetic storms.
In this work, we use EUHFORIA (ver 2.0) as installed on the wICE cluster of the Vlaams Supercomputer Centrum (http://www.vscentrum.be).Two nodes of this supercomputer with 72 cores per node (144 parallel processes) were utilized.The computational mesh has a radial resolution of 0.0037 au (corresponding to 0.798 R ⊙ ) for 512 cells in the radial direction between 0.1 -2 au, and an angular resolution of 2 • in the latitudinal direction between ±80 • and 2 • in the longitudinal directions extending between 0-360 • .This resolution is used at space weather modeling centers like CCMC and VSWMC.EUHFORIA typically provides global 3D MHD output (with hourly cadence, i.e., temporal resolution is one hour) and 1D time series (10-minute cadence) output at any given point in the simulation domain.In the perspective of in situ measurements, the observed solar wind plasma parameters typically have a cadence in the order of a minute, and the magnetic field components data have a cadence in the order of seconds.

Open Geospace General Circulation Model (OpenGGCM)
The OpenGGCM (Raeder, Wang, & Fuller-Rowell, 2001;Raeder, Wang, Fuller-Rowell, & Singer, 2001) is a global model of Earth's magnetosphere, ionosphere, and thermosphere.It is originally a coupling between the magnetosphere and the global ionospherethermosphere model called Coupled Thermosphere Ionosphere Model (CTIM; Fuller-Rowell & Rees, 1983;Fuller-Rowell et al., 1996;Raeder, Wang, & Fuller-Rowell, 2001).The outer magnetosphere is modeled by solving the semi-conservative form of MHD equations on a 3D stretched Cartesian grid.An explicit second-order predictor-corrector finite difference scheme is used for time stepping.The divergence-free condition is preserved by employing the Constrained Transport method.The MHD grid used in this study is based on Cramer et al. (2017).It contains 481 cells spanning −35 R E and 5000 R E in GSE x coordinate, and 180 cells covering −48 R E to 48 R E in both GSE y and z coordinates.The minimum cell sizes in the inner magnetosphere are down to 0.167 R E in x and 0.25 R E in y and z.The electrodynamic coupling of magnetosphere-ionosphere (MI) happens at the inner boundary radius of the ionosphere at 2.1 R E .This coupling is necessary to self-consistently model the closure of field-aligned currents (FACs) generated by the interaction of the solar wind and the magnetosphere, in the ionosphere (Raeder et al., 1998).MI is coupled to the default version of the CTIM, wherein, MI provides the electron precipitation parameters and magnetospheric electric field.CTIM, in turn, provides the conductances and dynamo currents to drive the ionosphere potential in the MI model.As the physics of the inner magnetosphere is not entirely represented by MHD formalism, an additional ring current model is employed for self-consistent feedback.RCM is the ring current model that solves the motion of the plasma flux tubes using the ionosphere potential, magnetic induction, and drift forces.The RCM bounds are at 10 R E on the dayside, 13 R E on the nightside, and ±12 R E on the flanks.The ring current plasma population is initialized by MHD density and pressure from OpenGGCM.After the execution of RCM, the calculated density and pressure in the RCM domain are fed back into the magnetosphere model to make it converge toward the RCM values.The MHD model calculates the precipitation and FACs normally, but the RCM influence has the effect of modifying the precipitation and FAC values indirectly.The calculated precipitation and FACs in the RCM domain are combined with those computed by the MHD magnetosphere model, to in turn drive the CTIM model.The physical quantities determined by the RCM model are used at the low latitudes and the MHD quantities at higher latitudes.A simplified schematic representation of the various couplings can be found in Fig. 2. We use the default setup of OpenGGCM (v5.0.ccmc) with the above specifications, as installed on the Marvin cluster at the University of New Hampshire.We generate 3D MHD output files with an hourly cadence and 2D files at a minute cadence.The ionosphere files are also obtained at a oneminute cadence for computing ionospheric geomagnetic indices.The advantage of using this global model is that both the magnetospheric and ionospheric effects of geomagnetic storms can be modeled.The ring current induced in the magnetospheric domain can be used to compute the Dst index.The short timescale changes in the ionosphere at different stations can be used to compute the A-indices, a proxy for substorms.

Coupling EUHFORIA and OpenGGCM
The global simulation of OpenGGCM is subjected to a time-dependent initialization using solar wind plasma and magnetic field measurements obtained from in situ solar wind satellites.In this study, we couple the results of the heliospheric part of EUHFORIA to OpenGGCM, i.e. we provide the output of EUHFORIA at Earth as boundary conditions to drive OpenGGCM.The physical quantities that are supplied are velocity (v, km/s), magnetic field ( B, nT), proton number density (n p , particles/cc), and thermal pressure (P , pPa).The magnetosphere shape and size are driven by the solar wind pressure gradient acting along the normal to the magnetosphere coming from the dynamic pressure (ρu 2 , ρ is the mass density and u is the bulk velocity), the plasma pressure (nk B T , n is total number density, k B is the Boltzmann constant, and T is the plasma temperature), and the magnetic field pressure (B 2 /2µ 0 , B is the IMF magnetic field strength and µ 0 is the magnetic permeability) in the interplanetary space.In addition to the normal pressures, tangential stresses transfer momentum and energy across the magnetopause and energize the magnetosphere (Chapter 4; Bothmer & Daglis, 2007).All the above physical parameters that are essential to drive the processes in the magnetosphere are provided by the EUHFORIA simulations.
The standard relaxation time considered for OpenGGCM simulation is around 10-12 hours before the main phase of the storm.Hence, it suffices to start feeding EUHFORIA output to OpenGGCM around 5-6 hours before the CME shock arrival which serves as a relaxation time of around 12 hours for OpenGGCM simulation until the magnetic cloud arrives at Earth.The time cadence of EUHFORIA output is 10 minutes.Given that EUHFORIA time series do not have any fluctuations in the temporal resolution below the conventional time resolution of EUHFORIA (i.e., 10 minutes), we do not increase the input resolution of OpenGGCM.The output of the coupled model is used specifically for the computation of geomagnetic indices, Dst index, and A-indices.

Geomagnetic indices -measured and modeled
Geomagnetic indices are obtained from ground magnetic field measurements and are proxies for currents in the magnetosphere and the ionosphere, which would be difficult to measure directly.They also have the advantage that they are readily available and that they cover long time periods.In fact, geomagnetic observatories have existed since the 1850s.In this section, the data sources of the geomagnetic indices (Dst and A-indices) and their processing are discussed.We also discuss the methodology to obtain the modeled geomagnetic index data by processing different outputs of the OpenGGCM simulations.

Disturbance storm time index (Dst)
Southward IMF (B z ) reconnects with Earth's northward magnetic field and deposits solar wind energy into the magnetosphere.This process can create storm conditions, wherein magnetic disturbances are induced in the magnetosphere and on the ground (Chapman, 1918).The Dst index is most sensitive to Earth's ring current, which flows in the inner magnetosphere.It is obtained from 4 stations near the equator that are nearly evenly spaced in longitude (https://isgi.unistra.fr/indicesdst.php).Dst index data has been available since 1957, but the drawback is its 3-hour cadence.An improved version of this index is the so-called SYM-H index, which is similarly obtained, but with 6 stations and a much better cadence.Data is openly available on the World Data Center (WDC), Kyoto website (https://wdc.kugi.kyoto-u.ac.jp/wdc/Sec3.html).In the following, we use SYM-H and Dst interchangeably.It can be shown that the Dst index is also a measure of the plasma energy content of the inner magnetosphere (Dessler & Parker, 1959;Turner et al., 2001).However, the Dst index is also sensitive to other magnetosphere currents, most importantly the currents on the magnetopause.Therefore, the raw Dst index is usually corrected as follows (Burton et al., 1975): where Dst * is the corrected Dst.Here, P dyn is the solar wind dynamic pressure, a is a constant representing the contribution for the ground-induced currents, b [nT/ √ nPa] is a scaling constant, and c [nT] accounts for the quiet-time contribution.These corrections and other important properties of these indices are discussed in Wanliss and Showalter (2006).Note that an enhanced ring current in the main phase of the storm corresponds to negative Dst until it reaches a minimum.After that, the recovery phase begins and the Dst values start increasing to restore back to pre-storm conditions (close to zero).
The model Dst is obtained from integrating the currents in the inner magnetosphere using Biot-Savart's law to obtain the perturbations at Earth's surface (Cramer et al., 2017).Since the Dst computed in OpenGGCM considers purely the contribution of the ring current (i.e.Dst * ), we invert Eq.1 to obtain Dst from the simulation data to compare it with the measured Dst.The empirical constants (a, b, and c) are taken from AK2 model O' Brien and McPherron (2000) and their values are a = 1, b = 7.26, c = 11, since their model resulted in the least root mean square error between the measured and empirically modeled Dst.

The auroral indices
The auroral indices or A-indices (AU, AL, and AE, Davis & Sugiura, 1966) are widely used to quantify ionospheric disturbances.They are obtained from 12 stations at auroral latitudes (around 70 degrees magnetic latitude) with more or less even spacing in longitude (See for more details, https://wdc.kugi.kyoto-u.ac.jp/wdc/pdf/AEDst version def v2.pdf) .The AU index is the upper envelope of the 12 time traces of the north-south perturbations, and, correspondingly, the AL index is the lower envelope.The AE index is the difference, AU-AL.AL is considered to be a measure of the westward electrojet, and thus a good indicator of substorms.By contrast, AU is considered an indicator of eastward currents, which occur mainly on the dayside.Unlike Dst, there is no known correlation between the A-indices and the actual currents.The traditional A-indices are no longer computed on a regular and reliable basis, mostly because data from Russian stations are lacking, which reduces the reliability of the observations used in this work.The SuperMAG magnetometer site, consolidating most worldwide magnetometer measurements, provides an alternative with their SML/SMU/SME indices, which are derived from a different set of stations (Newell & Gjerloev, 2011;DeJong et al., 2018) Convection (SMC) events (see, for example, Bergin et al. ( 2020) ), which hardly register in Dst data.From the model, we compute the A-indices using Biot-Savart integration over the ionosphere currents.Details of that procedure can be found in Raeder, Wang, and Fuller-Rowell (2001); Raeder, Wang, Fuller-Rowell, and Singer (2001).

Validation sample
In this section, we focus on two events to validate the EUHFORIA-OpenGGCM coupling.Event 1 is a single CME event where the CME was initiated on 12 July 2012 and caused a moderate geomagnetic storm.Event 2 is a more complex event resulting from the interactions between three CMEs that erupted between 4-6 September 2017, resulting in an intense geomagnetic storm at Earth.With Event 1, we aim to validate the capability of OpenGGCM to reproduce a single prolonged main-phase storm, and with Event 2, we extend its capabilities to model multi-step consecutive storms.

Event 1: 12 July 2012
This is a textbook geoeffective CME event extensively studied in the solar physics and space weather community (H.Hu et al., 2016;Gopalswamy et al., 2018;Scolini et al., 2019).The CME erupted on 12 July 2012 from the NOAA AR 11520 located at S17E06.It was associated with an intense GOES X1.4 flare.Based on the source region observations and 3D reconstruction of the CME, it was not expected to impact Earth (Webb & Nitta, 2017).However, a fast (an average projected speed of 885 km s −1 ), non-interacting halo CME arrived at Earth on 14 July 2012 and resulted in a prolonged negative B z signature (minimum B z = −18 nT), causing a moderate storm with the minimum Dst = −122 nT on 15 July 2012.The AE index reached a maximum of ∼1600 nT with fluctuations close to 1200 nT during the main phase of the storm.The planetary K (Kp) index reached above 5 during the main phase of the storm (https://www.swpc.noaa.gov/products/planetary-k-index).A detailed analysis of the geomagnetic indices threshold (including Dst) by Palacios et al. (2018) characterizes the severity of this storm towards the end of the moderate category, and close to an intense storm.Due to the prolonged nature of the negative B z , the main phase of the storm started around 10:00 UT on 15 July 2012 and maintained a Dst < −100 nT until around 09:00 UT on 16 July 2012.Comparing the in situ measurements from ACE (Chiu et al., 1998), Wind (Ogilvie & Desch, 1997), and OMNI database (https://omniweb.gsfc.nasa.gov/)during 14-17 July 2012, we found that Wind had the best data, i.e., without any gaps or anomaly.The in situ observations of speed, number density, and B z are plotted in Fig. 3, along with the demarcation of the time of arrival of shock and magnetic cloud.We choose data sources without gaps and anomaly to facilitate better input for the OpenGGCM simulations.As per the Wind ICME catalog (https://wind.nasa.gov/ICMEcatalog/), the CME shock was recorded on 14 July at ∼17:40 UT.The sheath region starting after the shock is characterized by an increased magnetic field, speed, and density.The passage of the magnetic ejecta lasts for almost 2 days between 15 July at ∼06:14 UT and 17 July at ∼03:22 UT, corresponding to low plasma beta (ratio of the plasma pressure to the magnetic pressure) and proton temperature.The B z component fluctuates between ±20 nT and stays southward for almost 20 hours, potentially responsible for loading the magnetosphere with energy and reconnecting with Earth's northward magnetic field and erosion.Raeder, Wang, and Fuller-Rowell (2001) re-emphasize that B z is most affecting in the OpenGGCM simulations.
EUHFORIA simulation was carried out for this event using the magnetized spheromak model (Verbeke et al., 2019).The initial parameters for inserting the CME into the EU-HFORIA heliospheric domain at 0.1 au are taken from Scolini et al. (2019).In this work, we simulate the CME by modifying the magnetic axis orientation, by aligning the magnetic axis of the spheromak with the observed polarity inversion line (PIL) orientation at the source region, to obtain a better fit to the in situ observations as compared to Scolini et al. (2019).The CME parameters are provided in Table 1.This modification was adopted to improve the modeling of the B z component.We have a better match with the minimum negative B z as compared to the results presented in Scolini et al. (2019).

Event 2: 4-6 September 2017
This event is associated with the interaction of the three successive CMEs that erupted between 4-6 September 2017 resulting in a two-step geomagnetic storm.It significantly impacted Earth by disrupting telecommunication during hurricane mitigation (Redmon et al., 2018), enhancing radiation dosages (Kataoka et al., 2018;Berger et al., 2018) and ground-level currents (Cohen & Mewaldt, 2018), and impacting the ionosphere, GNSS, and radio wave propagation (Yasyukevich et al., 2018).The geoeffective signatures were observed at Earth between 6-9 September 2017.The first CME (hereafter CME1) was associated with an M1.7 flare that had an average speed of about 600 km s −1 .CME1 was overtaken by the faster second CME (hereafter CME2) at ∼10 R ⊙ .CME2 was associated with an M5.5 flare that erupted from the same active region with an average speed of 1420 km s −1 .CME1 and CME2 merged in the corona, an hour after the eruption of CME2, to drive a single shock.The final full-halo CME (hereafter CME3) erupted on 6 September 2017, associated with an X9.3 flare, and propagated with a projected speed of about ∼ 1570 km s −1 .
The best in situ observations for this event obtained from the OMNI database (data gaps present in ACE and Wind data) are plotted in Fig. 6, along with the demarcation of the time of arrival of shock, sheath and magnetic cloud as reported in the Wind ICME catalog.The first shock (hereafter S1) was recorded on September 6 at around 23.00 UT followed by a 0.25 au wide sheath and then a magnetic ejecta associated with the merged CME1 and CME2 (hereafter E1) resulted in a minimum B z of around −32 nT on September 7, 23.00 UT.E1 was intercepted by the shock of CME3 (hereafter S2) which could have amplified the southward B z in E1.This enhancement led to the formation of the first and the strongest dip in the Dst index (−144 nT) on September 8 at ∼1.00 UT.There are two patches of magnetic ejecta associated with the passage of CME3 (minimum B z = −16 nT on September 8, 11:00 UT) which led to the development of the second dip in the Dst index (minimum of −124 nT) on September 8 around 13:00 UT.We perform the EUHFORIA simulation of this event using the solar wind and CME parameters (spheromak model) as used by Scolini et al. (2020).The CME parameters are listed in Table 1.The first dip in B z is captured by the EUHFORIA simulation, whereas the second dip modeled by the EUHFORIA simulation is underestimated compared to observations.

Analysis
In this section, we describe the results of the OpenGGCM simulations of the two storms described in Section 4. For each of the events, we perform two OpenGGCM simulations -one with the in situ observations recorded at L1 and the other with modeled solar wind output obtained from EUHFORIA simulation.They are named Event'n'-obs and Event'n'euh, respectively, where 'n' is the event number (either 1 or 2 for Event 1 or Event 2, respectively).We then compare the geomagnetic activity predicted by OpenGGCM, mainly the Dst and the AE index using the different sets of input.The event-specific Dst and A-indices predictions using OpenGGCM are plotted in Fig. 3 and Fig. 6.For the efficient usage of the computational resources, we choose a narrow window for running OpenGGCM, i.e., an interval corresponding to 6-12 hours before the arrival of the CME until the passage of the magnetic cloud through Earth.For assessing the prediction of the Dst index, we choose to evaluate a quite extended period before and after the Dst dip, instead of just evaluating its single minimum value.To do that, we adopt the dynamic time warping (DTW, Górecki & Luczak, 2013;Keogh & Pazzani, 2001, and references therein) technique, a method that can help us assess the overall performance of the predicted Dst time series compared to the observed data.DTW measures the similarity between two sequences that have similar patterns but differ in time.The sequences are warped in a non-linear manner to match each other.The optimal alignment of time-dependent sequences is achieved by finding an optimal path through the DTW cost matrix such that the cumulative cost is minimal compared to the other possible paths (Keogh & Pazzani, 2001).This method has been adopted in the space weather community for comparing solar wind characteristics (Samara et al., 2022) and Dst index as well (Laperre et al., 2020).For the application of DTW, we ensured that the initial and final points of the two sequences were aligned and verified the forward mapping of temporally advancing points.We further checked that there were no data gaps in our time series as required for the correct application of the DTW algorithm.We smoothed the higher cadence observed data to make them consistent with the low cadence modeled data.Windowing is applied to restrict the matching of the points within a certain time interval in order to minimize singularities, i.e., a scenario when a single data point in one sequence is mapped to a larger subsection of points in the other sequence.The window length can vary depending on the temporal spread of the features being aligned in a particular event.We first apply DTW between the measured Dst and the predicted Dst time series modeled by OpenGGCM.To do that, we calculate the DTW cost matrix based on the following equation:

Run
where D(i, j) is the cumulative DTW cost or distance, and δ(s i , q i ) = |s i − q i | corresponds to the Euclidean distance between the point s i from one time series and the point q i from the other time series.The first element of the array D(0, 0) is equal to δ(s 0 , q 0 ).The last element of this cost matrix is called the DTW score, and can be presumed as a quantification of alignment between the two time series.
To properly evaluate the performance of the predicted Dst time series, we compute the sequence similarity factor (SSF), a skill score to evaluate the performance of the predicted Dst sequences based on the ideal (observations) and a reference scenarios.In our case, the reference prediction scenario reflects a worst-case scenario in which no Dst prediction was made by the model, namely, Dst is null (i.e., no geomagnetic disturbances at all).The SSF is the ratio between the DTW score of the observed and modeled Dst time series, and the DTW score between the observed and reference (null) scenario time series.Namely, it is defined as: where O, M , and N represent the observed, modeled, and null Dst scenario (reference case) respectively.Comparing the SSF of Event'n'-obs and Event'n'-euh, the lower the SSF, the better the prediction is, compared to observations.We also estimate the time differences (∆t) and amplitude differences (∆Dst) based on the points that were aligned by DTW between the observed and modeled Dst time series (∆ refers to the difference between Observed and Modeled values).Their distributions are presented in Figs. 4, 5, 7 and 8.

Event 1
The inputs, as described in Section 3, and the results (Dst and A-indices predictions) of the OpenGGCM simulations are presented in Fig. 3.The sudden commencement (Dst starts to become negative), main (rapid decrease in Dst) and recovery (Dst increases to restore back to normal) phases of the storm are qualitatively modeled in both Event1-obs and Event1-euh.Two important aspects of comparing the Dst profile are the magnitude of the minimum Dst and the time of the beginning of the main phase.In Event1-euh the minimum Dst value is overestimated as compared to Event1-obs, and the main phase is initiated ∼2 hours earlier as compared to the observations.This corresponds to about 2 hours of early modeling of minimum B z in the EUHFORIA simulation compared to the observed profiles (panel 3 in Fig. 3).For the application of DTW, a window of 50 minutes was considered based on the time differences between features of the modeled and observed Dst for Event 1.The alignments between the predicted Dst of Event1-obs and Event1-euh with observations is shown in Fig. 4(a) and 5(a), respectively.For both Event1-obs and Event1-euh, most of the alignments showed a ∆t of +50 minutes between observed and modeled time series (see Fig. 4(b) and 5(b)).This means that, in most of the cases, the prediction of the Dst profile was earlier compared to observations.Figures 4(c) and 5(c) show the histograms of Dst amplitude difference between observed and modeled Dst time series.In both plots we notice that the ∆Dst is mostly positive, ranging between [-20, 85] nT and [-40, 105] nT, for Event1-obs and Event1-euh, respectively.This indicates that the predicted Dst, for most of the alignments, was more negative compared to observations.The SSF for Event1-obs is slightly lower than Event1-euh (cf.Table 3 for more details).This implies that the OpenGGCM predictions made with the measured solar wind properties can reproduce the observed Dst slightly better, although not significantly better, as compared to the predictions using EUHFORIA simulated data.
The AE index is a proxy of the response of the ionosphere to the substorms which are quite stochastic and have high temporal variability.The predicted AU, AL, and AE indices match the order of magnitude and show rapid variations during the period of the storm.Although there is a reasonable agreement in the order of magnitude of the AE index, it is not straightforward to align the peaks of the predicted AE indices to their measured counterparts to assess the performance of OpenGGCM.The extreme Dst, AU, AL, and AE indices measured and predicted by Event1-obs and Event1-euh, are provided in Table 2.

Event 2
The comparison of the predicted Dst and AE indices with the measurements is presented in Fig. 6.For this case, the Dst prediction of Event2-euh is qualitatively closer to the measurements as compared to Event2-obs.The two-step feature of the storm is well reproduced by Event2-euh.Event2-obs underestimates the first dip of the Dst.It is a possibility that the first measured Dst dip (8 September at 01:00) is predicted too early in Event2-obs on 7 September around 18:00.Otherwise, the first dip in Event2-obs Dst could be the overestimated Dst due to the slight negative B z occurring on 7 September around ∼07:00 -10:00.The DTW window of 120 minutes was applied based on the time differences between features of the modeled and measured Dst for Event2.The DTW alignments of Event2-obs and Event2-euh with the measured Dst values are shown in Fig. 7 to 120 minutes, meaning that the Dst patterns were mostly predicted earlier than observed.The opposite is recorded for Event2-euh.Namely,Fig. 8(b) shows that for the majority of alignments, ∆t is negative and equal to −110 to −115 minutes (most Dst patterns were predicted later than observed).Figures 7(c) and 8(c) present the histograms of the Dst amplitude difference between measured and predicted Dst, for Event2-obs and Event2-euh, respectively.Same as in Event 1, we notice that in both plots ∆Dst is positive for most of the alignments, meaning that the predicted Dst index was mostly overestimated (i.e., more negative) compared to observations.The SSF for Event2-euh is clearly lower than Event2-obs indicating better Dst prediction by OpenGGCM using EUHFORIA input (cf.Table 3 for more details).
The AE index has two distinct regions of enhancements on 8 September 2017, first around 00:00 UT and the other around 16:00 UT.Qualitatively, the first peak is modeled by both Event2-obs and Event2-euh, whereas the second peak is better modeled by Event2obs.Due to high temporal variability, we could not perform a sequence alignment with DTW on the AE time data for this event for the same reasons as mentioned for Event 1.

Sources of error
The first source of error enters through the process of constraining the CME and solar wind parameters for the EUHFORIA simulations.3D reconstruction of CMEs is subject to human bias (Verbeke et al., 2022).Methods of extracting magnetic field parameters for eruptions are dependent on its source region signatures which can be missing in the cases of "problem storms" (Dodson & Hedeman, 1964).The errors in CME arrival time and magnetic field predictions are also associated with the ambient solar wind through which the CME propagates.The solar wind maps are limited by the availability and the quality of the synoptic magnetogram, and the semi-empirical WSA model which is one of the simpler coronal models.All these factors determine the accuracy of the predicted B z profile.The second source of error is associated with the sensitivity of OpenGGCM to the first B z value provided as initial condition.The occurrence of rapidly and unnaturally varying features in the modeled Dst profile midway through the main phase of the storms and continuing all the way up to the recovery phase suggested the commencement of the accumulation of the instabilities in the simulations.So, we made a conscious strategic effort to choose the starting point such that B z was low but non-zero to avoid abrupt gradients in the simulation domain.Hence, the accuracy of Dst predictions of OpenGGCM is strongly correlated with the B z predictions of EUHFORIA.The final errors could be associated with the coupled ring current model with OpenGGCM which makes the Dst calculations.

Summary and Discussion
An overview of the methodologies adopted in this study is as follows: • We successfully demonstrate the coupling between the solar wind and CME evolution model, EUHFORIA, and the coupled magnetosphere-ionosphere-thermosphere model, OpenGGCM.
• The geomagnetic indices are derived by post-processing the output of this newly coupled model.The Dst index was computed using the ring current model, RCM, coupled to OpenGGCM.The A-indices were quantified by averaging the magnetic field perturbations in the ionosphere.• In addition, we have presented the validation of this coupling with the widely studied moderate-to-intense geomagnetic storms of July 2012 (Event 1) and September 2017 (Event 2).Two OpenGGCM simulations were performed for each event -one employed the in situ observations (Event'n'-obs), and the other used the EUHFORIA output (Event'n'-euh) for initializing OpenGGCM.• We then assessed the predicted Dst profiles of Event'n'-obs and Event'n'-euh using the DTW technique to evaluate the whole range of Dst predictions, not just the minimum Dst point.
In the case of Event 1 (CME on 12 July 2012) Dst prediction based on the in situ observations (Event1-obs) slightly outperformed the prediction made by EUHFORIA input (Event1-euh).The predicted AE indices, in both Event1-obs and Event1-euh, although qualitatively matched the order of magnitude of the observations during the period of the storm, a one-to-one correspondence between the extreme points was not attained.For Event 2 (CME on 4-6 September 2017), the Dst predictions by OpenGGCM with EUH-FORIA input outperformed the one with in situ observations.In this two-step storm, the enhancements in the predicted AE indices were clustered around the two peaks as in the observations.AE indices were captured better for Event 2. We conclude that EUHFORIAdriven OpenGGCM simulations, although overestimating the Dst, are still in reasonable agreement in order of magnitude with the in situ driven OpenGGCM simulations.The Aindices predicted by the EUHFORIA-driven OpenGGCM simulations match qualitatively the overall periods and the magnitudes of the enhancements.
The two event studies were merely used to showcase the validation of the coupling.It is not possible to comment on why OpenGGCM predictions were better using EUHFORIA inputs over the actual L1 observations, in Event 2 but not in Event 1.Further studies need to be undertaken to understand the effect of different solar wind magnetic fields and flow speed profiles on the variability of the OpenGGCM predictions.Parametric studies with synthetic EUHFORIA simulation inputs can improve our understanding of the drivers of storms and substorms.EUHFORIA simulations can be performed as fast as 25 times the physical time of the process.With sufficient resources OpenGGCM can be run easily as fast as 10 times real time, making long-range predictions (typically several days) possible with the coupled models.This study is one of the efforts to couple MHD models of the heliosphere with CMEs and global magnetosphere to improve space weather forecasting, which lays the groundwork for future improvements.The simulations discussed in this work show a potential to obtain forecasts 1-2 days in advance, with a similar accuracy as the current nowcasts.

Figure 1 .
Figure 1.Schematic showing the different models within EUHFORIA and the outputs they produce.On the left, the corona is modeled using the semi-empirical Wang-Sheeley-Arge (WSA) model, driven by the synoptic magnetogram maps (e.g., top left figure).WSA model employs the potential field source surface (PFSS) model in the low corona and the Schatten current sheet (SCS) in the upper corona to extrapolate the magnetic field lines up to 0.1 au and employs empirical relations to compute MHD parameters at 0.1 au (e.g., bottom left figure).The output of the coronal model is provided as a boundary condition to the 3D time-dependent ideal MHD model of the heliosphere, shown on the right.Figure on the right is an example showing the radial speed (vr) profile in the EUHFORIA domain depicting a propagating CME on a relaxed solar wind background (equatorial and meridional planes containing Earth).

F10Figure 2 .
Figure 2. Flowchart demonstrating the coupling of the MHD magnetosphere model of OpenG-GCM with CTIM and RCM, adapted from Cramer et al. (2017).B, n, and T are the magnetic field, density, and temperature of the plasma.All other abbreviations and the meanings of the symbols are provided on the left side of the figure.

Figure 3 .Figure 4 .Figure 5 .
Figure 3. Characteristics and predicted geomagnetic indices of Event 1. Panels 1-3 show the plasma parameters -speed (v), proton number density (np), and the magnetic field parameter -z-component of magnetic field (Bz) as obtained from the Wind spacecraft in situ observations (in black) and the EUHFORIA simulation of the event based on modified Scolini et al. (2019) (in blue), respectively.The horizontal blue line in Panel 3 corresponds to Bz = 0. Panels 4-7 show the geomagnetic indices -Dst index, AU index, AL index, and AE index as measured in Earth's magnetosphere and ionosphere (in red), and as obtained from OpenGGCM simulations using input from the Wind (in black) and EUHFORIA simulation (in blue).The magenta and green vertical solid lines depict the arrival of the CME shock and the beginning of the magnetic cloud passage at Earth, respectively.

Figure 6 .
Figure 6.Characteristics and predicted geomagnetic indices of Event 2. Panels 1-3 show the plasma parameters -v, np, and Bz as obtained from the in situ observations in the OMNI database (in black) and the EUHFORIA simulation of the event based on Scolini et al. (2020) (in blue), respectively.Panels 4-7 show the geomagnetic indices -Dst index, AU index, AL index, and AE index as measured in Earth's magnetosphere and ionosphere (in red), and as obtained from OpenG-GCM simulations using input from the OMNI database (in black) and EUHFORIA simulation (in blue).The magenta solid and dashed lines depict the arrival of two shocks (S1 and S2) associated with this event.The two green solid lines depict the boundary of the passage of the magnetic ejecta E1 at Earth and the dashed lines correspond to the boundary of E2 at Earth.These interplanetary CME features are detected 40 minutes later at Earth (as per OMNI database) relative to L1 (as per WIND ICME catalog).

Figure 7 .Figure 8 .
Figure 7. DTW analysis of the Dst index.(a) DTW alignment between the observed (blue) and predicted/modeled (red) time series, (b) Histogram of time differences between the aligned points, and (c) Histogram of Dst differences between the aligned points for Event2-obs.Each time element in (a) corresponds to 5-minutes.

Table 1 .
. Unlike Dst, the A-indices also respond to much weaker forcing, i.e., substorms, pseudo-breakups, and Steady Magnetosphere Parameters of the single CME erupting on July 12, 2012 (Event 1) and of the three CMEs erupting between September 4-6, 2017 (Event 2) used for initializing CMEs at the EUHFO-RIA heliospheric inner boundary at 0.1 au.

Table 2 .
Specifications of the EUHFORIA runs of Event 1 and Event 2 used for validating the EUHFORIA and OpenGGCM coupling.In column 1, Run ID defines an identifier of each simulation, and the input column specifies the source that drives OpenGGCM.The extrema of

Table 3 .
EventDTW score (O, M ) DTW score (O, N ) SSF DTW scores for each OpenGGCM simulation and their sequence similarity factor (SSF)