Simulations of ⁷Be and ¹⁰Be with the GEOS-Chem global model v14.0.2 using state-of-the-art production rates

The cosmogenic radionuclides ⁷Be and ¹⁰Be are useful tracers for atmospheric transport studies. Combining ⁷Be and ¹⁰Be measurements with an atmospheric transport model can not only improve our understanding of the radionuclide transport and deposition processes but also provide an evaluation of the transport process in the model. To simulate these aerosol tracers, it is critical to evaluate the influence of radionuclide production uncertainties on simulations. Here we use the GEOS-Chem chemical transport model driven by the Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2) reanalysis to simulate ⁷Be and ¹⁰Be with the state-of-the-art production rate from the CRAC:Be (Cosmic Ray Atmospheric Cascade: Beryllium) model considering realistic spatial geomagnetic cutoff rigidities (denoted as P16spa). We also perform two sensitivity simulations: one with the default production rate in GEOS-Chem based on an empirical approach (denoted as LP67) and the other with the production rate from the CRAC:Be but considering only geomagnetic cutoff rigidities for a geocentric axial dipole (denoted as P16). The model results are comprehensively evaluated with a large number of measurements including surface air concentrations and deposition fluxes. The simulation with the P16spa production can reproduce the absolute values and temporal variability of ⁷Be and ¹⁰Be surface concentrations and deposition fluxes on annual and sub-annual scales, as well as the vertical profiles of air concentrations. The simulation with the LP67 production tends to overestimate the absolute values of ⁷Be and ¹⁰Be concentrations. The P16 simulation suggests less than 10 % differences compared to P16spa but a significant positive bias (∼18 %) in the ⁷Be deposition fluxes over East Asia. We find that the deposition fluxes are more sensitive to the production in the troposphere and downward transport from the stratosphere. Independent of the production models, surface air concentrations and deposition fluxes from all simulations show similar seasonal variations, suggesting a dominant meteorological influence. The model can also reasonably simulate the stratosphere–troposphere exchange process of ⁷Be and ¹⁰Be by producing stratospheric contribution and ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{7}}\mathrm{Be}$ ratio values that agree with measurements. Finally, we illustrate the importance of including the time-varying solar modulations in the production calculation, which significantly improve the agreement between model results and measurements, especially at mid-latitudes and high latitudes. Reduced uncertainties in the production rates, as demonstrated in this study, improve the utility of ⁷Be and ¹⁰Be as aerosol tracers for evaluating and testing transport and scavenging processes in global models. For future GEOS-Chem simulations of ⁷Be and ¹⁰Be, we recommend using the P16spa (versus default LP67) production rate.

1 Introduction

The naturally occurring cosmogenic radionuclide ⁷Be (half-life of 53.2 d) is monitored worldwide and has been recognized as a useful tracer in atmospheric dynamic studies (Aldahan et al., 2001; Hernández-Ceballos et al., 2016; Terzi et al., 2019; Liu et al., 2016). In particular, ratios of radionuclides concentrations with very different half-lives, such as the ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{7}}\mathrm{Be}$ ratio, have become powerful tools (e.g., Liu et al., 2022b; Raisbeck et al., 1981) to disentangle the influence of transport and deposition since both ⁷Be and ¹⁰Be in the troposphere are mainly removed by wet deposition. In this paper, we aim to improve the utility of ⁷Be and ¹⁰Be as tracers for atmospheric transport using state-of-the-art production rates in a global 3-D chemical transport model.

⁷Be and ¹⁰Be are produced through interactions between atmospheric atoms (mostly oxygen and nitrogen) and incoming cosmic rays in the atmosphere (Lal and Peters, 1967, referred to as LP67 hereafter; Poluianov et al., 2016, referred to as P16 hereafter). Due to the atmospheric depth profile of fluxes of primary cosmic rays, the formed secondary particles, and their energy, ⁷Be and ¹⁰Be production rates reach their maxima in the lower stratosphere (Poluianov et al., 2016). About two-thirds of ⁷Be and ¹⁰Be are produced in the stratosphere, while the rest is produced in the troposphere (Poluianov et al., 2016; Heikkilä and Smith, 2013; Golubenko et al., 2022). Once produced, ⁷Be and ¹⁰Be rapidly attach to aerosol particles and get transported and deposited with their carrier aerosols by wet and dry deposition (Delaygue et al., 2015; Heikkilä et al., 2013). ¹⁰Be has a half-life of 1.39×10⁶ years (Chmeleff et al., 2010), and its decay is thus negligible compared to its average atmospheric residence time (about 1–2 years) (Heikkilä et al., 2008b). During transport away from the regions of their production, the ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{7}}\mathrm{Be}$ ratio increases because ⁷Be decays. The ratio ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{7}}\mathrm{Be}$ therefore could indicate the path-integrated age of the air mass. Due to different aerosol residence times in the stratosphere (more than 1 year) and troposphere (∼ weeks), the ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{7}}\mathrm{Be}$ ratio is higher in the stratosphere than in the troposphere. Hence, the ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{7}}\mathrm{Be}$ ratio can be used to detect the stratosphere–troposphere exchange.

Many studies have focused on understanding the signals in surface ⁷Be measurements from worldwide monitoring stations (e.g., Hernandez-Ceballos et al., 2015; Rodriguez-Perulero et al., 2019; Uhlar et al., 2020; Ajtić et al., 2022; Burakowska et al., 2021). Due to the cosmogenic origin of ⁷Be, surface air ⁷Be concentrations are found to be connected to the 11-year cycle of solar modulation (Leppänen et al., 2010; Zheng et al., 2021b). In addition, ⁷Be concentrations in the surface air are affected by different meteorological processes depending on locations, such as stratospheric intrusions (Jordan et al., 2003; Pacini et al., 2015; Yamagata et al., 2019), scavenging by precipitation (Chae and Kim, 2019; Kusmierczyk-Michulec et al., 2015), vertical transport in the troposphere (Aldahan et al., 2001; Ajtić et al., 2018; Zheng et al., 2021b), and large-scale atmospheric circulations (Hernández-Ceballos et al., 2022; Terzi and Kalinowski, 2017).

The ability of general circulation models (GCMs; e.g., GISS ModelE, ECHAM5-HAM, and EMAC) and chemical transport models (CTMs; e.g., GEOS-Chem and GMI) to capture the main characteristics in ⁷Be and ¹⁰Be transport and deposition has been demonstrated in previous studies (e.g., Heikkilä et al., 2008b; Koch and Rind, 1998; Field et al., 2006; Usoskin et al., 2009; Brattich et al., 2021; Spiegl et al., 2022; Liu et al., 2016; Sukhodolov et al., 2017). For example, Usoskin et al. (2009) found that the influence of the solar proton-induced ⁷Be production peak at the surface in early 2005 is small through the comparison of GISS ModelE simulations and surface air measurements. Heikkilä et al. (2009) showed that stratospheric ¹⁰Be contribution is dominant in the global ¹⁰Be deposition by tracing tropospheric and stratospheric ¹⁰Be separately in the aerosol–climate model ECHAM5-HAM. Spiegl et al. (2022) used the EMAC climate model to investigate the transport and deposition process of ¹⁰Be produced by the extreme solar proton event in 774–775 CE. They suggested that the downward transport of ¹⁰Be from the stratosphere is mainly controlled by the Brewer–Dobson circulation in the stratosphere and cross-tropopause transport. By comparing the measurements with GEOS-Chem simulations over January–March 2003, Brattich et al. (2021) found that increased ⁷Be values in surface air samples in northern Europe in early 2003 were associated with the instability of the Arctic polar vortex. They also showed that, while the model generally simulates the month-to-month variation in surface ⁷Be concentrations well, it tends to underestimate the observations (see their Table 2), partly due to the use of the default LP67 production rate for a solar maximum year (1958) in the GEOS-Chem model (Liu et al., 2001). Using the GMI CTM driven with four different meteorological datasets, Liu et al. (2016) showed that the observational constraints for ⁷Be and observed ⁷Be total deposition fluxes can be used to provide a first-order assessment of cross-tropopause transport in global models. In comparison to GCMs with or without nudged winds (e.g., Golubenko et al., 2021; Heikkilä et al., 2008b; Spiegl et al., 2022) which involve simulating the entire global circulation and climate, the “offline” CTMs are driven by archived meteorological datasets, either from output of GCMs or from atmospheric data assimilation systems. For example, GEOS-Chem can be driven by the GEOS assimilated meteorology (e.g., MERRA-2 reanalysis data; Gelaro et al., 2017a) or output from the GISS GCM (e.g., Murray et al., 2021).

