Numerical coupling of aerosol emissions, dry removal, and turbulent mixing in the E3SM Atmosphere Model version 1 (EAMv1), part I: dust budget analyses and the impacts of a revised coupling scheme

An earlier study evaluating the dust life cycle in the Energy Exascale Earth System Model (E3SM) Atmosphere Model version 1 (EAMv1) has revealed that the simulated global mean dust lifetime is substantially shorter when higher vertical resolution is used, primarily due to significant strengthening of dust dry removal in source regions. This paper demonstrates that the sequential splitting of aerosol emissions, dry removal, and turbulent mixing in the model's time integration loop, especially the calculation of dry removal after surface emissions and before turbulent mixing, is the primary reason for the vertical resolution sensitivity reported in that earlier study. Based on this reasoning, we propose a simple revision to the numerical process coupling scheme, which moves the application of the surface emissions to after dry removal and before turbulent mixing. The revised scheme allows newly emitted particles to be transported aloft by turbulence before being removed from the atmosphere, and hence better resembles the dust life cycle in the real world. Sensitivity experiments are conducted and analyzed to evaluate the impact of the revised coupling on the simulated aerosol climatology in EAMv1.


Introduction
The numerical modeling of aerosol-climate interactions is a research topic with high levels of uncertainty (Seinfeld et al., 2016;Bellouin et al., 2020;Smith et al., 2020).State-of-the-art global aerosol-climate models (e.g., Stier et al., 2005;Mann et al., 2010;Liu et al., 2012;Zhang et al., 2012;Wang et al., 2020) often consider a fairly large number of physical and chemical processes that affect aerosol life cycles.Such processes include emissions, new particle formation, particle growth and aging, transport by winds at local to global scales, as well as removal from the atmosphere due to gravitational settling, boundary layer processes, formation of clouds and precipitation, and wet removal by precipitation.Typically, the numerical representations of different processes are separately developed by subject matter experts and subsequently assembled to form global models.
Different aerosol-climate models are known to show substantial discrepancies in the simulated magnitudes and spatiotemporal variations of the source and sink processes (Textor et al., 2006;Gliß et al., 2021).Continuous efforts are being made to attribute such discrepancies to the physical assumptions and computational methods used for individual processes, while investigations on the impact of numerical process coupling (e.g., Wan et al., 2013) are relatively rare.
The Energy Exascale Earth System Model ( E3SM) is an Earth system model developed by the U.S. Department of Energy for addressing science questions related to the prediction of Earth system dynamics and climate change (Leung et al., 2020;E3SM Project, 2018).The E3SM Atmosphere Model version 1 (EAMv1, Rasch et al., 2019) is an atmospheric general circulation model that includes a comprehensive representation of the life cycles of various natural and anthropogenic aerosol species.A recent study by Feng et al. (2022) evaluated the dust life cycle in EAMv1 simulated with different horizontal and vertical resolutions.It was shown that an EAMv1 simulation using 1 • horizontal grid spacing and 72 vertical layers produced a global mean dust lifetime of 1.85 days, which was substantially shorter than the lifetime of 2.6 days reported by Liu et al. (2012) and Scanza et al. (2015) who used the predecessor model CAM5, the Community Atmosphere Model version 5 (Neale et al., 2012), with 2 • horizontal grid spacing and 30 layers.Feng et al. (2022) also showed that reverting EAMv1's vertical resolution from 72 layers to 30 layers would lengthen the global mean dust lifetime from 1.85 days to 2.4 days (see Table 4 therein), primarily through the weakening of dry removal in the dust source regions.This paper investigates the vertical resolution sensitivities of dust lifetime and dry removal reported in Feng et al. (2022).
Sections 2.1 and 2.2 provide a brief overview of EAMv1 and its parameterizations of aerosol emissions, dry removal, and turbulent mixing.Section 2.3 describes EAMv1's numerical scheme used for coupling the aerosol-related processes.Section 3 analyzes the simulated dust mass budget and compares simulations conducted with 72 or 30 layers to reveal weaknesses of the numerical process coupling in the publicly released EAMv1 (which we refer to as the original or default EAMv1 in this paper).A simple revision to the numerical process coupling is proposed in Sect.3.3 and its impacts on the simulated aerosol climatology are evaluated in Sect. 4. The conclusions are drawn in Sect. 5.
As is shown in Sect.4, the revised process coupling significantly weakens the dust dry removal in the source regions and substantially reduces sensitivities of the simulated dust dry removal rate and lifetime to vertical resolution.These changes appear to be desirable given the deficiencies of the EAMv1 results pointed out in Feng et al. (2022).On the other hand, the revision also results in large, global and systematical increases in the mass burden of dust and sea salt, meaning the top-of-the-atmosphere energy fluxes become out of balance.While such aerosol burden and energy flux changes can be offset by retuning uncertain parameters used in the dust and sea salt emission parameterizations, one might ask whether the revised coupling provides better results (in terms of dry removal and dust lifetime) for the right reasons and whether the retuning is worthwhile.
Other than presenting in Sect. 3 some process-level analyses of EAMv1 results to motivate the revision and confirming in Sect. 4 that the effects of the revised coupling on EAM's aerosol climatology meet our expectation, it is not straightforward to obtain additional evidence from EAM simulations to show that the revised coupling is an improvement in the numerical sense.
This challenge has to do with the fact that the current EAM code does not allow for convergence testing of numerical process coupling between aerosol emissions, dry removal, and turbulent mixing without changing the timestep sizes and coupling frequencies of other physical processes in the model, for example resolved fluid dynamics and parameterized cloud processes.
To address this challenge, the companion paper by Vogl et al. (2023) presents a proof from an applied mathematics perspective that the numerical error associated with dust dry removal caused by process splitting at each timestep is smaller when the revised scheme is used, hence further justifying the adoption of the revised scheme in EAMv1.

The EAMv1 model and simulations
The EAMv1 configuration described in Rasch et al. (2019) and used in this study is a global hydrostatic atmospheric model that simulates the spatial distribution and time evolution of air temperature, pressure, winds, humidity, clouds, and precipitation.
In addition, the model has 50 prognostic variables corresponding to the mixing ratios of aerosol particles of different sizes (diameters), chemical compositions, and attachment states (interstitial or cloud-borne).The mixing ratios of a few chemical gas species that are precursors of aerosols are also simulated using prognostic equations.

