A new model describing Forbush Decreases at Mars: combining the heliospheric modulation and the atmospheric influence

Forbush decreases are depressions in the galactic cosmic rays (GCRs) that are caused primarily by modulations of interplanetary coronal mass ejections (ICMEs) but also occasionally by stream/corotating interaction regions (SIRs/CIRs). Forbush decreases have been studied extensively using neutron monitors at Earth; recently, for the first time, they have been measured on the surface of another planet, Mars, by the Radiation Assessment Detector (RAD) on board the Mars Science Laboratory's (MSL) rover Curiosity. The modulation of GCR particles by heliospheric transients in space is energy‐dependent; afterwards, these particles interact with the Martian atmosphere, the interaction process depending on particle type and energy. In order to use ground‐measured Forbush decreases to study the space weather environment near Mars, it is important to understand and quantify the energy‐dependent modulation of the GCR particles by not only the pass‐by heliospheric disturbances but also by the Martian atmosphere. Accordingly, this study presents a model that quantifies — both at the Martian surface and in the interplanetary space near Mars — the amplitudes of Forbush decreases at Mars during the pass‐by of an ICME/SIR by combining the heliospheric modulation of GCRs with the atmospheric modification of such modulated GCR spectra. The modeled results are in good agreement with measurements of Forbush decreases caused by ICMEs/SIRs based on data collected by MSL on the surface of Mars and by the Mars Atmosphere and Volatile EvolutioN (MAVEN) spacecraft in orbit. Our model and these findings support the validity of both the Forbush decrease description and Martian atmospheric transport models.


Introduction and Motivation
Galactic Cosmic Rays (GCRs) are energetic charged particles, comprised of 2% electrons and 98% atomic nuclei, with the latter made up (in number density) of 87% protons, 12% helium, and about 1% heavier nuclei (Z ≥ 3) (Simpson, 1983). GCRs are omnipresent in the heliosphere with a flux that is nearly isotropic. However, their intensity varies as a result of modulation in the heliosphere due to levels of solar magnetic activity, which evolve both in the long term, following the 11-year solar cycle, and in the short term, typically during solar eruptions and recurring coronal holes.
Forbush decreases (FDs) are temporary reductions in the intensity of the GCRs, first discovered by Forbush (1937) using groundbased measurements at Earth. The rapid and temporary depression of an FD is typically followed by a comparatively slower recovery phase, generally lasting for a few days. The causes of FDs can be interplanetary counterparts of coronal mass ejections (ICMEs) and/or the stream interaction regions (SIRs) that result from interaction of the fast and slow solar wind streams (Cane, 2000;Richardson, 2004). FDs with sporadic character are generally caused by ICMEs, which are often associated with a shock front followed by ejecta, both of which modulate the intensity of GCRs in the interplanetary space. On the other hand, recurrent FDs are related to SIRs, which often occur periodically as the solar coronal holes generating high speed streams that may exist through several solar rotational periods. The frequency with which sporadic or recurrent FDs occur depends on the solar activity that evolves over the solar cycle (Melkumyan et al., 2019).