In comparison with the LP67 production rate using an empirical approach (Lal and Peters, 1967; Liu et al., 2001; Brattich et al., 2021), the recent production models apply full Monte Carlo simulations of the cosmic-ray-induced atmospheric nucleonic cascade (e.g., Poluianov et al., 2016; Masarik and Beer, 1999). LP67 shows the highest ⁷Be and ¹⁰Be production rates compared to other production models (Elsässer, 2013). P16 suggests that LP67 overestimates the ⁷Be production rate by 30 %–50 % compared to their production model (Poluianov et al., 2016). Furthermore, the LP67 production rate implemented in GEOS-Chem is only validated for the year 1958, a year with a high solar modulation function (i.e., high solar activity) of 1200 MeV (Herbst et al., 2017). This highlights the problem of quantitatively comparing these uncorrected model outputs with measurements from other periods. Some studies (e.g., Koch et al., 1996; Liu et al., 2016) have applied a scale factor to account for this solar modulation influence on LP67 production rate. However, this correction is not ideal as the influence of varying solar modulation is latitudinally and vertically dependent. In earlier studies, the ¹⁰Be production rate in GEOS-Chem was simply scaled to the ⁷Be production rate based on the ratio estimated from the surface measurements (Koch and Rind, 1998). In addition, ¹⁰Be as simulated by GEOS-Chem has not been evaluated so far. It is hence necessary to update the ⁷Be and ¹⁰Be production rates in GEOS-Chem and assess the corresponding impacts on model simulation results.

In this study, we incorporate global ⁷Be and ¹⁰Be production rates from the recently published CRAC:Be (Cosmic Ray Atmospheric Cascade: Beryllium) model (Poluianov et al., 2016) into the GEOS-Chem model. We simulate ⁷Be and ¹⁰Be using GEOS-Chem with the following three production scenarios:

• Scenario I, a production rate derived from the CRAC:Be model considering realistic geomagnetic cutoff rigidity (P16spa production rate);

• Scenario II, a production rate derived from the CRAC:Be model considering an approximation of geomagnetic cutoff rigidities using a geocentric axial dipole (P16 production rate);

• Scenario III, a default production rate in GEOS-Chem using an empirical approximation (LP67 production rate).

Scenario I is treated as the standard simulation, while the other two are sensitivity tests that also enable comparison to earlier studies. This paper is organized as follows. Section 2 introduces the GEOS-Chem model and three different ⁷Be and ¹⁰Be production rates, discusses the methodology and experiment design, and describes the observational data for model evaluations. In Sect. 3, we first investigate the differences between three different production scenarios (Sect. 3.1). Then, we evaluate model simulations of ⁷Be and ¹⁰Be with several published datasets of ⁷Be and ¹⁰Be measurements, in terms of absolute values (Sect. 3.2–3.3), vertical profiles (Sect. 3.4), and seasonal variations (Sect. 3.6). The budgets and residence times of ⁷Be and ¹⁰Be are given in Sect. 3.5. We also examine the ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{7}}\mathrm{Be}$ ratio in the model to assess its ability in capturing the stratosphere–troposphere exchange (Sect. 3.7). Finally, we investigate the influence of including solar-induced production rate variability on ⁷Be simulations (Sect. 3.8). A summary and conclusions are given in Sect. 4.

2 Models and data

2.1 GEOS-Chem model