EAMv1 overview
EAMv1's dynamical core solves the primitive equations of the global atmospheric flow using a continuous Galerkin spectralelement method on a cubed-sphere horizontal mesh (Dennis et al., 2012;Taylor et al., 2009).The vertical discretization uses a semi-Lagrangian scheme and a pressure-based terrain-following coordinate (Lin, 2004).The resolved-scale tracer transport uses the discretization method of Lin (2004) but has been adapted to the cubed sphere.The transport algorithm ensures local conservation of tracer and air masses as well as moist total energy (Taylor, 2011).
For the parameterization of unresolved processes, the transfer of solar and terrestrial radiation is calculated with the Rapid Radiative Transfer Model for General Circulation Models (RRTMG, Iacono et al., 2008;Mlawer et al., 1997).Deep convection is parameterized with the scheme of Zhang and McFarlane (1995) with modifications by Neale et al. (2008) and Richter and Rasch (2008).Shallow convection, turbulence, and stratiform cloud macrophysics are represented by the higher-order closure parameterization named Cloud Layers Unified By Binormals (Larson, 2017;Larson and Golaz, 2005;Golaz et al., 2002;Larson et al., 2002).Stratiform cloud microphysics is represented with a two-moment parameterization with prognostic equations for the mass and number concentrations of cloud droplets, ice crystals, rain, and snow (Gettelman and Morrison, 2015;Gettelman et al., 2015;Morrison and Gettelman, 2008).More detailed descriptions of the parameterization suite can be found in Section 2 of Rasch et al. (2019) and the references therein.

Aerosol processes in EAMv1
The Modal Aerosol Module (MAM, Wang et al., 2020;Liu et al., 2016Liu et al., , 2012;;Ghan and Easter, 2006;Easter et al., 2004) is a suite of parameterizations developed for global climate modeling, aiming at providing a simplified yet sophisticated representation of aerosol life cycles as well as their interactions with clouds and precipitation.Seven aerosol species are considered in EAMv1: sulfate, black carbon (BC), primary organic aerosols (POA), secondary organic aerosols (SOA), marine organic aerosols (MOA), dust, and sea salt.Mass and number concentration changes are predicted for aerosol particles of two attachment states, interstitial and cloud-borne, which refer to the particles found outside and within cloud droplets, respectively.The two attachment states correspond to two populations of aerosol particles.In the 4-mode configuration of MAM used in EAMv1 (i.e., MAM4), the particle size distribution in each population is represented by one coarse mode and three fine-particle modes, see Fig. 2 in Wang et al. (2020).Within a mode, the particle size distribution is represented by a lognormal function of particle diameter, assuming the particles are spherical and the geometric standard deviation of the lognormal function is fixed.Under these assumptions, the mass mixing ratios of different species in different modes, as well as the particle number mixing ratios of the different modes, are predicted.
Aerosol processes currently considered in MAM4 include emissions as well as particle mass and number changes resulting from new particle formation (aerosol nucleation), condensation and evaporation of chemical species, water uptake, coagulation, aqueous chemistry in cloud droplets, aerosol activation (cloud nucleation), resuspension from evaporating cloud droplets and precipitation, in-cloud and below-cloud wet removal, sub-grid vertical transport by deep convection and turbulence, gravitational settling, and turbulent dry deposition.Descriptions of MAM4 and its predecessors can be found in Rasch et al. (2019), Wang et al. (2020), Liu et al. (2016), Liu et al. (2012), Ghan andEaster (2006), andEaster et al. (2004).Here we only briefly summarize the processes that are the foci of this study.

Emissions
Since MAM considers a variety of sizes, species, and origins of aerosol particles, a comprehensive set of assumptions and treatments are needed to specify aerosol emissions, see Section S1.1.1 in Liu et al. (2012).In MAM4, the prescription of anthropogenic aerosol mass emissions follows protocols of model intercomparision projects and published studies.The partitioning of the mass emissions into MAM's lognormal modes and the calculation of aerosol number emissions in those modes are based on assumed emission size distributions.Natural aerosol emissions are calculated online (i.e., during a simulation) using emission parameterizations.The scheme from Zender et al. (2003) is used for dust.The calculation of sea salt emissions follows Mårtensson et al. (2003) for particle diameters from 20 nm to 2.5 µm and Monahan et al. (1986) for particle diameters from 2.5 µm to 10 nm.For MOA, the parameterization of Burrows et al. (2014) is used to calculate the mass portion of MOA in the emitted sea spray aerosols; this information is used in combination with the predicted sea salt emissions to determine Table 1.Global annual mean sources of aerosol mass (unit: Tg yr −1 ) from surface emissions, elevated emissions, and production inside the atmosphere.Also presented are the relative contributions (unit: %) of these sources to the total source of the corresponding species.The "-" symbol means not applicable.The results were derived from a wind-nudged simulation of the year 2010 conducted with the default EAMv1 configuration using 1 • horizontal grid spacing and 72 layers (simulation v1_ori_L72, see Sect.2.4).

Mass source of species
Surface emissions Elevated emissions Atmospheric production MOA emissions, assuming that MOA emissions add to the sea spray aerosol emissions, and that the emitted MOA is internally mixed with other aerosol species (Burrows et al., 2022).
Dust aerosols are assumed to be emitted only at the Earth's surface.The emission fluxes are parameterized as a function of various properties of the Earth's surface (e.g., soil moisture and erodibility) and atmospheric conditions (e.g., friction velocity and 10 m wind speed).Further details can be found in Section 2.4 of Zhang et al. (2016) and the references therein.In mathematical terms, the emission of dust aerosols results in source terms in the evolution equations of dust mixing ratio that are non-zero only at the bottom boundary of an atmosphere column.While this is different from the treatment in the IMPACT global chemistry and transport model which injects dust and some other aerosols uniformly in the boundary layer (see, e.g., Liu et al., 2005;Wang et al., 2009), the assumption of non-zero dust emissions only at the Earth's surface is commonly used in many global aerosol-climate models (e.g., Gong et al., 2003;Stier et al., 2005;Mann et al., 2010;Zhang et al., 2010).
Table 1 shows the global annual mean magnitudes of the surface and elevated emissions of the 7 aerosol species considered in MAM4.To put the numbers into context, the mass sources caused by physical and chemical production inside the atmosphere are also shown, together with the percentage contribution of each source.In Table 2, the surface emissions of species that have no atmospheric production (i.e., only emissions as sources) are broken into contributions from the 4 lognormal modes.The numbers were derived from a wind-nudged simulation of the year 2010 conducted with the default EAMv1 configuration using 1 • horizontal grid spacing and 72 layers.The simulation setup can be found in Sect.2.4.