GCRs and FDs Measured at Earth and Mars
Cosmic particles arriving at the top of the atmosphere of Earth or Mars experience ionization energy loss and also generate secondary particles via fragmentation and spallation processes as they interact with the atmosphere, making the particle spectra and energy at the surface of the planet completely different from those in space. Due to Earth's thick atmosphere, the cutoff energy at the planet's magnetic poles is approximately 450 MeV (Clem and Dorman, 2000). (Strictly speaking this cutoff energy is not a constant, as it depends on the atmospheric depth and solar cycle; it is somewhat lower at Earth's south pole than at the north pole because of the south's greater atmospheric depth.) At locations where the geomagnetic field is not negligible, there is an additional shielding of the magnetosphere, described by the local magnetic cutoff rigidity at different latitudes of Earth; for primary protons, the cutoff can reach up to a few GeV. FDs have been measured routinely at the surface of Earth for a few decades, e.g. using neutron monitors, muon telescopes (Arunbabu et al., 2015), and even water-Cherenkov detectors (Dasso et al., 2012). Studies of terrestrial FDs based on different geomagnetic cutoff rigidities have shown that there is an energy dependence of the FD amplitude (Lingri et al., 2016) and its recovery time (Usoskin et al., 2008) due to the energy-dependent modulation of the GCRs by heliospheric disturbances. Simply speaking, the modulation of GCRs is stronger for lower-energy particles than for higher-energy ones. Therefore the same FD has a larger amplitude when measured by stations at locations of lower magnetic cutoff rigidity.
Since August 2012, GCR fluxes and corresponding FDs therein have been measured on the surface of Mars (Guo J et al., 2018b), at Gale Crater, by the Radiation Assessment Detector (RAD, Hassler et al., 2012), on board the Mars Science Laboratory's (MSL) rover Curiosity. Different from Earth, Mars lacks a global magnetic field; thus there is no magnetic cutoff energy for particles measured at Mars. The Mars atmosphere is composed of 95% CO 2 ; the surface pressure measured at Gale Crater averages approximately 840 Pascals (Haberle et al., 2014), which is less than 1% of that at the sea level of Earth. Nevertheless, this thin atmosphere is able to slow down primary GCR particles, especially the lower-energy ones, and can cause high-charged particles to fragment in the atmosphere, resulting in an altered surface particle environment. Modeling results show that the Martian atmosphere shields away GCR protons below an energy of about 140 to 190 MeV at Gale Crater (Guo J et al., 2018a), depending on the atmospheric depth, which varies seasonally by as much as 25%.