GEOS-Chem is a global 3-D chemical transport model (http://www.geos-chem.org, last access: 23 November 2023) that simulates gases and aerosols in both the troposphere and stratosphere (Eastham et al., 2014; Bey et al., 2001). It is driven by archived meteorological data. We use version 14.0.2 (https://wiki.seas.harvard.edu/geos-chem/index.php/GEOS-Chem_14.0.2, last access: 23 November 2023) to simulate the transport and deposition of atmospheric ⁷Be and ¹⁰Be. We drive the model with the Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2) meteorological reanalysis (http://gmao.gsfc.nasa.gov/reanalysis/MERRA-2/, last access: 23 November 2023; Gelaro et al., 2017b). MERRA-2 has a native resolution of 0.5^∘ latitude by 0.667^∘ longitude, with 72 vertical levels up to 0.01 hPa (80 km). Here the MERRA-2 data are re-gridded to 4^∘ latitude by 5^∘ longitude for input to GEOS-Chem for computational efficiency.

GEOS-Chem includes a radionuclide simulation option (²²²Rn–²¹⁰Pb–⁷Be–¹⁰Be), which simulates transport (advection, convection, boundary-layer mixing), deposition, and decay of the radionuclide tracers (e.g., Liu et al., 2001, 2004; B. Zhang et al., 2021; Yu et al., 2018). The model uses the TPCORE algorithm of Lin and Rood (1996) for advection, archived convective mass fluxes to calculate convective transport (Wu et al., 2007), and the non-local scheme implemented by Lin and McElroy (2010) for boundary-layer mixing. As mentioned in the Introduction, the standard GEOS-Chem model uses the LP67 ⁷Be and ¹⁰Be production rates. After production, ⁷Be and ¹⁰Be attach to ambient submicron aerosols ubiquitously, and their behavior becomes that of aerosols until they are removed by wet-deposition (precipitation scavenging) and dry-deposition processes. Note that neither is the process of attachment explicitly represented, nor is the aerosol size distribution considered in the model. In addition, the decay process is included for the short-lived ⁷Be with a half-life of 53.2 d. The decay is minor for the long-living ¹⁰Be, which has a half-life of 1.39×10⁶ years (e.g., Chmeleff et al., 2010).

Wet deposition includes rainout (in-cloud scavenging) due to stratiform and anvil precipitation (Liu et al., 2001), scavenging in convective updrafts (Mari et al., 2000), and washout (below-cloud scavenging) by precipitation (Wang et al., 2011). Scavenged aerosols from vertical layers above are allowed to be released into the atmosphere during the re-evaporation of precipitation below the cloud. In the case of partial re-evaporation, we assume that half of the corresponding fraction of the scavenged aerosol mass is released at that level because some of the re-evaporation of precipitation are due to partial shrinking of the raindrops, which does not release aerosol (Liu et al., 2001). MERRA-2 fields of precipitation formation and evaporation are used directly by the model wet deposition scheme. Dry deposition is based on the resistance-in-series scheme of Wesely (1989). The process of sedimentation is not included in the model.

To quantify the stratospheric contribution to ⁷Be and ¹⁰Be in the troposphere, we separately transport ⁷Be and ¹⁰Be produced in the model layers above the MERRA-2 thermal tropopause (i.e., stratospheric ⁷Be and ¹⁰Be tracers). This approach was previously used to study cross-tropopause transport of ⁷Be in GEOS-Chem (Liu et al., 2001; Brattich et al., 2021) and Global Modeling Initiative chemical transport models (Liu et al., 2016; Brattich et al., 2017). The stratospheric fractions of ⁷Be and ¹⁰Be are defined as the ratio of the stratospheric ⁷Be and ¹⁰Be concentrations to the ⁷Be and ¹⁰Be concentrations.

2.2 ⁷Be and ¹⁰Be production models

The GEOS-Chem currently uses the LP67 production rates of ⁷Be and ¹⁰Be (Lal and Peters, 1967). These production rates are calculated using an analytically estimated rate of nuclear disintegration (stars) in the atmosphere (stars per gram of air per second), multiplied by the mean production yield of 0.045 atoms per star for ⁷Be and 0.025 atoms per star for ¹⁰Be (Lal and Peters, 1967). These rates are represented as a function of latitude and altitude for the year 1958 and are not time-varying.

Here we update the atmospheric ⁷Be and ¹⁰Be production rates in GEOS-Chem with the latest production model: the CRAC:Be model by P16 (Poluianov et al., 2016) using the solar modulation function record by Herbst et al. (2017). The solar modulation function record is based on the local interstellar spectrum by Herbst et al. (2017), which was also used in the production model. Given spatially and temporally resolved geomagnetic cutoff rigidities, the P16 model allows the calculation of 3-dimensional, temporally variable ⁷Be and ¹⁰Be production rates, which are necessary for input to atmospheric transport models. The P16 production model is regarded as the latest and one of the most accurate production models for ⁷Be and ¹⁰Be and was used in recent general circulation model simulations (e.g., Golubenko et al., 2021; Sukhodolov et al., 2017).

The production rates of ⁷Be and ¹⁰Be are calculated as Eq. (1) by an integral of the yield functions of ⁷Be and ¹⁰Be (Y_i, atoms g⁻¹ cm² sr) and the energy spectrum of cosmic rays (J_i, (sr s cm²)⁻¹) above the cutoff energy E_c:

$$\begin{array}{}\text{(1)}& Q\left(\mathrm{\Phi },h,P_{\mathrm{c}}\right)=\sum _i\underset{E_{\mathrm{c}}}{\overset{\mathrm{\infty }}{\int }}{Y}_i\left(E,h\right)J_i\left(E,\mathrm{\Phi }\right)\mathrm{d}E,\end{array}$$

where i refers to different types of primary cosmic ray particles (e.g., proton, alpha, and heavier particles). For modeling the contribution of alpha and heavier particles to the total production, their nucleonic ratio in the local interstellar spectrum was set to 0.353 (Koldobskiy et al., 2019). The yield function Y_i is a function of height (h) and kinetic energy per incoming primary nucleon (E) and is directly taken from P16. The energy spectrum of cosmic rays J_i is a function of the kinetic energy (E) and depends on the solar modulation function (Φ) (Herbst et al., 2017). E_c is calculated as Eq. (2) as a function of the local geomagnetic rigidity cutoff (P_c):

$$\begin{array}{}\text{(2)}& E_{\mathrm{c}}=E_{\mathrm{r}}(\sqrt{\mathrm{1}+{\left({\displaystyle \frac{Z_i{P}_{\mathrm{c}}}{A_i{E}_{\mathrm{r}}}}\right)}^{\mathrm{2}}}-\mathrm{1}),\end{array}$$

where Z_i and A_i are the charge and mass numbers of particles, respectively. E_r is the rest mass of a proton (0.938 GeV).

The geomagnetic rigidity cutoff P_c is a quantitative estimation of the Earth's geomagnetic field shielding effect (Smart and Shea, 2005). Cosmic ray particles with rigidity (momentum per unit charge of the particle) higher than the geomagnetic cutoff rigidity value can enter the Earth's atmosphere. In several model simulations of ⁷Be and ¹⁰Be (e.g., Field et al., 2006; Koch et al., 1996; Liu et al., 2001), the production is calculated with a P_c simplified as a function of the geomagnetic latitude and geomagnetic dipole moment, called the vertical Störmer cutoff rigidity equation (see Eqs. 5.8.2–2 in Beer et al., 2012). However, this is different from the real geomagnetic cutoff rigidity inferred from the trajectories of particles with different energies using real geomagnetic field measurements (e.g., Copeland, 2018), which also includes non-dipole moments of the field (Beer et al., 2012) (Fig. S1 in the Supplement). Earlier studies suggested that using the simple centered dipole models (e.g., Störmer cutoff rigidity) for cutoff rigidity approximation is limited as they can significantly distort the cutoff rigidity for some regions (e.g., low-latitude regions) (Pilchowski et al., 2010; Nevalainen et al., 2013).

Here we take the geomagnetic cutoff rigidity from Copeland (2018) that provides the cutoff rigidity at a fine interval (1^∘) in both latitude and longitude. This production rate is denoted as P16spa. To investigate the effect of this more realistic representation of cutoff rigidity on ⁷Be and ¹⁰Be simulations, we also perform simulations where the cutoff rigidities are approximated by the Störmer equation (denoted as P16). The influence of the geomagnetic field intensity variations can be considered negligible on annual and decadal timescales and are ignored here (e.g., Muscheler et al., 2007; Zheng et al., 2020). It should be mentioned that the LP67 production is based on an ideal axial dipole cutoff rigidity similar to the P16 production model.

2.3 GEOS-Chem model experiments and evaluations

An overview of the performed simulations is shown in Table S1 in the Supplement. The simulation with the P16spa production rate is considered the standard simulation, while the simulations with the P16 and LP67 production rates are sensitivity tests. The simulation with the P16 production rate is conducted to evaluate the influence of a simplified approximation of cutoff rigidities resulting from a geocentric dipole. In earlier studies, the LP67 production rate was used for global model simulations of ⁷Be (e.g., Liu et al., 2016; Brattich et al., 2017; Liu et al., 2001; Koch et al., 1996). The purpose of performing the simulation with the LP67 production rate is to evaluate to what extent model simulations are biased when applying the default LP67 production. Since the LP67 production rate applies only for the year 1958 (with a solar modulation function of about 1200 MeV) and does not consider the influences of the solar variations (e.g., 11-year solar cycle), it underestimates the production rate for the period of 2008–2018 that has an average solar modulation function of 500 MeV. To correct for this solar modulation influence, we follow the previous studies (e.g., Liu et al., 2016; Koch et al., 1996) by multiplying the model results by a scale factor of 1.39. It should be noted that this correction is not ideal as the effects of a varying solar modulation on cosmogenic radionuclide production rates depend on altitude and latitude. All simulations are performed from 2002 to 2018, with the first 6 years for spin-up to make sure the ¹⁰Be nearly reaches equilibrium in the atmosphere and the 2008–2018 period (11 years) for analysis. The simulations are conducted using a 4^∘ latitude × 5^∘ longitude resolution for computational efficiency (e.g., Liu et al., 2016, 2004).

To evaluate the model's ability to reproduce the variabilities in the observations, we use the statistical parameters: Spearman correlation coefficients and root mean square error (RMSE) (Chang and Hanna, 2004). Spearman's rank correlation (R) (Myers et al., 2013) is used as it does not make any assumptions about the variables being normally distributed. It is less sensitive to outliers in the data compared to the commonly used Pearson correlation. The fraction of modeled concentrations within a factor of 2 of observations (FA2) is calculated, i.e., for which $\mathrm{0.5}<X_{\mathrm{model}}/X_{\mathrm{observation}}<\mathrm{2}$. Usually, if the scatter plot of the model and measurements is within a factor of 2 of observations, the model is considered to have a reasonably good performance (e.g., Heikkilä et al., 2008b; Brattich et al., 2021). For model comparison with surface air concentrations, the model value from the bottom grid box closest to the corresponding measurement site is selected.

2.4 ⁷Be and ¹⁰Be observational data for model validation

The annual mean ⁷Be surface air concentration and deposition measurements are taken from a compilation by F. Zhang et al. (2021). The compilation includes a total of 494 annual mean values for surface air ⁷Be concentrations and 304 for ⁷Be deposition fluxes. For the deposition measurements, most of them include both wet and dry deposition, while a few are collected only during rainfall events and thus include only wet deposition. It includes the data from the following:

• the Environmental Measurements Laboratory (EML; https://www.wipp.energy.gov/namp/emllegacy/index.htm, last access: 23 November 2023) Surface Air Sampling Program (SASP), which began in the 1980s;

• the ongoing international monitor program Radioactivity Environmental Monitoring (REM) network (e.g., Hernandez-Ceballos et al., 2015; Sangiorgi et al., 2019);

• the International Monitoring System (IMS) organized by the Comprehensive Nuclear-Test-Ban Treaty Organization (CTBTO) (e.g., Terzi and Kalinowski, 2017);

• some additional datasets in publications not included in the above programs.

We only include the data covering more than 1 year to reduce the influence of inherent seasonal variations. We further include several recently published data for ⁷Be surface air concentrations and deposition fluxes records that cover more than 1 year (Burakowska et al., 2021; Liu et al., 2022b; Kong et al., 2022).

The data used for investigating the seasonality of ⁷Be surface air concentrations are mainly taken from a multiyear compilation dataset of IMS from Terzi and Kalinowski (2017). The seasonal ⁷Be deposition data are taken from Courtier et al. (2017), Du et al. (2015), Dueñas et al. (2017), Hu et al. (2020), Lee et al. (2015), and Sangiorgi et al. (2019). The vertical profile of ⁷Be concentrations is taken from the Environmental Measurements Laboratory (EML) High Altitude Sampling Program (HASP), spanning the years of 1962–1983. It should be noted that, different from surface air measurements, the vertical air samples were usually collected during single-day flight campaigns.

There are fewer ¹⁰Be measurements compared to ⁷Be. Here we compiled two datasets of published ¹⁰Be surface air measurements (Table S2) (Aldahan et al., 2008; Liu et al., 2022a; Yamagata et al., 2019; Padilla et al., 2019; Rodriguez-Perulero et al., 2019; Huang et al., 2010; Méndez-García et al., 2022; Elsässer et al., 2011; Dibb et al., 1994) and deposition fluxes (Table S3) covering more than 1 year, to validate the model performance. The air samples are continuously collected by filters using a high-flow aerosol sampler. The sampling volume is approximately 700 m³ of air for daily samples (e.g., Liu et al., 2022a) and between 3000 and 5000 m³ for weekly samples (e.g., Yamagata et al., 2019). The deposition data include the precipitation samples (wet deposition) (Graham et al., 2003; Monaghan et al., 1986; Somayajulu et al., 1984; Heikkilä et al., 2008a; Raisbeck et al., 1979; Maejima et al., 2005) and ice core samples (wet and dry deposition) that cover the recent period (Heikkilä et al., 2008a; Zheng et al., 2021a; Pedro et al., 2012; Baroni et al., 2011; Aldahan et al., 1998; Berggren et al., 2009; Auer et al., 2009; Zheng et al., 2023a). The ¹⁰Be vertical profile measurements are mainly taken from Dibb et al. (1992, 1994) and Jordan et al. (2003).

3 Results and discussions

3.1 ⁷Be and ¹⁰Be production rates

Figure S2 shows the comparison between ⁷Be_P16 and ⁷Be_LP67 production rates for the year 1958. Generally, the ⁷Be_P16 production rate shows a similar production distribution to the ⁷Be_LP67 production rate, with a maximum ⁷Be production over the polar stratosphere (∼100 hPa). The ⁷Be_LP67 production rate shows, on average, about a 72 % higher production rate compared to ⁷Be_P16 in the stratosphere and about 38 % in the troposphere (Fig. S2c; Table S4). On a global average, the ⁷Be_LP67 production rate is about 60 % higher than that of ⁷Be_P16, as shown in previous studies (Poluianov et al., 2016). The stratospheric production contributes about 67 % to the total production for the ⁷Be_LP67 production rate, while it is about 62 % for the ⁷Be_P16 production rate for the year 1958.

The ¹⁰Be_LP67 production rate in the GEOS-Chem model uses the identical source distribution as ⁷Be with a scaling factor based on the estimates from surface air measurements (Koch and Rind, 1998). This leads to a constant ${}^{\mathrm{10}}{\mathrm{Be}}_{\mathrm{LP}\mathrm{67}}{/}^{\mathrm{7}}{\mathrm{Be}}_{\mathrm{LP}\mathrm{67}}$ production ratio (0.55) throughout the entire atmosphere. However, as shown in many ⁷Be and ¹⁰Be production models (e.g., Poluianov et al., 2016; Masarik and Beer, 2009), ⁷Be and ¹⁰Be have different altitudinal production distributions. The P16 production shows an increasing ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{7}}\mathrm{Be}$ production ratio from higher altitudes (0.35) to lower altitudes (0.6) (Fig. S3). Using a constant ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{7}}\mathrm{Be}$ production ratio may thus result in large errors in the modeled ¹⁰Be concentrations as well as ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{7}}\mathrm{Be}$ ratios. The stratospheric production contributes about 67 % of the total production with ¹⁰Be_LP67, while it is about 58 % with the ¹⁰Be_P16 production for the year 1958 (Table S4).