Turbulent dry deposition and gravitational settling
In this paper, we use the term "dry removal" to refer to both the gravitational settling and the turbulent dry deposition of aerosol particles, as the two processes are handled together in the MAM4 parameterization suite embedded in EAMv1.Gravitational settling is the downward movement of particles under the action of gravity, whereas turbulent dry deposition refers to the loss of particles to the Earth's surface through Brownian diffusion, impaction, interception, etc.For both processes, the downward aerosol mass or number fluxes per unit time across a unit area at the Earth's surface are calculated using the formula Here, M i,sfc is the downward flux of the i-th aerosol tracer, ρ b is air density in the bottom layer of the atmosphere (i.e., the lowest model layer above the Earth's surface), and q i,b is the mixing ratio of aerosol tracer i in the bottom layer.v i,b is the downward deposition velocity of aerosol tracer i calculated using the aerosol properties and ambient conditions in the bottom layer.
Like in many other models (see, e.g., Mann et al., 2010;Zhang et al., 2012), MAM4 in EAMv1 assumes gravitational settling of aerosols can occur throughout the atmosphere column.The resulting fluxes at altitudes above the Earth's surface are parameterized as where z is geopotential height; ρ(z), q i (z), v grav i (z), and M grav i (z) are air density, aerosol mixing ratio, gravitational settling velocity, and gravitational settling flux, respectively, at altitude z.The calculation of the settling velocity is based on the Stokes' law, assuming particles reach their terminal velocities instantly.The settling velocity of a single particle is calculated with Eqs. ( 2) and (3) in Zhang et al. (2001).Correction factors are included to account for the impact of the lognormal size distribution.
The turbulent dry deposition velocity is parameterized using Eq. ( 21) in Zender et al. (2003), for which the calculation of the quasi-laminar layer resistance follows Sect. 2 in Zhang et al. (2001).

Turbulent mixing
Turbulent mixing of aerosols in MAM4 is parameterized using the eddy diffusivity approach (see, e.g., Garratt, 1994), which gives aerosol concentration tendencies in the form of where ρ is air density, q i is mixing ratio of aerosol tracer i, z is geopotential height, and K h is eddy diffusion coefficient.
In EAMv1, K h is calculated by the turbulence and cloud parameterization CLUBB, while the turbulent mixing of aerosol mass, aerosol number, and cloud droplet number is treated separately (outside CLUBB) in conjunction with aerosol activation and resuspension from evaporating cloud droplets.In other words, MAM4's parameterization of turbulent mixing and aerosol activation-resuspension solves differential equations in the form of where the last term in Eq. ( 4) is the rate of change caused by aerosol activation and resuspension from cloud droplets.Additional information about the parameterization can be found in Section S1.1.8 in the supplementary materials of Liu et al. (2012) and in Ghan and Easter (2006).
For clarification, we note that MAM also accounts for the resuspension of aerosols from evaporating precipitation as a part of aerosol wet removal.In this paper, when resuspension is mentioned in conjunction with activation, we are referring to the resuspension from cloud droplets.

Numerical process coupling in EAMv1
The schematic in Fig. 1 depicts the sequence in which the various atmospheric processes are calculated within a timestep of 30 minutes in the 1 • EAMv1 simulations.The schematic also shows where the coupling between EAM and the other components of E3SM (e.g., land and ocean) occurs during the 30-minute timestep.EAMv1 uses primarily the sequential splitting method for the numerical coupling of aerosol-related processes and most other parameterizations.With this method, the rates of change (i.e., tendencies) of aerosol mixing ratios caused by a process (or a process group) are calculated and then applied immediately to obtain new (updated) values of the mixing ratios.The updated mixing ratios are then passed on to 1. Emission fluxes of aerosols are calculated based on atmospheric and Earth surface conditions or derived from emission datasets.Dust emissions are calculated in E3SM's land model and passed to EAM by the coupler.Emissions of the other species are calculated (or read in from data files, partitioned to MAM's lognormal modes, and mapped to EAM's vertical grid) within the atmosphere model.2. Gas-phase chemistry, aqueous-phase chemistry, and aerosol microphysics parameterizations (gas-aerosol mass transfer, new particle formation, inter-mode transfer due to particle growth, coagulation, and aging of primary carbon particles to accumulation mode) are calculated for a 30-minute time window.Elevated emissions of aerosols and precursor gases are also applied, as well as wet removal of the gases.Rates of conversion between aerosol and gas species and between different aerosol modes are calculated, and the corresponding mixing ratio tendencies are used to update the aerosol and gas mixing ratios.
3. The surface fluxes of tracers (not including water vapor) stored in the Fortran variable cam_in%cflx are converted to mixing ratio tendencies in the lowest model layer.These tendencies are applied over a 30-minute timestep to update the corresponding mixing ratios.We note that the surface fluxes applied at this location of the time loop are the emission fluxes of aerosols and the net fluxes (emission minus turbulent dry deposition) of precursor gases.
4. Dry removal of aerosols is calculated for a 30-minute time window, and the aerosol mixing ratios are updated.Gravitational settling affects all model layers where aerosols are present, whereas turbulent dry deposition affects only the lowest model layer.
5. The large-scale transport scheme updates mixing ratios over two 15-minute vertical remapping timesteps, each of which consists of three 5-minute sub-cycles of horizontal advection.
6.The deep convection parameterization changes EAM's temperature, humidity, and wind profiles as well as mixing ratios of hydrometeors but not yet the aerosol mixing ratios.
7. The parameterization of turbulent mixing and activation-resuspension of aerosol particles is calculated.The cloud-borne and interstitial aerosol mass and number mixing ratios are updated within this parameterization.The tendencies of cloud droplet number mixing ratio are passed on to the stratiform cloud microphysics parameterization.In default 1 • EAM simulations, the turbulent mixing and activation-resuspension of aerosol particles, ice nucleation, the stratiform cloud microphysics, and the turbulence and shallow convection parameterization CLUBB are sub-cycled together using 5minute timesteps (see Sect. 2.1 in Wan et al., 2021 andSect. 2 in Santos et al., 2021).Within each sub-cycle, CLUBB handles the turbulent transport of heat, water, and precursor gases; CLUBB also provides the eddy diffusion coefficient, turbulent kinetic energy, and cloud fraction to the aerosol mixing and activation-resuspension parameterization.The equation of turbulent mixing of aerosols is solved not by CLUBB but in conjunction with aerosol activation-resuspension, using an explicit time-stepping method with dynamically determined stepsizes.
8. After the turbulence and cloud microphysics sub-cycles, the processes of aerosol water uptake, aerosol in-cloud and below-cloud wet removal, and the vertical transport of aerosols by deep convection are calculated for a timestep of 30 minutes, and the corresponding mixing ratios are updated.9.The mixing ratios, together with many other physical quantities describing the simulated atmospheric state, are recorded and passed to the software infrastructure that handles model output.In other words, the mixing ratios in EAM's output files that the model developers and users typically analyze are from this location in the time loop.(The mixing ratios at other locations presented later in the paper were obtained using the online diagnostic tool CondiDiag (Wan et al., 2022); those additional mixing ratio values were included in the output files under different variable names following the convention described in Wan et al. (2022).)10.The mixing ratios and other properties of aerosol particles are passed to the radiation parameterization where the aerosol impact on the atmospheric energy budget is calculated.Radiation does not directly affect aerosol mixing ratios.
After radiation and some additional diagnostics, the atmosphere model exchanges information with other components of E3SM, and the calculation goes back to step 1 listed above.The sequence of calculations described above is used by the original EAMv1.In the revised process coupling scheme discussed in Sects.3.3 and 4, the surface emissions of aerosols and the net surface fluxes (emissions minus turbulent dry deposition) of precursor gases are applied after deep convection and before the cloud macrophysics-microphysics sub-cycles.
That is, box 3 in Fig. 1 is moved to the location indicated by the dashed pink box (between boxes 6 and 7) when the revised coupling scheme is used.

Simulations
To analyze features of the simulated dust life cycle in the 1 • configuration of EAMv1, we present in later sections 15-month simulations started from October 2009 using the default timestep settings documented in Wan et al. (2021) and Santos et al. (2021).The results presented in this paper were derived from model output over the year 2010.
A total of 4 simulations were conducted, two of which used EAMv1's original process coupling scheme described in Sect.2.3, and the other two used the revised process coupling described in Sect.3.3.For each of the coupling schemes, we conducted one simulation with EAMv1's default L72 vertical grid and another simulation with the L30 grid following the earlier studies of Feng et al. (2022) and Zhang et al. (2018).The simulation short names are listed in Table 3.A comparison between the near-surface layers in the L30 and L72 grids is shown in Fig. 2.
The horizontal winds were nudged to 6-hourly ERA-Interim reanalysis (Dee et al., 2011) with a relaxation timescale of 6 hours following Sun et al. (2019).In this study, wind nudging was applied to model levels above 850 hPa approximately, (i.e., above level 58 in the L72 grid and level 24 on the L30), as previous studies have shown that nudging winds near the Earth surface can result in lower wind speeds compared to the free-running simulations, and consequently affect the simulated dust emissions (Timmreck and Schulz, 2004;Sun et al., 2019).The external forcings, e.g., the sea surface temperature and sea ice extent as well as the emissions of anthropogenic aerosols and their precursors, followed the CMIP6 protocol of the twentieth-century transient simulations.During a simulation, the evolution of dust mixing ratios within 30-minute timesteps and the tendencies associated with the colored boxes in Fig. 1 were tracked with the online diagnostic tool of Wan et al. (2022).3 Motivation for a revised coupling scheme This section presents budget analyses for the total mass mixing ratio of interstitial dust particles (Sect.3.1) and discusses the weaknesses of the original coupling scheme (Sect.3.2).A revised coupling scheme is described in Sect.3.3 and evaluated in the next section.

Dust mass budget
Geographical distributions of the annual mean dust mass emission fluxes and interstitial dust burden from the default EAMv1 configuration (i.e., simulation v1_ori_L72) are presented in panels (a) and (b) of Fig. 3.The mass burden shown here is the total mass mixing ratio summed over MAM4's accumulation mode and coarse mode, multiplied by air density, and vertically inte-grated over the atmosphere column.As expected, the dust emissions are highly inhomogeneous in space while the distribution of the dust burden is much smoother due to transport by winds.
To get an overview of the balance between different physical processes, we selected the magenta box in Fig. 3a to cover the major dust sources in Asia, Europe, and North Africa.Within the magenta box, the grid columns with non-zero annual mean emissions were marked as the source regions and the grid columns with no emissions were marked as the source-vicinity regions.The brown box in Fig. 3a was selected to represent regions far away from dust emissions (i.e., the remote regions).
The annual mean dust mass mixing ratio tendency profiles were averaged over the source, vicinity, and remote regions and are shown in Fig. 3c-h.
The characteristic tendency magnitudes in the three types of regions as well as at different altitudes can be seen in Fig. 3c-h.
The upper row is the results in the lower troposphere with the lowest model layer excluded, while the lower row shows the 10 model layers closest to the Earth surface, including the bottom layer.Furthermore, the tendencies in panels c and e-h have been multiplied by factors of 100 to 100,000 in order to be plotted with the same x-axis range as in panel d.These factors of multiplication are noted at the top of each panel.A key feature revealed by these vertical profiles is that the annual and regional mean tendencies in the lowest model layer in the dust source regions are 2 orders of magnitude stronger than the tendencies in the upper layers within the same regions.Furthermore, the tendencies in the lowest model layer in the source regions are 2-4 orders stronger than the tendencies in any layer in the source-vicinity or remote regions.The dominant source and sink terms in the global-scale annual mean dust mass budget are the following processes in the lowest model layer in the dust source regions: (1) surface emissions, (2) dry removal, and (3) turbulent mixing and aerosol activation-resuspension.
To demonstrate that similar contrasts in the magnitudes of process rates can be seen in individual grid columns and at individual timesteps, Fig. 4 and Fig. 5 present results derived from 6-hourly instantaneous model output in the 90-day period of 2010-01-29 to 2010-04-28. Figure 4 is a histogram of the dust mass emission fluxes in the magenta box in Fig. 3a, derived from grid columns and timesteps with non-zero dust emissions.The instantaneous emission fluxes span more than 10 orders of magnitude.About 82% of the emitted mass is attributable to events with mass fluxes between 2×10 −8 kg m −2 s −1 and 2×10 −6 kg m −2 s −1 ; about 13% of the emitted mass is contributed by a large number of weak emission events and the remaining 5% by a very small number of very strong events (Fig. 4).
For the middle portion of the histogram, Fig. 5 shows the statistical distributions of the instantaneous dust mixing ratio tendencies in the lowest 10 model layers of the L72 grid, where the different panels correspond to different physical processes.
The results for the other two portions of the histogram as well as for the source-vicinity and remote regions are shown in

Weaknesses of the original coupling scheme
Since sequential splitting is used in the default EAMv1 for the various aerosol processes discussed above, and the ordering is emissions, dry removal, resolved transport, turbulent-mixing-activation-resuspension, and wet removal, we show in Fig. 6 the statistical distributions of instantaneous dust mixing ratios after each of these processes has been calculated and the tendencies have been applied (i.e., the mixing ratios have been updated).The results in this figure (Fig. 6) correspond to the same grid columns and timesteps as shown in Fig. 5. Given the nature of the sequential splitting method and the magnitudes of tendencies seen in Fig. 5, it is understandable that large mixing ratio spikes are seen in the lowest model layer after the surface emissions are applied (Fig. 6a).The mixing ratio profiles in Fig. 6 also indicate that although dry removal is a strong sink in the lowest model layer and resolved transport can be a significant sink, neither is sufficiently strong to remove the near-surface mixing ratio spikes.In contrast, the turbulent mixing and aerosol activation-resuspension processes are very effective in vertically smoothing the profiles.Since the dust source regions are typically dry and the lowest model layer is thin (about 20 m on average), we do not expect aerosol activation to occur there frequently; hence the smoothing effect is most likely attributable to turbulent mixing.Using the online diagnostic tool of Wan et al. (2022), we also analyzed the dust mixing ratios after each of the 5-minute sub-cycles used for the parameterizations of turbulence mixing and aerosol activation-resuspension.We found that the spikes of dust mixing ratio were typically eliminated after 1 or 2 sub-cycles (i.e., 5 to 10 minutes).This is consistent with the results from Wang et al. (2011), who showed that particles injected from ships traveling below marine stratocumulus can be lofted and mixed through the cloud-topped boundary layer within minutes (see Fig. 1c and Sect.3.1 in Wang et al., 2011).
While the sub-timestep evolution of dust mixing ratio shown in Fig. 6 is to be expected from the use of sequential splitting and the characteristic magnitudes of process rates in the default EAMv1, the figure provides a clue that the ordering of the various process has severe deficiencies.In the real world, dust emissions typically occur in turbulent atmospheric environments where the turbulent eddies are efficient in transporting the emitted particles aloft.In EAM, when the emission fluxes are applied to update the dust mixing ratios in the lowest model layer, the effect can be interpreted as immediately mixing the emitted particles throughout the lowest model layer but also temporarily trapping the particles within that layer, resulting in high concentrations being passed to the next process in the time loop, which is dry removal in the default EAMv1.Since the dry removal fluxes at the Earth's surface are proportional to the mixing ratios in the lowest model layer (see Sect. 2.2.2, Eq. 2), the high mixing ratios depicted in Fig. 6a are expected to lead to strong dry removal fluxes.If one chose, instead, to calculate turbulent mixing immediately after the emissions are applied, and then calculate dry removal afterwards (i.e., after turbulent mixing), then the input to the dry removal parameterization would look similar to Fig. 6d (i.e., without a spike in the bottom layer), hence resulting in mush weaker dry removal fluxes at the surface.Following this logic, we expect the EAM simulations to be sensitive to the ordering of the emission, turbulent mixing, and dry removal processes.
The reasoning above also provides an explanation for the strong vertical resolution sensitivity of dust dry removal reported in Feng et al. (2022).As mentioned above, dust is emitted only into the lowest model layer in EAMv1, meaning the particles are temporarily trapped below the upper interface of the layer.The lowest layer in the L72 grid is about 1/5 in thickness compared to the lowest layer in the L30 grid (Fig. 2).Therefore, given the same emission fluxes, the temporary increases in dust mixing ratio after the emissions are applied in a simulation using the L72 grid are expected to be 5 times as high as in another simulation that uses the L30 grid (see schematic in Fig. 7), which in turn can lead to significantly stronger dry deposition in the L72 simulation.This expectation is confirmed by Fig. 8 which shows results from simulation v1_ori_L30.In the L30 simulation, both the tendencies and the mixing ratios spikes in the lowest model layer are substantially weaker than in v1_ori_L72.
It is worth clarifying that the process coupling issue identified here is not specific to EAMv1.In the predecessor model CAM5, although the atmospheric turbulence parameterization used the scheme of Park and Bretherton (2009) and was calculated before aerosol dry removal and after surface emissions, the aerosol tracers were not mixed by the Park and Bretherton (2009) parameterization, but rather by the same turbulent mixing and aerosol activation-resuspension parameterization as in EAMv1.In other words, as far as the aerosol tracers are concerned, the sequence of calculation in CAM5 was the same as in EAMv1.We expect that if CAM5 simulations are performed with EAM's L72 grid, significantly stronger dry removal and shorter dust lifetime will result as well.
Furthermore, we note that the same process coupling issue also exists for other aerosol species in EAM that have surface emissions, although the magnitude of the impact depends on the relative importance of the surface emissions as well as the typical sizes of the particles.This is further discussed in Sect.4.2 using EAM results.
The precursor gases in EAMv1, on the other hand, do not suffer from this coupling issue, because the splitting and ordering of the gas-related processes are different.Precursor gases in the model are assumed to experience turbulent dry deposition but no gravitational settling.The surface dry removal fluxes of gases are calculated after gas-phase chemistry inside the box marked with "2" in the schematic in Fig. 1 and then, instead of being used to update gas mixing ratios, the dry removal fluxes are subtracted from the surface emission fluxes, and the residual (i.e., the net flux) is saved in the Fortran variable cam_in%cflx.
When cam_in%cflx is used to update tracer mixing ratios in the lowest model layer in box 3 in the schematic, it is the residual of emission and dry deposition (i.e., the net flux) that is applied.Turbulent mixing of precursor gases is handled by CLUBB, i.e., inside box 7 in the schematic in Fig. 1.

Upper layers
Surface emission L30 L72

Surface emission
Figure 7.A schematic depicting the impact of layer thickness on the dust mixing ratio in the lowest model layer after surface emissions are applied.In a numerical model where the surface emissions are treated as a separate physical process and coupled with other processes using sequential splitting, given the same emission flux and same timestep size, a thinner bottom layer means a higher layer-mean mixing ratio will result after the surface emissions are applied.The left and right parts of the schematic correspond to the L30 and L72 vertical grids depicted in Fig. 2, respectively.A discussion of this schematic can be found in Sect.3.2.

A simple revision to the original coupling scheme
The discussions above motivates an alternate numerical scheme that couples the surface emissions of dust (as well as other aerosol species) more tightly with the turbulent mixing.From a mathematical perspective, the ideal approach would be to provide the surface emission fluxes as a boundary condition to the equations of turbulent mixing so that the two processes (emissions and mixing) can be solved together numerically.This approach would require a significant amount of code modifications and is deferred to future work.Here we take a simpler and admittedly less optimal method that requires the least amount of code modifications, namely to move the update of aerosol and precursor gas mixing ratios using cam_in%cflx from the original location to right before the cloud macrophysics-microphysics sub-cycles.In other words, box 3 in Fig. 1 is moved to the location indicated by the pink box with dashed outline, after box 6 (deep convection) and before box 7 (in which turbulent mixing is calculated).The emission fluxes, all other parameterizations, and the resolved dynamics are calculated at their original locations in the time loop.For aerosols, this simple modification still uses sequential splitting between emissions and turbulence-mixing-activation-resuspension, but the dry removal processes are calculated before the surface emissions are applied, and the wet removal processes are calculated after turbulent mixing; hence neither the dry removal nor the wet removal uses mixing ratios with spikes in the lowest model layer.

Impact of the revised process coupling on aerosol climatology in EAMv1
We now compare EAMv1 simulations conducted with the original and revised coupling schemes.The goal is twofold: (1) to verify the reasoning in the previous section about the features of the two schemes, and (2) to evaluate the impacts of the revision on the simulated aerosol climatology at regional and global scales.

Dust aerosols
From the discussions in Sect.3, we expect the direct effect of the revised coupling scheme to be a weakening of dry removal in grid cells and timesteps where dust emissions occur.Based on the understanding of process interactions in EAM, we can reason how the other aerosol processes can be affected.In the revised scheme, turbulent mixing is the first aerosol process calculated after surface emissions are applied.Since the newly emitted particles have not gone through dry removal, more (compared to the case in the original scheme) particles are available for upward transport by turbulence.After turbulent mixing, more aerosol particles can be expected in upper model layers in the source regions.These particles can go through wet removal, or get advected out of the atmosphere column by resolved winds.More transport from source to non-source regions can increase aerosol load in the non-source regions and consequently lead to more removal there.These expectations are confirmed by the EAM results shown below.
The weakening of dry removal in the lowest model layer in dust source regions can be seen in Fig. 9 which shows annual mean results in North Africa as an example.The emission-induced dust mixing ratio tendencies are shown in panel (a) to help identify locations with emissions.The changes in dry removal tendencies caused by the revision of coupling show spatial   patterns that closely match the emissions (Fig. 9c vs. 9a).The relative differences shown in Fig. 9d suggest that decreases exceeding 50% can be found in the majority of the North African dust source regions.
The process rate changes in dust source regions in model layers above the lowest can be seen from the tendencies profiles shown in Fig. 10.The figure is essentially the same as Fig. 3c but with the results from the revised coupling added as open circles.Like in Fig. 3c, the tendencies shown here are annual mean values averaged over dust source regions located in the magenta box in Fig. 3a.The comparison in Fig. 10 shows that when the revised coupling is used, large and systematic increases in magnitude are seen in turbulent mixing, resolved transport, and dry removal in the lower troposphere above the lowest model layer.Between 900 hPa and 1000 hPa, the relative magnitude increases are on the order of 100%.
The global impacts of the revised coupling can be seen in Fig. 11 which shows the annual mean dust load in the atmosphere.
The first row corresponds to the interstitial dust burden integrated from the Earth's surface to model top.The second and third rows are pressure-longitude cross-sections of dust mixing ratios meridionally averaged over the two latitude bands of 0-60 • N and 40 • S-10 • S, respectively.The burden or mixing ratio changes shown in the middle and right columns of the figure indicate In Sect.2.2.1, it has been explained that 5 out of the 7 aerosol species in MAM4, namely dust, sea salt, MOA, BC, and POA, have large portions (50% or more) of their sources coming from surface emissions.The other two species, sulfate and SOA, are predominantly or 100% produced in the atmosphere, as can be seen in the attribution of mass sources shown in Table 1.Therefore, we expect the revised coupling to have significant impacts on the first 5 aerosol species mentioned above but negligible direct impacts on sulfate and SOA.
In Table 5, we show the vertically integrated and annually averaged dry removal rates for regions with or without surface emissions as well as for the entire globe, for the 5 aerosol species with substantial surface sources.The relative differences listed in the table reveal consistent results among the different species, namely dry removal becomes weaker in regions with surface emissions but stronger in regions without surface emissions.The global averages, which are dominated by the results in the source regions, show a general trend of weaker dry removal when the revised coupling is used.

22
Figure 12.Comparison of the zonal and annual mean interstitial aerosol mass mixing ratios (unit: kg kg −1 ) simulated with EAMv1 using the original or revised process coupling and with the L72 vertical grid.The different rows correspond to the 5 aerosol species listed in Table 2.
Left column shows the L72 simulation using the revised coupling scheme.Middle column shows the differences caused by the revision in process coupling.Right column shows the relative differences with respect to the original model.For the reference differences plots, the locations with mean mixing ratios lower than the first contour level in the left panel are masked out.
It is worth noting that among those 5 aerosol species that have substantial surface emissions, dust and sea salt emissions introduce particles primarily to MAM4's coarse mode (where the mode's geometric mean diameter range is 1 µm ⩽ D p ⩽ 4 µm), while MOA, BC, and POA emissions add particles mostly to MAM4's primary carbon mode (10 nm ⩽ D p ⩽ 100 nm) and accumulation mode (53.5 nm ⩽ D p ⩽ 440 nm).This can be seen from the emission attribution shown in Table 2.Although there are various mechanisms for different aerosol species to mix and for aerosol particles to grow in the atmosphere, MOA, BC, and POA are predominantly found in aerosol particles of submicron diameters.Since smaller particles are less susceptible to dry removal, we expect to see smaller impacts of the revised coupling on the atmospheric load of MOA, BC, and POA than on dust and sea salt.This reasoning is confirmed by the global annual mean mass burdens shown in Table 6 and the zonally and annually averaged mass mixing ratios shown in Fig. 12.The global mean burdens increase by 39% and 52% for dust and sea salt, respectively, but less than 20% for MOA, BC, and POA.The relative changes in zonal mean mixing ratio shown in the rightmost column in Fig. 12 further confirm the large impacts on dust and sea salt and the smaller impacts on the predominantly submicron species (MOA, BC, POA).As for sulfate, SOA, and precursor gases in EAMv1, the changes in zonal and annual mean mixing ratios are typically less than 10% and hence not shown.