Motivation of the Current Work
In order to study the propagation of ICMEs/SIRs in the interplanetary space where in situ measurements are scarce, it is important to use, as much as possible, any information available from different heliospheric and planetary missions. As the FD onset time matches very well with the arrival of the corresponding solar wind disturbance structures, FDs have been used successfully as a proxy to determine ICME or SIR arrival times, especially when no plasma or magnetic field measurements are available (e.g., Möstl et al., 2015;Witasse et al., 2017;Wang YM et al., 2018;Guo J., 2018c;Winslow et al., 2018;Freiherr von Forstner et al., 2018;Papaioannou et al., 2019;Freiherr von Forstner et al., 2019).
On the other hand, the magnetic field in ICMEs generally weakens due to its expansion during the prop-agation (e.g., Démoulin and Dasso, 2009) and the shock also evolves as it propagates into the interplanetary space. Winslow et al. (2018) suggest that FDs are steeper and deeper closer to the Sun, based on depressions in the GCR flux caused by an ICME measured at various locations by three planetary missions: by MErcury Surface, Space ENvironment, GEochemistry, and Ranging (MESSENGER) at Mercury; by Lunar Reconnaissance Orbiter (LRO) at the Moon; and by MSL at Mars.
However, such studies remain very rare despite the fact that almost every heliospheric and planetary mission carries a high-energy particle detector or a radiation monitor measuring the background GCR flux. One major barrier is the different energy ranges measured by different instruments chosen for different scientific purposes in different environments. For instance, RAD measures GCR dose rate with an atmospheric cutoff energy of 140 to 190 MeV (see Section 1.1); but even above this energy, the contribution to the surface radiation is not constant over energy. The largest contribution at the surface comes from the secondary particles produced in the atmosphere by primary GCR protons in the ~1 to 3 GeV energy range (Guo J et al., 2019b). On the other hand, the MESSENGER Neutron Spectrometer is designed for the detection of neutrons generated by GCRs through nuclear spallation reactions with Mercury's surface (Feldman et al., 2004); but the actual energy range of the primary GCRs from such detection is again a complex result of the incoming particles interacting with the planetary environment. While with the same ICME strength, the modulation of GCR particles and the corresponding FD amplitude is already energy-dependent (see Section 1.1), it is therefore very challenging to derive the evolution of the ICME property during its propagation based on FD amplitudes measured by different instruments recording in different GCR energy ranges.
In order to employ the abundant particle measurements in space to describe FD and ICME properties, we need to address two issues: 1) proper evaluation of the energy range of the primary GCR particles detected by each instrument in use and 2) better description of the energy-dependent modulation by ICMEs of the GCRs and corresponding FDs. Part 1) requires detailed studies of the GCR interaction with the environment in which measurements were made, and a thorough understanding of each detector's response function. In this study, we show how this can be achieved for the case of RAD, employing a model describing the particle interaction with the Martian atmosphere. Part 2) needs reliable models/measurements of the modulation of FDs, which has been a long-standing problem in the field of heliospheric studies.
There have been successful analytic models describing the formation and profile of FDs. For instance, the diffusive barrier model (Wibberenz et al., 1998) describes the diffusion process of GCR particles into the ICME sheath region, which is the part between the leading shock and the proceeding magnetic ejecta that is characterized by enhanced solar wind speed and turbulence. Luo X et al. (2017) built a full 3D time-dependent GCR transport model that includes diffusion, drift, solar wind convection, and adiabatic cooling. The reduced diffusion mechanism is utilized to model a 3D propagating barrier, i.e., the sheath region, for simulating FD profiles. The ForbMod model  takes into account the diffusion of GCR particles into a cylindrical magnetic flux rope, which is typical of an ICME's magnetic ejecta. There is a concurrent process of flux rope expansion and increased diffusion of GCRs into the rope as the ICME propagates further; this process drives the evolution of the FD profiles. However, both models depend on the energy-dependent diffusion coefficient, which is unknown; therefore the modeled FD is often scaled to an arbitrary amplitude.
Recently, Usoskin et al. (2015) has used a force-field parametrization approximation to describe the energy-dependent modulation of GCRs and the corresponding FDs (more details of the force-field approximation appear in Section 2.2). In particular they used a force-field model to fit the hydrogen and helium spectra obtained by the Payload for Antimatter Matter Exploration and Light-nuclei Astrophysics (PAMELA, Adriani et al., 2011) in December 2006 before and during a significant FD caused by a major ICME. They found good agreement between the fitted force-field parameter and the one obtained from the neutron monitor network, and concluded that the force-field solution is a reasonable mathematical description of the GCR energy spectra, in the energy range from a few hundred MeV up to 100 GeV, in different phases of the FD.
However, the force-field assumptions based on steadiness of the solar wind are likely violated on the short time scales and disturbed heliospheric conditions of an ICME event. Besides, the simple parametrization to describe the heliospheric transport and modulation of the GCR spectrum is questionable, and the modulation potential derived from the force-field approximation has limited physical meaning. Recently, the force-field formalism has been improved by Corti et al. (2016) and Gieseler et al. (2017) in order to take into account the rigidity dependence of the diffusion tensor. These models updated the classical one-parameter force field solution, with two different Φ values at low and high rigidities and a transition function in between. The new approach improves the fitting at both the low and the high energy parts of the spectrum.
However, the classical one-parameter force-field solution does allow investigation of long-and short-term variations in limited energy ranges, i.e., from 1 to 50 GV, as suggested by recent Alpha Magnet Spectrometer (AMS 02) measurements (Aguilar et al., 2018). This is the energy range of primary particles that contribute mostly to the Terrestrial or Martian surface particle environment (Clem and Dorman, 2000;Guo J et al., 2019b) where GCRs and FDs are measured. Therefore, we adopt a similar parametrization approximation to describe the energy-dependent modulation of GCRs.
In the current paper, we develop a model that combines the heliospheric modulation of GCRs (i.e., Part 2 above) and the Martian atmospheric modification of such modulated GCR spectra (i.e., Part 1 above) to quantify the amplitudes of the Forbush decreases at Mars during the pass-by of ICMEs -both in the interplanetary space near Mars and on the planet's surface. We also review a statistical study of Forbush decreases measured on the surface of Mars by MSL and by the Mars Atmosphere and Volatile EvolutioN (MAVEN) spacecraft in orbit, and compare these empirical observations to our modeled results. The comparison shows good quantitative agreement, suggesting that both the Forbush decrease description and the Martian atmospheric transport models are valid.