Figure 1 Spatial distribution of (a) ⁷Be_P16spa and (b) ⁷Be_P16 production rates at 825 hPa over the period 2008–2018. (c) Relative differences (%), i.e., (⁷Be_P16spa-⁷Be${}_{\mathrm{P}\mathrm{16}}){/}^{\mathrm{7}}$Be_P16×100 %, between production rates with and without considering the detailed spatial cutoff rigidity. (d) Relative differences (%) of the zonal mean production rates between P16spa and P16 at 30^∘ N.

Figure 1 shows the comparison between ⁷Be_P16 and ⁷Be_P16spa production rates for the period 2008–2018. The global production is similar for P16spa and P16 (Table S4). However, considering non-dipole moment influence on geomagnetic cutoff rigidity, ⁷Be_P16spa and ¹⁰Be_P16spa production rates in the Southern Hemisphere show production rates that are ∼11 % higher compared to the Northern Hemisphere (Table S4). This difference is not present when an axial dipole is assumed. Compared to the P16 production rate, the ⁷Be_P16spa production rate is 30 %–40 % lower over eastern Asia and southeastern Pacific but 40 %–50 % higher over North America and from subtropical South Atlantic to Australia (Fig. 1). ¹⁰Be_P16spa shows similar results to ⁷Be_P16spa. These differences are not constant throughout the atmospheric column but generally increase with altitude (Fig. 1d).