Sensitivity to bottom layer thickness
In Sect.3.2, we attributed the strong vertical resolution sensitivity of dust dry removal reported in Feng et al. (2022) to the original sequential splitting scheme that calculates dry removal after surface emissions and before turbulent mixing.The revised coupling scheme allows the emitted particles to be vertically mixed by turbulence before other processes are calculated.Since turbulent mixing is very efficient in reducing the vertical gradients and resulting in similar mixing ratios in model layers near the Earth's surface (Fig. 6), we expect the simulated dry removal rates to be less sensitive to the thickness of the lowest model layer when the revised coupling is used, and we expect the same qualitative results to be seen in all of the 5 aerosol species that have significant surface emissions.These expectations are confirmed by the results shown in Table 7.
Furthermore, the global and annual mean aerosol lifetimes shown in Table 8 indicate that when the revised coupling is used, the relative differences between results from the L30 and L72 simulations are reduced to less than 10% in contrast to 32% for dust and 47% for sea salt in the original EAMv1.

Summary, conclusions, and outlook
The earlier work by Feng et al. (2022) pointed out that various aspects of the dust aerosol life cycle simulated by EAMv1, in particular the dry removal fluxes and lifetime, were sensitive to the use of 72 versus 30 layers for the vertical grid.In this paper, we investigated the resolution sensitivity by carrying out detailed budget analyses for the interstitial dust mass mixing ratio and by tracking the evolution of the mixing ratio within the model's time integration cycle.Dust emissions in the real world typically occur in turbulent ambient conditions, hence the emitted particles can be efficiently distributed to a significant depth of the atmosphere column.In EAMv1, the numerical coupling method treats surface emissions as a separate process which adds aerosol particles to the lowest model layer, while dry removal and turbulent mixing are calculated after the emissions are applied, using the sequential splitting method.This ordering of processes results in high, unrealistic temporary values of aerosol mixing ratios to be provided as input to the parameterization of dry removal, causing large numerical errors in the simulated dry removal rates.The situation gets worse when the model's bottom layer is thinner, as the same emission fluxes will result in higher temporary mixing ratios in the bottom layer when emissions are applied in isolation.
Based on this reasoning, we proposed a simple revision to the numerical process coupling in EAMv1, namely to move the application of emissions to the location right before turbulent mixing in the model's time integration loop.The revision allows the newly emitted aerosol particles to be vertically distributed by turbulence before experiencing other processes considered in the model, hence giving a numerical coupling scheme with closer resemblance to the process interactions in the real world.
The revised coupling scheme was implemented in EAMv1, and wind-nudged simulations were conducted for the year 2010 using 72 or 30 grid layers with 1 • horizontal grid spacing.As expected, the revision substantially weakened dry removal and strengthened vertical mixing of dust in its source regions, strengthened the large-scale transport from source to non-source regions, strengthened dry removal outside the source regions, and strengthened activation and wet removal globally.When using 72 grid layers without retuning of uncertain parameters, the revised process coupling was found to cause a 39% increase in the global annual mean dust burden and an increase of dust lifetime from 1.9 days to 2.6 days.
The revised process coupling was implemented for all aerosol species in EAMv1 and was found to cause the same qualitative changes also for other aerosol species that have substantial surface emissions, i.e., sea salt, MOA, BC, and POA.Quantitatively, the changes in global mean burden and lifetime were considerably smaller for the predominantly submicron species (MOA, BC, POA), with burden and lifetime increases less than 20%, whereas the increases in dust and sea salt burden and lifetime exceeded 30% and 50%, respectively (Table 6).The impact on sulfate, SOA, and precursor gases were found to be very small.
The revised process coupling proposed in this study allowed us to use a minimal amount of code changes to obtain significant improvements in the EAMv1 results.However, both the surface emissions and dry removal of aerosols are still sequentially split from turbulent mixing using relatively long coupling timesteps of 30 minutes.This is likely a significant contributor to the remaining 17% difference in the annual mean source-region mean dust dry removal rate between the L72 and L30 simulations when the revised coupling is used (Table 7).In the next steps, it will be worth exploring additional revisions to further reduce the numerical errors caused by process splitting.For example, since emissions in the real world are typically accompanied by turbulent mixing, we plan to implement and evaluate a scheme that combines the time integration of emissions and turbulent mixing.Since turbulent dry deposition is also closely related to turbulent mixing, it is worth considering combining these two processes as well.In contrast, it is less obvious which coupling methods will be cost-effective in reducing the splitting errors related to gravitational settling.Different options will be explored in future work.We do not expect that vertical resolution sensitivities in the simulated aerosol life cycles will be completely eliminated merely by further revising the numerical coupling of aerosol processes discussed in this paper, because discretization errors in the individual parameterizations, as well as discretization errors and numerical coupling errors associated with other processes in EAM (e.g., clouds), can also cause vertical resolution sensitivities, and some of the sensitivities can be advantageous as they demonstrate the benefits of using higher resolutions.Nevertheless, further work on numerical coupling will likely make useful contributions to the goal of reducing numerical errors in EAM simulations.3).The mixing ratio tendencies in the upper and lower rows have been multiplied by 500 and 50,000, respectively, in order to be plotted with the same x-axis range as in Fig. 5.