The GCR and FD Measurement at Mars
2.1.1 FDs observed on the Martian surface On the surface of Mars, RAD measures a mix of primary GCRs and secondary particles generated in the atmosphere including both charged  and neutral Guo J et al., 2017b) particles. RAD also measures dose in two detectors: the silicon detector B and the plastic scintillator E. Radiation dose is defined as the energy deposited by radiation per unit mass with a unit of J/kg (or Gy). Dose rate is a key quantity used to evaluate the energetic particle environment. Detector B measures a lower dose rate because the ionization potential of silicon is larger than that of plastic, as shown in Figure 1a. Detector E measures dose over a 4π geometric angle with a much larger geometric factor than detector B, resulting in a better signal-to-noise ratio. Therefore we use the dose rate measured in the plastic scintillator E as a good proxy for GCR fluxes and for detecting FDs. Figure 1a shows a typical FD measured by RAD, which lasted for about 14 days. The originally surface dose rate data show a diurnal cycle that has been found to be correlated with the Martian daily thermal and pressure cycle at Gale crater (Rafkin et al., 2014;Guo J et al., 2017a). We have filtered out this daily oscillation using a notch filter tuned to remove all the multiples of 1 sol harmonic, retainingthe lower and higher frequencies of the signals in the data (Guo J et al., 2018b). Both the daily average and solfiltered (1 sol = 1 Mars day ≈ 1.03 Earth day) dose rates are also over-plotted on the original data. RAD-seen FDs are then identified based on the sol-filtered data. Figure 1b shows the histogram of FD magnitudes at Mars' surface seen at MSL for all FDs from October 2014 to September 2016 (Guo J et al., 2018b). The FD magnitude is in units of drop percentage (%) and the bins in the x-axis are in logarithmic scale. The yaxis shows the normalized histogram (counts per %) and the red curve is the result of a power-law fit. The legends show 1) the fitted power-law index and its uncertaintiy as −2.08 ± 0.32; 2) the correlation coefficient between the FD magnitude and its distribution, which is −0.77; and 3) the mean value and standard deviation of the FD distribution, which are 4.35% and 2.56%, respectively. The maximum amplitude of the FDs detected over these two years is 18.10% ± 1.38% on 2014 October 17, which was caused by an intense ICME event that passed STEREO-A, Mars, comet 67P, Saturn, and New Horizons en route to Pluto, ranging up to 31.6 AU (Witasse et al., 2017).

FDs observed in Mars' orbit
On top of the Martian atmosphere, the MAVEN spacecraft has been continuously monitoring the local space weather conditions around Mars since October 2014, following its arrival at the red planet. The solar energetic particle (SEP) instrument on board was designed to measure energetic solar particles and pickup ions at Mars within an energy range of 20 to 6000 keV for ions and 20 to 1000 keV for electrons, as shown in Fig. 6 in Larson et al. (2015). These are not really suitable energy ranges for identifying FDs, which are depressions in the highly energetic (100 MeV/nuc-GeV/nuc) GCR fluxes. However, in a not-so-straightforward manner, SEP data can be selected to measure particles penetrating through the whole detector set and ions with energies larger than 100 MeV/nuc, with a contribution by electrons with energies ≥ 600 keV during solar events. Since GCRs are composed of mainly protons, it is reasonable to assume that FTO penetrating channels measure mainly GCR protons during solar quiet times. More detailed explanations of how the MAVEN/SEP integrated count rate of protons above 100 MeV can be employed to study FD properties on top of the Martian atmosphere can be found in Guo J et al. (2018b). fewer FDs than those seen on the surface (panel b) because MAVEN/SEP detects more electrons in the SEP events which (partly) coincide with the ICME passage making it impossible to identify FDs in the data. The distribution has a correlation coefficient of −0.79 and a mean value of 6.37% and the fitted powerlaw index is −1.84 ± 0.55. The maximum amplitude of the FDs detected over this period is 26.93% ± 1.94%, on 2014 October 17 as also seen by MSL/RAD. Statistically speaking and as expected, the FDs detected in orbit have larger amplitudes than those seen on the surface and the distribution has a flatter power-law shape.
The comparison of the FD amplitude above and below the atmosphere is better shown in Figure 1d where each point corresponds to one event detected simultaneously by MAVEN and MSL. An error bar of 1.38% for FDs at MSL and 1.94% for FDs at MAVEN has been added to each data point, accounting for the uncertainties resulting from statistical errors and manual determinations of FD onsets and the nadir points (Guo J et al., 2018b). The fitted linear slope representing the averaged ratio between each FD's amplitude measured in orbit and measured at the planetary surface is about 1.33 ± 0.14. In other words, an FD event seen on the   surface of Mars by MSL with ~10% drop magnitude corresponds to ~13 % drop magnitude of GCRs with energies ≥ 100 MeV outside Mars' atmosphere. This difference is related to the energy-dependent modulation of the GCR particles by both the ICME/SIR and the Martian atmosphere, and has been explained, qualitatively, by Guo J et al. (2018b). We do not have sufficient plasma and magnetic field data to identify each FD as being ICME or SIRcaused. But the FD magnitude and the ratio fitted by observations ( Figure 1) do not show two distinct populations; moreover, since plasma in SIRs is similar to that in the turbulent sheath region of the ICMEs, we assume that the energy-dependent modulations during FDs are similar during ICMEs and SIRs. Besides, based on Ulysses measurement, Paizis et al. (1999) have shown that the amplitude of FDs decreases with energy which is in qualitative agreement with what is observed in the case of ICMEs (Lingri et al., 2016, e.g.). In this study, we make a closer examination of the concurrent and energy-dependent modulation of the GCRs by both heliospheric disturbances and the planetary atmosphere. We analyze this phenomenon quantitatively using a combination of the heliospheric modulation model (Section 2.2) and the Martian particle transport code (Section 2.3); results are presented in Section 3.