Figure 2 Scatter plot of modeled ⁷Be_P16spa versus observed ⁷Be surface air concentrations (a) and deposition fluxes (b). The model values are averaged over the years of 2008–2018. The dashed lines are the factor of 2 of 1:1 line (straight lines). The “FA2” label indicates the fraction of modeled concentrations within a factor of 2 of observations, while “RMSE” indicates the root mean square error.

3.2 ⁷Be surface air concentrations and deposition fluxes

Figure 2 compares the simulated ⁷Be_P16spa averaged over 2008–2018 with the measurements. Due to the data availability, the measurements do not necessarily cover the same period as model simulations. The model deposition fluxes here include both dry and wet deposition. About 93.7 % of modeled air ⁷Be_P16spa concentrations agree within a factor of 2 with the observed values. The model also shows reasonable agreement with the measured deposition fluxes (60.9 % within a factor of 2), although the discrepancy between the modeled and observed deposition fluxes is larger than that for surface air concentrations. The deposition fluxes are usually less well monitored compared to the air ⁷Be samples and cover only shorter periods (e.g., 1 or 2 years). Further, the limited model resolution applied here may not be able to capture meteorological conditions on local scales (e.g., precipitation, convection, and tropopause folding) in some sites (e.g., Yu et al., 2018; Spiegl et al., 2022), especially for coastal regions when the sub-grid-scale orographic precipitation is important.

Figure 3 (a) Modeled ⁷Be_P16spa surface air concentrations (mBq m⁻³) and (b) deposition fluxes (Bq m⁻² yr⁻¹) averaged over the period 2008–2018. Color-coded dots denote ⁷Be measurements. Zonal mean of (c) observed ⁷Be surface air concentrations and (d) deposition fluxes (black lines, for each 5^∘ latitude bin) compared with the model simulation using the P16spa production rate (blue lines). Dots are individual measurements. The error bars indicate 1 standard deviation. The outliers, defined as more than 3 scaled median absolute deviations (MADs) away from the median, are excluded from the calculation. The observations are averaged over the years available.

Figure 3 shows the spatial distribution and zonal mean of measurements in comparison with the model-simulated ⁷Be_P16spa surface air concentrations and deposition fluxes. Generally, the model captures the spatial distribution of ⁷Be air concentrations and deposition fluxes. The “latitudinal pattern” of surface air ⁷Be concentrations differs from that of the ⁷Be production rate, reflecting the effects of atmospheric transport and deposition processes. The model suggests high ⁷Be air concentrations, mainly over the dry regions (Fig. 3a), due to low wet deposition rates (e.g., desert regions over northern Africa, Arabian Peninsula, central Australia, and central Antarctica), and over high-altitude regions (e.g., Tibetan Plateau). The model captures the observed latitudinal peaks in surface air concentrations over the subtropics and mid-latitudes (Fig. 3c around 30–40^∘ N and 30–40^∘ S). These peaks are consistent with the high stratospheric contribution (25 %–30 %) at mid-latitudes (Fig. S4). The model overestimates ⁷Be air concentrations over the Arctic (70–90^∘ N; Fig. 3c) by about 30 %–40 %. By contrast, high ⁷Be deposition fluxes are observed at mid-latitudes due to the influence of the high precipitation (wet deposition) and strong stratosphere–troposphere exchange (Fig. 3d). In the Northern Hemisphere, the model-simulated deposition fluxes peak at a lower latitude (∼30^∘ N) relative to the observations (∼45^∘ N). These modeled spatial distributions of the air concentrations and deposition rates of ⁷Be also agree generally well with previous model simulations (e.g., Heikkilä and Smith, 2012).

The modeled ⁷Be_P16spa air concentrations show better agreements (smaller RMSE and higher FA2 values) with the measurements in comparison to ⁷Be_LP67 (Fig. S5). ⁷Be_LP67 tends to overestimate the absolute values of ⁷Be concentrations. This is caused by (i) the overestimation of the ⁷Be production rate by LP67 for a given solar modulation function and (ii) the use of a simple scale factor to account for the solar modulation influence on the LP67 ⁷Be production rate.

Figure 4 Relative differences (percentage) of surface air concentrations (a) and deposition fluxes (b) between ⁷Be_P16spa and ⁷Be_P16 for the period 2008–2018, i.e., (⁷Be_P16spa-⁷Be${}_{\mathrm{P}\mathrm{16}}){/}^{\mathrm{7}}$Be_P16×100 %.

We also examine whether using the dipole approximation of the cutoff rigidity or real cutoff rigidity (P16 and P16spa, respectively) in the production model leads to significantly different results (Fig. 4). Although large regional differences (up to 40 %–50 %, Fig. 1) in the production model are observed between P16spa and P16 production rates, such differences are reduced in surface air concentrations and deposition fluxes due to transport and deposition processes, as expected. The ⁷Be_P16sap air concentrations show higher values (∼7 %) over 10–40^∘ S and lower values (∼12 %) over the East Asian region (Fig. 4a) compared to ⁷Be_P16. These differences are higher for the deposition fluxes, up to 10 % higher over the 10–40^∘ S and up to 18 % lower over the east Asian region (Fig. 4b). Since the total deposition flux reflects precipitation scavenging through the tropospheric column, it tends to be more sensitive to ⁷Be air concentrations at higher altitudes and downward transport of ⁷Be from the stratosphere. Indeed, model results suggest that deposition fluxes have a higher stratospheric fraction compared to surface air concentrations (Fig. S4), as previously shown by Liu et al. (2016). The ⁷Be_P16spa deposition fluxes show better agreement with measurements than those of ⁷Be_P16 (Fig. S5). The comparison for ¹⁰Be shows similar results to ⁷Be except with less than 10 % difference. For ¹⁰Be deposition fluxes in Antarctica and Greenland, this influence is less than 3 %. This is because the dominant contribution of ¹⁰Be is from the stratosphere where the hemispheric production differences are diminished by the long stratospheric residence time of ¹⁰Be. However, it does not suggest that the cutoff rigidity including the non-dipole influence could be ignored for ¹⁰Be deposition in polar regions, as the spatial pattern of cutoff rigidities was very different in the past, e.g., during the Laschamps geomagnetic field minimum around 41 000 years before the present (Gao et al., 2022). Further studies are warranted to investigate this spatial cutoff rigidity influence on ¹⁰Be in more detail.

Figure 5 Modeled annual mean ¹⁰Be_P16spa (a) surface air concentrations and (b) deposition fluxes averaged over 2008–2018 overplotted with measurements (color-coded dots). (c–d) Scatter plot between model results and measurements for (c) surface air concentrations and (d) deposition fluxes. The dots in (c–d) indicate measurements with careful examination of dust ¹⁰Be contributions or from the polar regions which are not influenced by dust ¹⁰Be. The crosses indicate the samples without examining dust contributions. The FA2 and RMSE are calculated only using the dust-free samples (dots). Blue and orange colors indicate the results using P16spa and LP67 production rates, respectively.

3.3 ¹⁰Be surface air concentrations and deposition fluxes