Figure A3
. As in Fig. 6 but showing statistical distributions of instantaneous dust mass mixing ratios in the source-vicinity regions (upper row) and remote regions (lower row).The mixing ratios in the upper and lower rows have been multiplied by 10 and 1,000, respectively, in order to be plotted with the same x-axis range as in Fig. 6.

Figure 1 .
Figure 1.A schematic showing the sequence of calculations in a 30-minute timestep in 1 • EAMv1 simulations.Rectangles are calculations that do not directly affect aerosol mixing ratios.Boxes with round corners are calculations that change aerosol mixing ratios.Colored boxes correspond to the physical processes for which numerical results are shown in this paper.The numbers in black squares correspond to the numbering in Sect.2.3.The dashed pink box between deep convection (box 6) and the cloud sub-cycles (box 7) is where box 3 in the originalEAMv1 is moved to when the revised coupling scheme is used.

Figure 2 .
Figure 2. Near-surface layers of the L30 and L72 vertical grids used by EAMv1.Each gray box with a solid outline and a number at the center represents one model layer.The numbers are layer indices.The dashed lines indicate locations of layer midpoints.The thick horizontal lines at 1000 hPa nominal pressure represent the Earth's surface.Here, nominal pressure refers to the air pressure at a layer interface or midpoint when the air pressure at the Earth's surface is 1000 hPa.