The GCR Modulation Model
2.2.1 GCR proton spectra A commonly used description of heliospheric modulation is the modulation parameter Φ (Moraal, 2013, and references therein), which is directly associated with solar activity and is proportional to the momen-tum/charge required for a particle to penetrate the heliosphere.
Considering the full equation of GCR transport in the heliosphere, the GCR distribution is determined by the outward convection of solar wind, the gradient and curvature drifts in the global heliospheric magnetic field, the diffusion through the irregular heliospheric magnetic field, and the adiabatic energy loss due to the divergence of the expanding solar wind (Parker, 1965). Assuming spherical symmetry and if only the convection and diffusion processes are considered (Moraal, 2013), the GCR distribution can be simplified as a force-field solution depending on a single modulation potential Φ, or two Φ values using the rigidity-dependent diffusion tensor (Corti et al., 2016;Gieseler et al., 2017).
Given its assumptions, the force-field approximation for describing the GCR spectra under different solar modulation condition has its own caveat. More physical than the force-field approximation, the Badhwar-O'Neill model (BON, O'Neill, 2010) describes the GCR spectra via the spherically symmetric Fokker-Planck equation that accounts for particle propagation in the heliosphere due to diffusion, convection, and adiabatic deceleration. The boundary condition is the constant energy spectrum, i.e., the Local Interstellar Spectrum (LIS), for each GCR element at the outer edge of the heliosphere (~100 AU). The current BON model ignores the drift motion of charged particles which depends on the solar cycle polarity. Data analysis by Munini et al. (2018) and the FD model of Luo X et al. (2017) show that the drift of GCRs in the heliosphere may result in different FD amplitudes and recovery times during different solar polarity cycles. However, the BON model depends on a single modulation parameter Φ and provides specie-and energy-dependent particle fluxes, making it easy and efficient to implement. An example of the proton and helium ion particle spectra (together contributing about 99% to the the GCR ion flux) is shown in Figure 2 under different solar minimum to solar maximum conditions: big (or small) values of Φ represent strong (or weak) solar modulation during solar maximum (or minimum).
We note that the values of the solar modulation parameter Φ represent strong (or weak) solar modulation during solar maximum (or minimum). used in the BON model -and in the force-field model of Usoskin et al. (2015) -may not be exactly the same during the same time, due to differences in physical assumptions. Besides, Φ= can be determined (in principle) by a number of methods, such as sunspot number, radio and X-ray flux, neutron monit- or rates, and spacecraft GCR measurements. The Φ value in Usoskin et al. (2015) is derived using data from the world network of neutron monitors validated against balloon-and space-measurements of the cosmic ray spectrum (Usoskin et al., 2011). Alternatively, the BON model uses the International Sunspot Number to determine the level of solar modulation. Since the heliosphere's response has a lag of 8-14 months to account for changes in solar activity, Φ determined from sunspot numbers tends to precede the GCR modulation. To compensate for this time lag, the BON 2010 model uses actual spacecraft data (IMP-8 and ACE) to further calibrate the sunspot number.
Since the measurement-calibrated analytic BON model is more physical and provides the full energy spectra of GCRs of different species, which can not be easily obtained from GCR measurements (with their limited energy range, limited particle species, and constrained time periods), we employ the BON model to approximate the GCR spectra driven by variation of solar and heliospheric activities. In contrast to the Ansatz used by Usoskin et al. (2015) who obtained the modulation potential during the FD, we carry out the modelling approach in an opposite direction: we assume that the modulation parameter Φ changes during the FD and thus calculate the amplitude of the FD using the GCR spectra from the BON model under different Φ values. We apply this approach to calculate the FD amplitudes on top of the Martian atmosphere and on the surface of Mars in order to compare our model results with observations, as shown in Section 2.1.