Figure 5 shows the comparison between modeled annual mean ¹⁰Be_P16spa surface air concentrations (or deposition fluxes) averaged over 2008–2018 and measurements. The ¹⁰Be_P16spa shows similar spatial distributions as ⁷Be_P16spa because both radionuclides share the same transport and deposition processes. The model underestimates the measured ¹⁰Be surface air concentrations and deposition fluxes at some sites (Fig. 5b, d). This may be attributed to the influence of resuspended dust with ¹⁰Be attached, which could typically contribute 10 %–35 % to the air ¹⁰Be concentrations (Monaghan et al., 1986). It should be mentioned that ⁷Be decays in the dust because of its short half-life and therefore does not contribute to the surface air ⁷Be concentrations. Indeed, data for which a careful examination of the recycled dust ¹⁰Be in samples was conducted (e.g., Monaghan et al., 1986) or from locations that are less influenced by recycled dust ¹⁰Be (e.g., polar regions; dots in Fig. 5b–d) show better agreement with the model simulations. This suggests the importance of considering the dust contribution when measuring the air ¹⁰Be samples. The model also shows relatively good agreement with most ¹⁰Be deposition data from polar ice cores (marked as dots in Fig. 5d) within a factor of 2.

Figure 6 Comparison of the vertical profile between measurements (circles) and model zonal mean ⁷Be_P16spa and ¹⁰Be_P16spa concentrations for each latitudinal band (15^∘) over the period 2008–2018. The ⁷Be (circle with error bar) observations (from the EML/HASP) are averaged for the altitude band of every 2 km where more than five samples are available. We exclude the outlier from the calculation, which is defined as more than 3 scaled median absolute deviations (MADs) away from the median. The ¹⁰Be profile measurements are mainly taken from Dibb et al. (1992, 1994) and Jordan et al. (2003).

3.4 Vertical profiles of ⁷Be and ¹⁰Be

Figure 6 shows the simulated annual zonal mean vertical profiles of ⁷Be_P16spa and ¹⁰Be_P16spa concentrations compared with those from aircraft measurements in the troposphere and stratosphere from the EML/HASP. The measurements cover different regions and specific meteorological conditions; hence they should only provide a range in which the model results should lie. Following previous modeling studies (Heikkilä et al., 2008b; Koch et al., 1996), we compare model zonal mean values in each 15^∘ latitude band with the corresponding observations.

The simulated ⁷Be_P16spa profiles agree well with the measurements, especially capturing the peaks at ∼20–22 km at mid-latitudes and low latitudes (e.g., Fig. 6c, e, h). The feature that ⁷Be increases with altitude without a peak at 22 km at northern high latitudes (60–75^∘ N) is also captured by the model (Fig. 6a). The ⁷Be_P16spa shows high concentrations in the polar stratosphere and low values over the equatorial stratosphere (Fig. S6), mainly reflecting the latitudinal distribution of the production. This “latitudinal structure” is modulated for ¹⁰Be_P16spa in the stratosphere as ¹⁰Be is better mixed than ⁷Be due to its slow decay, together with a relatively long residence time in the stratosphere (Waugh and Hall, 2002). Both ⁷Be and ¹⁰Be show very low concentrations in the tropical upper troposphere, reflecting the frequent injection of air from the lower troposphere in wet convective updrafts, where aerosols are efficiently scavenged (Fig. S6).

The model also reasonably simulated ¹⁰Be vertical profiles compared with observations, with a tendency to underestimate observations in the stratosphere (Fig. 6j–l). A previous general circulation model study by Heikkilä et al. (2008b) also showed model stratospheric ¹⁰Be that is too low compared to measurements. They attributed this underestimation to too short a stratospheric air residence time in the model, which prevents ¹⁰Be concentrations from sufficiently accumulating in the stratosphere. However, this may not be the case in our study, as the stratospheric air residence time in the MERRA-2 reanalysis agrees reasonably with the observations (Chabrillat et al., 2018). Another explanation is that the ¹⁰Be production rate may be underestimated in the stratosphere. ⁷Be is less affected by this process than ¹⁰Be because of its short half-life compared to its stratospheric residence time (Delaygue et al., 2015).

Figure 7 Seasonal cycle of simulated and measured surface air ⁷Be concentrations, MERRA-2 total precipitation (4^∘ × 5^∘, bar graph), and modeled stratospheric contributions to surface air. The plots are arranged based on the site latitudes. The model results using the LP67 production rate are normalized to the ones using the P16spa production rate.

3.5 Global budgets and residence time

Table 1 shows the global budgets for ⁷Be_P16spa and ¹⁰Be_P16spa over the period of 2008–2018. About 22.1 % of tropospheric ⁷Be_P16spa is lost by radioactive decay, 75.8 % by convective and large-scale precipitation, and 2.1 % by dry deposition. The wet deposition contributes to about 97 % of total deposition for ⁷Be_P16spa and ¹⁰Be_P16spa (Table 1; Fig. S7), which is slightly higher than the ∼93 % contribution in previous model studies (Heikkilä et al., 2008b; Koch et al., 1996; Spiegl et al., 2022). The global mean tropospheric residence time of ⁷Be_P16spa is about 21 d, which is comparable to that reported by previous model studies: 18 d by Heikkilä et al. (2008b) and 21 d by Koch et al. (1996) and Liu et al. (2001). This also agrees with the residence time of about 22–35 d estimated from the observed deposition fluxes and air concentrations at 30–75^∘ N (Bleichrodt, 1978). The averaged tropospheric residence time of ¹⁰Be_P16spa is about 24 d, which is consistent with the 20 d suggested by Heikkilä et al. (2008b).

Table 1 Global budgets of ⁷Be and ¹⁰Be averaged from 2008 to 2018 in GEOS-Chem using P16spa.

* Against deposition only.

3.6 Seasonality in ⁷Be and ¹⁰Be

The seasonality of ⁷Be is influenced by (a) the amount of precipitation, (b) the stratosphere–troposphere exchange processes, and (c) the vertical transport of ⁷Be in the troposphere. The roles of these factors may vary depending on location. We compare the seasonal variations of modeled ⁷Be_P16spa and ⁷Be_LP67 concentrations with measurements from a dataset compiled by Terzi and Kalinowski (2017), with the data covering more than 6 years (Fig. 7). It should be noted that the model ⁷Be results and MERRA-2 precipitation rates are averaged over the years of 2008–2018, while the measurements are based on the data availability over the period 2001–2015.

In the Southern Hemisphere from 25–40^∘ S, the ⁷Be concentration peak is observed in austral summer (December–February), resulting from the combined influence of stratospheric intrusions and strong vertical transport during this season (Villarreal et al., 2022; Zheng et al., 2021b; Koch et al., 1996). The summer peak is also observed at northern mid-latitudes. This “summer peak” feature is well simulated by the model at some sites (e.g., KWP40 (29.3^∘ N, 47.9^∘ E), AUP04 (37.7^∘ S, 145.1^∘ E) and AUP10 (31.9^∘ S, 116^∘ E) shown in Fig. 7) but not at others (e.g., GBP68 (37.1^∘ S, 12.3^∘ W) and PTP53 (37.7^∘ N, 25.7^∘ W) in Fig. 7). This may not be related to stratospheric intrusion in the model as the simulated stratospheric contributions (Fig. S4) agree fairly well with estimates inferred from measurements, i.e., ∼25 % on annual average at northern mid-latitude surface (Dutkiewicz and Husain, 1985; Liu et al., 2016). Hence this could be due to the errors in vertical transport (e.g., convection) during the summer season.