Figure 3 .
Figure 3. One-year averages of various quantities in simulation v1_ori_L72.Left column: surface emission flux (panel a) and burden (panel b) of dust aerosol mass.Other columns: tendencies (i.e., rates of change) of interstitial dust aerosol mass mixing ratio caused by various physical processes, averaged over dust source regions (panels c, d), source-vicinity regions (panels e, f), and remote regions (panels g, h).Source regions and source-vicinity regions are defined as grid columns within the magenta box in panel (a) in which the annual mean emission fluxes are non-zero and zero, respectively.Remote regions refer to grid columns inside the brown box in panel (a).Line plots in the upper row show results of the lower troposphere, with the bottom layer of the model's vertical grid (i.e., the grid layer closest to the Earth's surface) excluded.Line plots in the lower row are results in the lowest 10 model layers including the bottom layer.Note that the tendencies shown in panels (c) and (e)-(h) have been multiplied by factors of 100 to 100,000, as noted at the top of each panel, in order for the same x-axis range to be used for all profile plots.

Figs.
Figs.A1 and A2, respectively.These figures suggest the dominant sources and sinks of dust mass on a timestep-by-timestep and grid-column-by-grid-column basis are the same as what we saw in the regional and annual averages, namely the following processes in the lowest model layer in the dust source regions: (1) surface emissions, (2) dry removal, and (3) turbulent mixing and aerosol activation-resuspension.