The Martian Particle Transport Model
In order to calculate the GCR spectra and corresponding FD amplitude on the surface of Mars, we need to simulate particles going through the Martian atmosphere and regolith environment.
Hereby we employ the PLANETOCOSMICS model, which is a welldeveloped tool based on GEANT4, a Monte Carlo approach to simulating the interactions of particles as they traverse matter (Agostinelli et al., 2003). Different settings and features, e.g. the composition and depth of the atmosphere and the soil, can be used in the simulations. For the Martian environment, we use the Mars Climate Database (MCD), which has been created using different Martian atmospheric circulation models that are further compared and modified by observations based on data from past and current Mars missions (Lewis et al., 1999).
In the current model setup, we use the atmosphere profile from MCD above the ground at the location of Gale Crater, which is the landing site of the Curiosity rover. More detailed descriptions of features of the MCD implemented in the model used here can be found in Guo JN et al. (2018a, 2019a and Appel et al. (2018). Particularly relevant to the current study is that the GEANT4/PLAN-ETOCOSMICS model has been shown to be able to predict reliably the surface charged particle spectra and dose rates of the GCR-induced surface radiation (Matthiä et al., 2016) as actually measured by RAD.
In particular, we have developed some ready functions which can be used to convert quickly any given interplanetary GCR proton/helium ion spectra to the radiation dose on the surface of Mars (Guo J et al., 2019b). In the current work, we use these functions to calculate the surface dose rate and analyze how it changes during different solar modulation conditions, with GCR spectra obtained from the model described in Section 2.2.

Results and Discussion
As described in the last section, we use the BON model to approximate the GCR spectrum and its variation during different solar activities, such as the pass-by of ICMEs or SIRs. For instance, before the arrival of an event, the solar modulation condition is at about Φ = 600 MV (solar minimum); during the event the GCR flux is depressed, due to enhanced solar wind turbulence and the magnetic obstacle effect. If the GCR spectra during the solar eruption can be represented by the BON spectra of Φ = 700 MV, we can approximate the modulation of the ICME on the GCR spectra using an enhanced δ Φ of 100 MV. The total amplitude of the FD can be modeled as the drop ratio of the GCR flux/dose rate, integrated over measured energy ranges (Section 2.1), before the onset and at the nadir point of the FD.
At the top of the atmosphere or in deep space, we integrate the flux of protons and helium ions, shown in Figure 2, with energies above the cutoff energy E c before and during the FD to be C bkg and C min , respectively. We then calculate the amplitude of the FD to be: where F H and F He are the GCR proton and helium ion spectra, which are a function of particle energy and modulation parameter Φ; E c is 100 MeV/nuc as measured in the MAVEN/SEP penetrating counter during a period of time with no SEPs. We ignore the GCR flux contributed by heavier ions, as they contribute only about 1% of the total flux. (We do note that the actual measurement of the particle flux may have an instrumental response that is energy-dependent; this requires further careful investigation and inclusion of the yield function of each instrument. Such a defined FD amplitude is a function of Φ and δ Φ -A(Φ, δ Φ ) -has been calculated over the range of Φ between solar minimum (400 MV in the BON model) and solar maximum (1700 MV) conditions, and of δ Φ between 10 MV (for weak solar disturbance) and 300 MV (for extreme solar eruptions such as the event studied by Usoskin et al. (2015)), as shown in the top panel of Figure 3. The magenta area covers the range of FD amplitudes, up to about 27%, measured by MAVEN in orbit during the studied period (Section 2.1). As the solar modulation condition was mostly weak during the period of the studied data, with Φ values no larger than 700 MV (http://cosmicrays.oulu.fi/phi/), we also set a constraint of Φ < 700 MV on the highlighted area. Note that this is an approximate constraint of the observation, since the Φ values are not necessarily the same between the Neutron-monitor-based force-free solution and the BON model, as explained in Section 2.2.
Similarly, we calculate the amplitude of the FDs on the surface of Mars -B(Φ, δ Φ ). We use the Martian atmospheric and regolith condition at Gale crater (the location of the MSL rover), with a pressure of about 830 Pascals, close to the observed annual average pressure. Instead of using flux rate, we calculate the relative decrease of the surface dose rate induced by GCR spectra under different modulation conditions, described by the following equation: where D is the surface dose rate measured by MSL/RAD and is a function of Φ. Note that the dose measurement includes the total energy deposited by surface particles of different types (both primaries and secondaries of different species) and with different energies, which can be calculated via folding the primary GCR proton and helium ion spectra with the atmospheric dose rate functions (which are a function of the incoming particle energy and type) obtained by Guo J et al. (2019a); a mathematical description of D and numerical calculations of the functions can also be found in this reference. The middle panel of Figure 3 shows the calculated surface FD amplitude B(Φ, δ Φ ) under different Φ and δ Φ values. The magenta area covers the range of FD amplitudes, up to about 19%, measured by MSL on the ground with a constraint of Φ < 700 MV during the studied period. The points covered in this plot are almost the same as in the upper panel, supporting the consistency between the measured and modeled results both above and below the atmosphere.
The bottom panel of Figure 3 shows the ratio of A(Φ, δ Φ )/B(Φ, δ Φ ) for each case calculated. The ratios fall within the range of about 1.3 to 1.45, with the magenta area covering the range of points with the solar modulation (Φ < 700 MV) and observation constraints (as marked in the upper two panels). We calculate the av-erage of the A/B ratios within the magenta area to be about 1.36 which is in very good agreement with the average ratio derived from observations, 1.33 ± 0.14, further supporting the FD modeling approach employed here. From both observed and modeled results, we have found the relatively larger FD amplitude above the Martian atmosphere (equivalent to that in space) compared to that on the ground. This can be explained by taking into consideration the atmospheric modulation of the GCRs. As particles transverse through the atmosphere, they slow down or stop in the atmosphere due to ionization energy loss or generate secondaries via fragmentation or spallation interactions. Therefore the surface-observed GCR flux consists primarily of primary GCRs above a few hundreds of MeVs. Besides, the heliospheric modulation of GCRs is more effective on the low-energy particles, which are less likely to contribute to the surface radiation. As a result, the FDs on the surface of Mars are a better proxy of the modulation of 1-3 GeV GCRs, which are less depressed by the enhanced heliospheric activities (compared to particles with energy of a few hundreds of MeVs) and therefore the magnitude of the same FD on the Martian surface is smaller than that in orbit. One may argue that the FDs in the GCR flux rate can not be directly compared with those detected in the GCR dose rate and their difference could be due to the difference in quantifying the GCR flux. Therefore we also calculated the dose rate in space and the flux rate on the surface and compared the FD amplitudes in space and on the Martian surface, measured in the same unit of GCRs, dose rate, or flux rate. The results, shown in Figure 4, suggest the same trend as that shown in Figure 3. In fact the difference of the FD amplitude in Mars's orbit compared to that on Mars's surface is even bigger: the ratio ranges between 1.4 and 1.8, suggesting that the atmospheric effect on FD amplitudes is significant.   (1)) and on Mars' surface (Equation (2)), respectively. The bottom panel shows the ratio of the above two FD amplitudes. The magenta areas mark the points with FD amplitude and solar modulation values falling into the observational range. More explanations of the data shown in the figure can be found in the text of Section 3. Figure 3 also shows the change of the FD amplitude as a function of Φ and δ Φ . Under the same value of Φ, FD amplitude increases with δ Φ during stronger solar disturbances, as expected. Moreover, for a given "strength" of solar disturbance, i.e., a fixed δ Φ , if it happens during a less active solar cycle (smaller Φ), the FD amplitude is slightly bigger. This is because under weaker solar activity there are more lower-energy GCR particles (Figure 2), which are more easily modulated by the temporarily enhanced solar modulation. During solar maximum years (high values of Φ), the GCR spectra are already depressed, especially at low energies, and further modulation by the solar eruptions becomes more difficult, leading to a relatively smaller FD for the same δ Φ . However, this does not mean that the FD amplitude is generally bigger during solar minimum than during solar maximum because it depends also on the "strength" of the solar eruptions. As shown by Belov (2008), the overall FD amplitudes are positively-correlated with solar activity.
Concerning the ignored drift effect in the current BON model, Kadokura and Nishida (1986) and Luo X et al. (2017) have modeled FDs to have a larger amplitude under positive polarities than under negative polarities. During a negative polarity epoch, positively charged GCRs enter the inner heliosphere through the equatorial plane due to current-sheet drift. Consequently, the effect of the sheath diffusion barrier is weakened and the FD amplitude is smaller. However, such model predictions are difficult to verify from available observations because it is almost impossible to find ICMEs with the same shock strength, sheath property, and the same orientation and propagation paths during two different solar polarities, with the other heliospheric conditions being the same. In fact, Richardson et al. (1999) analyzed statistically the amplitudes of recurrent FDs caused by CIRs under different polarities and found that the observational pattern is inconsistent with a prediction from the drift model. Therefore the difference between FD amplitudes under opposite heliospheric polarity signs is still a topic under investigation. Alternatively, recent observations by Munini et al. (2018) seem to suggest that the recovery phase of the FD would depend more on heliospheric polarity, which is supported by models (Luo X et al., 2018) that predict that, during a positive epoch, the recovery time of FDs of GeV protons is shorter than that of electrons, and vice versa. As we are not addressing the recovery phase of the FD in the current study, this should not affect our conclusions.

Conclusions
In order to study the space weather environment near Mars using the ground-measured Forbush decreases, it is important to understand and quantify the energy-dependent modulation of the GCR particles not only by pass-by heliospheric disturbances but also by the Martian atmosphere. In this study, we develop a model that quantifies the amplitudes of the Forbush decreases at Marsboth in the interplanetary space near Mars and on the planet's surface -through combining the heliospheric modulation of GCRs and the further atmospheric modification of such modulated GCR spectra. The modeled results are in excellent agree- ment with our previous statistical study of Forbush decreases as measured by MSL on the surface of Mars and by the MAVEN spacecraft in orbit above the planet.
The good quantitative agreement between the modeled results and empirical measurements suggests that our model would, in principle, be a good approach for quantifying the amplitude of FDs during solar eruptions both in space and on the surface of Mars. However, there has not been a physical function describing how the modulation parameter would change, i.e., δ Φ , depending on the properties of the interplanetary transient, an ICME or an SIR. So the current status of our model cannot directly predict the FD amplitude from solar eruptions. But the main advantage of this approach is that we can quantify the relative FD amplitudes as caused by the same interplanetary transient measured by different instruments, which is important to account for when studying the evolution of propagating ICMEs.
On the other hand, it may be possible to derive δ Φ at certain locations based on measured FD amplitudes by different instruments, i.e., by using multiple y-axis values in Figure 3 to derive the x-axis value of the pass-by transient, and to study the δ Φ dependence on the ICME/SIR properties, if in situ plasma measurement is available. A statistical study may potentially give an empirical function of the enhanced modulation parameter and the ICME/SIR properties, which can be used to quantify FD amplitudes from ICME properties or vice verse.
In conclusion, corrections of FDs measured under different conditions to the same condition (same heliospheric modulation before the onset, without the local shielding influence and with the same cutoff energy) should be applied in order to compare the FD amplitudes measured at different locations by different instruments. To make this correction for FDs observed at Mars, our method takes into account the energy-dependent modulation of the GCR particles by the pass-by heliospheric disturbances and includes the local shielding and instrumental response functions.
Only after this correction we can use the FD measured at multiple locations in the heliosphere to describe the properties of the corresponding ICME. With the abundant particle measurements in space carried out by almost every heliospheric and planetary mission, we open a new window to better understand ICMEs and their evolution during their propagation in the interplanetary space.