The sites at northern high latitudes (>50^∘ N) show spring peaks that are well simulated by the model (e.g., ISP3 (64.1^∘ N, 21.9^∘ W)). This spring peak coincides with high stratospheric contributions, reflecting the influence of stratospheric intrusions. The influence of precipitation changes is also seen at several sites, especially in locations with high precipitation rates (e.g., monsoon regions). For example, two sites from Japan (JPP38 (36.3^∘ N, 139.1^∘ E) and JPP37 (26.5^∘ N, 127.9^∘ E) in Fig. 7) show summer minima coinciding with the high precipitation, even with relatively high stratospheric contributions in the same month.

The seasonal variation of stratospheric contribution is quite similar for the sites located in the Northern Hemisphere, with a high contribution in spring and a low contribution in fall. This is consistent with the estimates based on air samples that indicate stratospheric contributions varying from ∼40 % in spring to ∼15 % in fall at latitudes 38–51^∘ N (Dutkiewicz and Husain, 1985).

Generally, the model simulates the annual cycle of surface air ⁷Be concentrations well for most sites in terms of amplitude and seasonality (Fig. 7). For a few sites (e.g., DEP33 (47.9^∘ N, 7.9^∘ E)), the model captures the observed seasonality but not the correct absolute values. This could be partly due to the coarse resolution of the model. The ⁷Be_LP67 is normalized to ⁷Be_P16spa as we focus on the comparison of seasonal variability between these simulations. The very similar features (differences within 1 %) between all simulations using different production rates indicate a dominant influence of the meteorological conditions on the seasonal variations of the air ⁷Be concentrations.

Figure 8 Seasonal cycle of simulated (color lines) and measured (black line) ⁷Be deposition fluxes together with MERRA-2 total precipitation (4^∘ × 5^∘, bar graph). The model results using the LP67 production rate are normalized to the ones using the P16spa production rate.

Figure 8 compares model results with the seasonal ⁷Be deposition flux observations over the overlapping periods. Usually, high precipitation leads to high ⁷Be deposition fluxes (e.g., Du et al., 2015). Interestingly, low deposition fluxes are observed during the summer season in Taipei (Lee et al., 2015; Huh et al., 2006) coinciding with high precipitation. This feature is well captured in the model. Taipei has a typhoon season in summer when strong precipitation can occur in a very short period. The atmospheric ⁷Be could be removed quickly at the early stage of the precipitation event, while at the later stage there is little ⁷Be left in the air that can be removed (Ioannidou and Papastefanou, 2006).

Figure 9 Seasonal cycle of simulated ¹⁰Be deposition fluxes (2008–2018) and measured ¹⁰Be deposition fluxes in GRIP (1986–1990) and DSS (2000–2009) ice cores. The solid lines (grey) refer to seasonal variations of the measurements for each year. The solid black line indicates seasonal data of measurements in the year 1988. The dashed lines indicate the averaged seasonal variations of measured ¹⁰Be (black), ¹⁰Be_P16spa (blue), and ¹⁰Be_LP67 (red) concentrations.

To examine the ability of model to simulate ¹⁰Be in polar regions, we compare model results with two sub-annual ice core records (Fig. 9): the GRIP record from Greenland (1986–1990) (Heikkilä et al., 2008c) and the DSS record from Antarctica (2000–2009) (Pedro et al., 2011b). It should be noted that the direct measurements from ice cores are concentrations in the ice (atoms g⁻¹). To calculate deposition fluxes, the ice concentrations are multiplied by ice accumulation rates. However, for sub-annual accumulations, this bears large uncertainties. Therefore, we calculate the modeled ¹⁰Be concentrations for the selected sites using the model deposition fluxes at the selected sites, timed by ice density and then divided by the corresponding model precipitation rates.

Firstly, there is no consistent seasonal cycle in the GRIP ¹⁰Be measurement, indicating a strong role of local meteorology. The model does not reproduce the mean seasonal cycle, partly because the model was not run for the exact same period. However, we note that the measurements for the year 1988 show an annual cycle similar to that in the model, suggesting that the model ¹⁰Be seasonality falls within the range of the observations. For the DSS site, the model simulates the austral winter minima but not the austral fall maxima (February–April). These model biases could be due to the limited model resolution and local effects (e.g., ice redistribution due to wind blow) that are not resolved by the model. Such discrepancies were also reported by previous model studies using the ECHAM5-HAM general circulation model (2.8^∘ × 2.8^∘) over the overlap period (Heikkilä et al., 2008c; Pedro et al., 2011a). Global model simulations at higher resolutions or using a regional model could help improve the agreements between model results and measurements in Greenland and Antarctica. However, it should be kept in mind that local surface processes can cause a high degree of spatial variability in the impurity concentrations in ice cores, even for short distances (Gfeller et al., 2014), which cannot be resolved in climate models.

Figure 10 (a, b) Simulated ¹⁰Be${}_{\mathrm{P}\mathrm{16}\mathrm{spa}}{/}^{\mathrm{7}}$Be_P16spa ratio in spring (March–May) (a) and autumn (September–November) (b), averaged over the years 2008–2018. (c) Comparison between the annual averaged model ¹⁰Be${}_{\mathrm{P}\mathrm{16}\mathrm{spa}}{/}^{\mathrm{7}}$Be_P16spa ratios (lines) and those from measurements (circles; Jordan et al., 2003). The comparison is shown for the latitude bands of 60–75^∘ N and 45–60^∘ N, respectively.

3.7 ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{7}}\mathrm{Be}$ ratio

Figure 10 shows the modeled zonal mean ${}^{\mathrm{10}}{\mathrm{Be}}_{\mathrm{P}\mathrm{16}\mathrm{spa}}{/}^{\mathrm{7}}{\mathrm{Be}}_{\mathrm{P}\mathrm{16}\mathrm{spa}}$ ratios during boreal spring (March–May) and austral spring (September–November), respectively, when the stratosphere–troposphere exchange is strong in either of the two hemispheres. Also shown is the comparison of the altitudinal profile of the ${}^{\mathrm{10}}{\mathrm{Be}}_{\mathrm{P}\mathrm{16}\mathrm{spa}}{/}^{\mathrm{7}}{\mathrm{Be}}_{\mathrm{P}\mathrm{16}\mathrm{spa}}$ ratio with measurements from three aircraft missions (Jordan et al., 2003). The model ${}^{\mathrm{10}}{\mathrm{Be}}_{\mathrm{P}\mathrm{16}\mathrm{spa}}{/}^{\mathrm{7}}{\mathrm{Be}}_{\mathrm{P}\mathrm{16}\mathrm{spa}}$ ratio generally lies within the ranges of measurements (Fig. 10c). Due to the decay of ⁷Be and long residence time in the stratosphere, the ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{7}}\mathrm{Be}$ ratio is higher (>1.5) in the stratosphere and increases with altitude, with a maximum (>10) in the tropical stratosphere. During the period without strong stratospheric intrusion (e.g., the autumn season in Northern Hemisphere, Fig. 10b), the monthly ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{7}}\mathrm{Be}$ ratio near the surface is around 0.9–1. This surface ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{7}}\mathrm{Be}$ ratio could be up to 1.4 when the strong stratosphere–troposphere exchange happens (e.g., the spring season in Northern Hemisphere; Fig. 10a).