Figure 4 .Figure 5 .
Figure 4. Histogram of dust emission fluxes derived from 6-hourly instantaneous model output in the 90-day period of 2010-01-29 to 2010-04-28 in grid columns with non-zero dust emissions in the magenta box marked in Fig. 3a.

Figure 6 .
Figure 6.Similar to Fig. 5 but showing the statistical distributions of instantaneous dust mass mixing ratios after the tendencies shown in Fig. 5 have been applied.Results here correspond to the middle portion of the emission histogram in Fig. 4. Results for the source-vicinity and remote regions are shown in Fig. A3.

Figure 8 .
Figure 8. Similar to Figs. 5 and 6 but showing results from simulation v1_ori_L30 which used the L30 grid instead of L72.The upper row shows statistical distributions of the dust mass mixing ratio tendencies caused by different aerosol processes.The lower row shows the statistical distributions of dust mass mixing ratios after the corresponding tendencies in the upper row have been applied.Only 5 model layers are shown in this figure, as Fig. 2 suggests that the 5 lowest layers in the L30 grid cover roughly the same altitude range as the 10 lowest layers in the L72 grid.

Figure 9 .
Figure 9. (a) and (b): Annual mean interstitial dust mass mixing ratio tendencies (unit: kg kg −1 s −1 ) in the lowest model layer in simulation v1_ori_L72 caused by surface emissions and dry removal.(c) and (d):The differences (unit: kg kg −1 s −1 ) and relative differences (unit: fraction) in dust dry removal between simulations using the revised and original process coupling schemes.

Figure 10 .
Figure 10.Same as Fig. 3c but comparing results obtained with the original coupling scheme (filled circles) and the revised coupling scheme (open circles).Shown are annual mean results averaged over dust source regions located in the magenta box in Fig. 3a.

Figure 11 .
Figure11.Comparison of annual mean interstitial dust load simulated with the original or revised process coupling scheme in EAMv1 using the L72 grid.Top row: global distribution of vertically integrated dust burden (unit: kg m −2 ).Second and third rows: longitudepressure cross-sections of interstitial dust mass mixing ratio (unit: kg kg −1 ) averaged over the latitude bands of 0 • -60 • N and 40 • S-10 • S, respectively.Left column shows the L72 simulation using the revised coupling scheme (v1_rev_L72).Middle column shows the changes caused by the revised coupling (i.e., v1_rev_L72 minus v1_ori_L72).Right column shows the relative differences with respect to the original model.

Figure A2 .
Figure A2.As in Figs. 5 and A1 but showing results in the source-vicinity regions (upper row) and remote regions (lower row).The sourcevicinity regions are defined as grid columns and timesteps with zero dust emission in the 90-day period of 2010-01-29 to 2010-04-28 in the magenta box marked in Fig.3a.The remote regions are defined as the brown box in Fig.3a.The statistical distributions were derived from 6-hourly instantaneous output from the default EAMv1 (i.e., simulation v1_ori_L72 in Table3).The mixing ratio tendencies in the upper

Table 3 .
Simulations conducted in this study.Further details can be found in Sect.2.4.

Table 5 .
Similar to Table4but showing the impact of the revised process coupling on dry removal of aerosol species that have substantial surface emissions.Here, source region and non-source region refer to geographical locations with or without surface emissions of the corresponding aerosol species.Dry removal rates are shown in Tg yr −1 .Relative differences are shown in percentages.

Table 6 .
Similar to Table5but showing the impact of the revised process coupling on the global mean burden (unit: Tg) and lifetime (unit: day) of different aerosol species.Relative differences are shown in percentages.

Table 7 .
Sensitivity of the annual (year-2010) mean source-region mean aerosol dry removal rates (unit: Tg yr −1 ) to vertical resolution when using the original or revised process coupling in EAMv1.Results are shown for the 5 aerosol species in MAM4 that have significant surface emissions.

Table 8 .
Sensitivity of the global and annual (year-2010) mean aerosol lifetime (unit: day) to vertical resolution when using the original or revised process coupling in EAMv1.Results are shown for the 5 aerosol species in MAM4 that have significant surface emissions.