Figure 11 Comparison of monthly mean ⁷Be (a) and ¹⁰Be (b) concentrations and the ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{7}}\mathrm{Be}$ ratio (c) between model results with P16spa production and measurements for the Tokyo station over the period 2008–2014. Note that all ⁷Be and ¹⁰Be values are normalized to focus on variability. The dashed black line bridges the gap in measurements.

Figure 11 compares model surface air ⁷Be_P16spa and ¹⁰Be_P16spa concentrations and ${}^{\mathrm{10}}{\mathrm{Be}}_{\mathrm{P}\mathrm{16}\mathrm{spa}}{/}^{\mathrm{7}}{\mathrm{Be}}_{\mathrm{P}\mathrm{16}\mathrm{spa}}$ ratios with monthly mean observations in Tokyo (Yamagata et al., 2019) during the period of 2008–2014. Here we mainly focus on the relative variations, and ⁷Be and ¹⁰Be data are normalized. The model captures the observed variability in Tokyo well. ⁷Be and ¹⁰Be show a peak in early spring (March–May), while the ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{7}}\mathrm{Be}$ ratio shows a wider peak over March–July. The summer minima of ⁷Be and ¹⁰Be are due to strong scavenging associated with the monsoon/typhoon season precipitation. While the ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{7}}\mathrm{Be}$ ratio is independent of precipitation scavenging, the peaks of ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{7}}\mathrm{Be}$ coincide well with the enhancements of stratospheric contribution in the model. This indicates that the ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{7}}\mathrm{Be}$ ratio is a better indicator of the vertical transport and stratospheric intrusion influences than either tracer alone.

Figure 12 Comparison of annual mean model surface air ⁷Be concentrations with measurements from 2008–2018. Also shown is the model tropospheric ⁷Be production (green lines) at each station. All data are normalized by being divided by the mean over the first 5 years. The linear Spearman correlation coefficient R value is between ⁷Be_P16spa and measurements, while the values in brackets are between ⁷Be_LP67 and measurements.

3.8 Solar modulation influences

Here we examine the ability of the model to simulate the inter-annual variability of ⁷Be surface air concentrations, especially whether the model can simulate the solar modulation influence using the updated production model. Figure 12 shows the comparison of model-simulated annual mean surface air ⁷Be concentrations with measurements during 2008–2018 from four sites: Kiruna, Ljungbyhed, Vienna, and Hong Kong SAR (Kong et al., 2022; Zheng et al., 2021b). The tropospheric ⁷Be production rate from each site is also plotted for comparison as measured annual mean surface air ⁷Be concentrations are predominantly influenced by the local tropospheric ⁷Be production signal (Zheng et al., 2021b).

The model ⁷Be_P16spa surface air concentrations show a better agreement with annual ⁷Be measurements (higher R value) compared to ⁷Be_LP67 concentrations at all surface sites (Fig. 12). The variability in the measurements (Kiruna, Ljungbyhed, and Vienna) agrees well with the trend in production, suggesting a dominant influence of solar modulations during this period. This is further supported by strong deviations between ⁷Be_P16spa and ⁷Be_LP67 as no solar influence is considered in ⁷Be_LP67. This also emphasizes the importance of including solar modulation of the ⁷Be and ¹⁰Be production in modeling studies, especially for high-latitude regions. The mismatch of measurements and production at Kiruna from 2012 to 2015, together with the similar year-to-year variability between ⁷Be_P16spa and ⁷Be_LP67, suggests the meteorological influence is dominant at Kiruna for this period. This also suggests that meteorological influences can suppress the solar signal in the ⁷Be and ¹⁰Be observations.

4 Summary and conclusions

We have incorporated the ⁷Be and ¹⁰Be production rates derived from the CRAC:Be model, considering realistic spatial geomagnetic cutoff rigidities (P16spa), into the GEOS-Chem global chemical transport model, enabling the model output to be quantitatively comparable with the measurements. In addition to the standard simulation using the P16spa production rate, we further conducted two sensitivity simulations: one with the default production rate in GEOS-Chem based on an empirical approach (LP67) and one with the production rate from the CRAC:Be but considering only geomagnetic cutoff rigidities for a geocentric axial dipole (P16). On a global average, the LP67 production rate is 60 % higher compared to that of P16 and P16spa. The P16 production rate shows some regional differences (up to 50 %) compared to the P16spa production rate.

In comparison with a large number of air and deposition flux measurements, the model ⁷Be_P16spa shows good agreements with respect to surface air concentrations (93.7 % of data within a factor of 2) and reasonably good agreements regarding deposition fluxes (60.9 % of data within a factor of 2). The model simulates the surface air concentration peaks well in the subtropics, associated with strong downward transport from the stratosphere. This agreement is better than that using the default ⁷Be_LP16 production and the ⁷Be_P16 production with simplified axis symmetric dipole cutoff rigidity. The ⁷Be_LP67 simulation overestimates the absolute value of ⁷Be. The ⁷Be_P16 simulation tends to produce a positive bias (∼18 %) for the ⁷Be deposition fluxes in the East Asia region; nevertheless, no large bias is found for ⁷Be surface air concentrations. The surface deposition fluxes are more sensitive to the production in the mid-troposphere and upper troposphere and downward transport of ⁷Be from the stratosphere, due to the effect of precipitation scavenging throughout the troposphere.

For the first time, the ability of GEOS-Chem to simulate ¹⁰Be is assessed with measurements. The ¹⁰Be_P16spa model results agree well with ¹⁰Be observational data that were evaluated for dust influences or from the regions less influenced by dust (e.g., polar regions), while underestimating most samples that were not corrected for dust influences. This highlights the importance of examining the dust contribution to ¹⁰Be measurements when using these data to evaluate models.

Independent of the production models, surface ⁷Be and ¹⁰Be concentrations from all three simulations show similar seasonal variations, suggesting a dominant meteorological influence. The model generally simulates the annual cycle of ⁷Be surface air concentrations and deposition fluxes well at most sites in terms of amplitude and seasonality. The model fails to capture the “summer peak” in a few sites likely due to errors in convective transport during summer.

The model ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{7}}\mathrm{Be}$ ratios also lie within the measurements, suggesting the stratosphere–troposphere exchange process is reasonably represented in the model. The mismatch of the peaks between ⁷Be(¹⁰Be) and ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{7}}\mathrm{Be}$ ratios at the Tokyo site suggests that the ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{7}}\mathrm{Be}$ ratio is a better indicator of the vertical transport and stratospheric influences than either tracer alone as the ratio is independent of precipitation scavenging.

Finally, we demonstrate the value and importance of including time-varying solar modulation in ⁷Be and ¹⁰Be production rates for model simulations of both tracers. It significantly improves the agreement of interannual variations between the model and measurements, especially at surface sites from mid-latitudes and high latitudes. The mismatch of trends in the modeled ⁷Be production rate and observed air concentrations at Kiruna from 2012–2015 also suggests that the solar signal can be suppressed by meteorological influences.

In summary, we have shown that with the state-of-the-art P16spa production rate, the ability of GEOS-Chem to reproduce the ⁷Be and ¹⁰Be measurements (including interannual variability of ⁷Be) is significantly improved. While uncertainties in transport and deposition processes play a major role in the model performance, reduced uncertainties in the production rates, as demonstrated in this study, allow us to use ⁷Be and ¹⁰Be tracers as better tools for evaluating and testing transport and scavenging in global models. We recommend using the P16spa (versus default LP67) production rate for GEOS-Chem simulations of ⁷Be and ¹⁰Be in the future.