2.31. Dust Emission

Atmospheric dust is mobilized from the land by wind in the CLM. The most important factors determining soil erodibility and dust emission include the wind friction velocity, the vegetation cover, and the soil moisture. The latest CTSM allows users to choose between two dust emission schemes: One is Leung_2023 (Leung et al. 2023; Leung et al. 2024) which is the current default for the CLM6 physics or later, and the other is Zender_2003 (Mahowald et al. 2006) based on the DEAD (Dust Entrainment and Deposition model of Zender et al. (2003), which is the default for the CLM5 or older physics.

We here describe the Leung_2023 scheme based on Leung et al. 2023 and Leung et al. (2024) and document some differences in tuning in the latest CTSM. CTSM users can look for the previous documentation for CLM5 physics for a description of the Zender_2003 scheme.

2.31.1. Dust Emission Thresholds

2.31.1.1. Fluid Threshold

Dust emission modeling is a threshold parameterization of an aeolian (wind-driven) process for both Leung_2023 and Zender_2003. For Leung_2023, in any given timestep the soil surface wind velocity \(u_{*s}\) has to be greater than the fluid threshold friction velocity \(u_{*ft}\) to generate saltation and dust emission:

(2.31.1)\[u_{*ft} = u_{*ft0}(D_{p},\rho_{atm}) f_{m}(w)\]

where \(u_{*ft0}(D_{p},\rho_{a})\) is the dry fluid threshold without the soil moisture effect \(f_{m}\), as a function of median soil diameter \(D_{p}\) and air density \(\rho_{atm}\). In CTSM for Leung_2023, \(D_{p}\) is a globally uniform number of 130 \(\mu\) m. \(u_{*ft0}(D_{p},\rho_{atm})\) is given by Shao and Lu (2000):

(2.31.2)\[u_{*ft0}(D_{p},\rho_{atm}) = \sqrt{\frac{A(\rho_{p} g D_{p} + \gamma / D_{p}) }{\rho_{atm}} }\]

where g is the acceleration of gravity (Table 2.2.7), \(\rho_{p} = 2650\) kg m^-3 is typical soil particle density, and \(A = 0.0123\) and \(\gamma = 1.65 \times 10^{-4}\) kg s^-2 are empirical constants.

2.31.1.2. Impact of Soil Moisture

The soil moisture effect \(f_{m}(w)\) is a function of gravimetric soil moisture \(w\) (kg water / kg soil) at the topmost soil layer.

\(w\) is converted from the CLM volumetric soil moisture \(\theta\) and porosity (saturation moisture \(\theta_{sat}\)) at the topmost soil layer:

(2.31.3)\[w=\theta\frac{ \rho _{liq} }{\rho_{bulk} }\]

Note that \(w\) or \(\theta\) in CTSM is conventionally (Mahowald et al. 2006) treated as a sum of both liquid and ice/frozen soil moisture at the topmost soil layer, i.e., \(w_{1} = w_{liq,1} + w_{ice,1}\). We skip \(w_{1}\) and use \(w\) for simplicity. \(\theta\) is the volumetric soil moisture (water+ice) in the topmost soil layer (m^-3water m^-3 soil) (section 2.7.3), \(\rho _{liq}\) is the density of liquid water (kg m^-3) (Table 2.2.7), and \(\rho _{bulk}\) is the bulk density of soil in the top soil layer (kg m^-3) defined as in section 2.6.3 rather than as in Zender et al. (2003). \(\rho_{bulk}\) is given by

(2.31.4)\[\rho_{bulk} = (1 - \theta_{sat} ) \rho_{p}\]

Then, the soil moisture effect \(f_{m}(w)\) on increasing the fluid threshold is given by Fecan et al. (1999):

(2.31.5)\[\begin{split}f_{m}(w) =\left\{\begin{array}{l} {1{\rm \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; for\; }w\le w_{t} } \\ {\sqrt{1+1.21\left[100\left(w-w_{t} \right)\right]^{0.68} } {\rm \; \; for\; }w>w_{t} } \end{array}\right.\end{split}\]

where \(w_{t}\) is the minimum gravimetric moisture threshold required to increase the interparticle force and \(u_{*ft}\). \(w_{t}\) increases with clay fraction \(f_{clay}\):

(2.31.6)\[w_{t} =0.01a\left(17f_{clay} +14f_{clay}^{2} \right){\rm \; \; \; \; \; \; 0}\le f_{clay} =\% clay\times 0.01\le 1\]

where \(a=f_{clay}^{-1}\) for tuning purposes. Note that this is different from the paper (Leung et al. 2023; Leung et al. 2024) in which \(a=1\) was chosen. The coefficient 0.01 is used for converting \(w_{t}\) from % to fraction (kg water / kg soil). \(f_{clay}\) is the mass fraction of clay particles in the topmost soil layer and %clay comes from the surface dataset (section 2.2.3.3).

2.31.1.3. Impact Threshold

Another essential dust emission threshold is the impact/dynamic threshold \(u_{*it}\), which is the lowest friction velocity or wind stress to matintain saltation:

(2.31.7)\[u_{*it} = B_{*it} u_{*ft0}\]

where \(B_{*it}\) is a constant on Earth following Kok et al. (2012). In Leung_2023, \(u_{*it}\) does not depend on and increase with soil moisture. The above equations imply that \(u_{*it} \, < \, u_{*ft0} \le \, u_{*ft}\). This means that the winds need a bigger momentum to initiate saltation and dust emission but can reduce below \(u_{*ft}\) and still maintain a weak dust emission flux. The emission flux goes to zero when \(u_{*s}\) drops below \(u_{*it}\).

2.31.2. Dust Emission Flux

The total vertical mass emission flux of dust, \(F_{d}\) (kg m^-2 s^-1), from the ground into a transport mode/bin \(j\) of aerosol is based on Kok et al. (2014a):

(2.31.8)\[F_{d} = \eta C_{tune} C_{d} f_{bare} f_{clay'} \frac{ \rho_{atm} (u^2_{*s} - u^2_{it} ) }{ u^2_{it} } \left( \frac{ u^2_{*s} }{u^2_{it} } \right) ^\kappa\]

where \(C_{tune} = 0.05\) is a constant, \(\eta\) is the intermittency factor (we will derive it in section 2.31.4), and \(F_{d}\) is the total emission flux summed across modes/bins following a revised form of Kok et al. (2014b). The dust emission flux goes to zero when \(u_{*s} \, < \, u_{*it}\). \(\rho_{atm}\) is surface air density from CAM (the atm model). \(f_{clay'}\) is a modified clay fraction term appeared earlier in Zender et al. (2003). In Zender_2003 it is used to indicate the sandblasting efficiency. Zender limited this term to be capped at 0.2:

(2.31.9)\[f_{clay'} = \textnormal{max}(f_{clay}, 0.2)\]

\(f_{clay'}\) is later adopted by Kok et al. (2014a) to indicate the amount of fine dust particles for sandblasting. The same \(f_{clay'}\) is used in Leung et al. (2024). But we later further limit \(f_{clay'}\) to be within 0.1 and 0.2 to further reduce the impact of how the abundance of fine dust particles can scale the total dust emission flux, as part of the CESM tuning:

(2.31.10)\[f_{clay'} = \textnormal{max}(0.1+0.5f_{clay}, 0.2)\]

Then, \(\kappa\) is the fragmentation exponent, and \(C_{d}\) is the dust emission coefficient (or the soil erodibility coefficient):

(2.31.11)\[C_{d} = C_{d0} \exp{ (-C_{e} \frac{ u_{*st} - u_{st0} }{ u_{st0} } ) }\]

(2.31.12)\[\kappa = C_{\kappa} \frac{ u_{*st} - u_{st0} }{ u_{st0} }\]

where \(C_{\kappa} = 2.7\), \(u_{st0} = 0.16\) m s ^-3, \(C_{d0} = 4.4 \times 10^{-5}\), and \(C_{e} = 2.0\). \(F_{d}\) thus roughly scales with \(u^{2+\kappa}_{*s}\), where \(\kappa \sim 1\) over major deserts and \(\sim 3\) or higher over semiarid and nonarid regions. Since Kok et al. (2014a) has not measured \(\kappa > 3\) in their measurements, we cap \(\kappa\) at a maximum value (currently set as 2.5). \(u_{*st}\) is the standardized wet fluid threshold at a typical atmospheric surface air density (Kok et al., 2014):

(2.31.13)\[u_{*st} = u_{*ft} \sqrt{ \rho_{atm} / \rho_{0atm}}\]

where \(\rho_{0atm} = 1.225\) kg m^-3. As can be seen, \(u_{*st}\) scales with \(u_{*ft}\) and thus soil moisture \(w\). Therefore, moisture \(w\) decreases soil erodibility \(C_{d}\) but increases dust emission sensitivity \(\kappa\) to the winds.

Kok et al. (2014a) is different from many other dust emission parameterizations in the way that the soil erodibility \(C_{d}\) is not a time-invariant input data but is a transient function, with erodibility increasing with reducing \(u^2_{*ft}\) (and thus implicitly soil moisture). Similarly, the fragmentation exponent \(\kappa\) is also transient and increases with enhancing soil moisture.

The grid cell fraction of exposed bare soil suitable for dust mobilization \(f_{bare}\) is given by

(2.31.14)\[f_{bare} =\left(1-f_{lake} \right)\left(1-f_{sno} \right)\left(1-f_{v} \right)\frac{w_{liq,1} }{w_{liq,1} +w_{ice,1} }\]

where \(f_{lake}\) and \(f_{sno}\) are the CLM grid cell fractions of lake (section 2.2.3.3) and snow cover (section 2.8.1), all ranging from zero to one. Not mentioned by Zender et al. (2003), \(w_{liq,\, 1}\) and \(w_{ice,\, 1}\) are the CLM top soil layer liquid water and ice contents (mm) entered as a ratio expressing the decreasing ability of dust to mobilize from increasingly frozen soil. The grid cell fraction of vegetation cover,\(f_{v}\), is defined as

(2.31.15)\[0\le f_{v} =\frac{\mathrm{VAI}}{\mathrm{VAI_{thr}} } \le 1{\rm \; \; \; \; where\; } \mathrm{VAI_{thr}} =0.6{\rm \; m}^{2} {\rm m}^{-2}\]

where \(\mathrm{VAI}=\mathrm{LAI}+\mathrm{SAI}\) is the vegetation area index as a sum of the CLM leaf and stem area index values (m ² leaf m^-2 grid) averaged at the land unit level so as to include all the pfts and the bare ground present in a vegetated land unit. Currently, Leung_2023 in CTSM sets the areas with \(\mathrm{VAI_{thr}}\) smaller than or equal to 0.6 m ² m^-2 to be dust-emitting grids, different from the Leung et al. (2024) paper that set \(\mathrm{VAI_{thr}}\) to be 1 m ² m^-2. \(\mathrm{LAI}\) and \(\mathrm{SAI}\) may be prescribed from the CLM input data (section 2.2.1.4) or simulated by the CLM biogeochemistry model (section 2.21).

2.31.3. Drag Partition Effect On Reducing Wind Stress

On top of Kok et al. (2014a), Leung_2023 introduced the soil surface friction velocity \(u_{*s}\) as the friction velocity \(u_{*}\) from CLM corrected by the surface roughness due to the presented rocks and vegetation on the soil surface, encapsulated by the so-called drag partition factor \(F_{eff}\).

(2.31.16)\[u_{*s} = u_{*} F_{eff}\]

The Leung et al. (2023) paper uses an area-weighted averaging method to determine the mean drag partitioning for a grid cell: \(F_{eff}^3 = A_{rock} f_{rock}^3 + A_{veg} f_{veg}^3\). \(A_{rock}\) and \(A_{veg}\) are fractional area cover (in fraction) from the CLM-prescribed land use from the Land Use Harmonization 2 (LUH2; section 2.28). \(F_{eff}\) is thus a weighted mean of the rock drag partitioning and the vegetation drag partitioning in Leung et al. (2023). However, since CTSM has the privilege of supporting sub-grid patch-level simulations of dust emissions, we simply separate the calculations of dust emissions into the areas of bare soils and areas of the short vegetation. For a bare soil patch/PFT we use:

(2.31.17)\[ F_{eff} = f_{rock}\]

And for a patch/PFT with short vegetation (shrub, grass, crop) we use:

(2.31.18)\[ F_{eff} = f_{veg}\]

The rock drag partition factor scales with the surface roughness density of rocks, captured by the aeolian roughness length \(z_{0a}\) from the satellite-derived dataset from Prigent et al. (2005). The expression was first developed by Marticorena and Bergametti (1995), which is more valid for the low-roughness surfaces:

(2.31.19)\[f_{rock} = 1 - \frac{\ln(\frac{z_{0a}}{z_{0s}})}{\ln[b_{1}(\frac{X}{z_{0s}})^{b_{2}}]}\]

where \(X = 10\) m is the distance downstream the point of discontinuity in surface obstacle, \(b_{1} = 0.7\) and \(b_{2} = 0.8\) are coefficients (Darmenova et al., 2009), \(z_{0s}\) is the soil roughness length. \(z_{0a}\) is from Prigent et al. (2005) and should not be confused with the aerodynamic roughness length \(Z_{0}\) from the model. This equation only applies for gridcells with VAI smaller than the VAI threshold for dust emission.

The vegetation drag partition factor scales with vegetation density as captured by VAI following Okin (2008) and Pierre et al. (2014):

(2.31.20)\[f_{veg} = \frac{K + f_{0} c }{K + c}\]

(2.31.21)\[K = 2 (\frac{1}{f_{v}} - 1) = 2 (\frac{1}{\mathrm{VAI}/\mathrm{VAI_{thr}}} - 1)\]

where \(f_{0} = 0.33\) is the wind dissipation right behind the obstacle, \(c = 4.8\) m is the average e-folding distance for wind momentum to recover after dissipated by the plant, and \(K\) is defined as the average normalized gap length between two obstacles given the vegetation density quantified by VAI. This equation assumes that the gap length K goes to zero as VAI increases and approaches \(\mathrm{VAI_{thr}}\).

2.31.4. Emission Intermittency Due To Turbulent Wind Fluctuations

Then, the fraction of time within a timestep \(\Delta t\) (typically of 30 minutes for a ~1° horizontal grid) when dust emission is active is accounted for by Leung_2023, known as the intermittency factor \(\eta\). The fraction of time \(\eta \in [0,1]\) is due to the boundary-layer turbulence that causes sub-timestep turbulent wind fluctuations that cross the dust emission thresholds multiple times within a timestep, in turn initiating and shutting off emission at time within \(\Delta t\). The parameterization of \(\eta\) comes from Comola et al. (2019), in which a statistical substepping method was proposed to account for the temporary shutoff of dust emission fluxes.

The fraction of time \(\eta\) is parameterized using the surface winds and thresholds at the saltation height. Therefore, the friction velocities are translated using the log law of the wall to the saltation height, which was defined as \(z_{sal}\) = 0.1 m by Comola et al. (2019):

(2.31.22)\[u_{s} = \frac{u_{*s}}{k} \ln(z_{sal}/z_{0a})\]

(2.31.23)\[u_{ft} = \frac{u_{*ft}}{k} \ln(z_{sal}/z_{0a})\]

(2.31.24)\[u_{it} = \frac{u_{*it}}{k} \ln(z_{sal}/z_{0a})\]

where k is the von Karman constant (Table 2.2.7), and \(z_{0a}\), the aeolian roughness length, is set to be 10 ^-4 m here for simplicity. With saltation-height variables defined, the instantaneous wind \(\tilde{u}_s\) is assumed by Comola to follow a Gaussian distribution with a mean equal to the mean wind speed and the spread \(\sigma_{u_{s}}\) parameterized by the Similarity Theory (Panofsky et al., 1977):

(2.31.25)\[\tilde{u}_s \sim N(u_s, \sigma_{u_s})\]

And the fluctuation strength is parameterized by the similarity theory:

(2.31.26)\[\sigma_{u_s} = u_{*s} \left( 12 - 0.5 \frac{z_i}{L} \right)^{1/3} \quad \text{for } 12 - 0.5 \frac{z_i}{L} \ge 0\]

where \(z_i = 1000\) m is the planetary boundary-layer height set as a constant for now, and \(L\) is the Obukhov length scale. This means the instantaneous wind's fluctuation comes from both a shear contribution and a buoyancy contribution.

Then, the total fraction of time \(\eta\) when saltation is active within a model timestep is then formulated as

(2.31.27)\[\eta = 1 - P_{ft} + \alpha \left( P_{ft} - P_{it} \right)\]

where \(P_{it}\) is the cumulative probability that the instantaneous wind \(\tilde{u}_s\) does not exceed the impact threshold \(u_{it}\), and \(P_{ft}\) is the cumulative probability that \(\tilde{u}_s\) does not exceed the fluid threshold \(u_{ft}\). The fluid threshold crossing fraction \(\alpha\) is defined as the rate of \(\tilde{u}_s\) sweeping across \(u_{ft}\) divided by the rate of sweeping across \(u_{it}\) and \(u_{ft}\) summed up. The detailed physical interpretation of this formula is in Leung et al. (2023). Here we just document the equations for each term:

(2.31.28)\[\alpha \approx \left\{ \exp\left[ \frac{u_{ft}^2 - u_{it}^2 - 2u_s(u_{ft}-u_{it})}{2 \sigma^2_{u_s} } \right] + 1 \right\}^{-1}\]

Then, the fraction of time in \(\Delta t\) when \(\tilde{u}_s\) is above \(u_{ft}\) is given by \(1 - P_{ft}\), where

(2.31.29)\[P_{ft} = \frac{1}{2} \left[ 1 + \operatorname{erf} \left( \frac{u_{ft} - u_s}{\sqrt{2} \sigma_{u_s}} \right) \right]\]

And, the fraction of time in \(\Delta t\) when \(\tilde{u}_s\) is dropping from above \(u_{ft}\) to within the hysteresis zone between \(u_{it}\) and \(u_{it}\) is given by \(\alpha \left( P_{ft} - P_{it} \right)\), where

(2.31.30)\[P_{it} = \frac{1}{2} \left[ 1 + \operatorname{erf} \left( \frac{u_{it} - u_s}{\sqrt{2} \sigma_{u_s}} \right) \right]\]

And so the fraction of time \(\eta\) within \(\Delta t\) with active emission is determined for (2.31.8).

2.31.5. Emitted Dust Size Distribution And Dust Transport In Atmosphere

The total vertical mass emission flux of dust, \(F_{d}\) (kg m^-2 s^-1) is then computed with all the above equations and terms. The emission flux is then passed through the coupler to the atmospheric model (CAM) to simulate dust aerosol transport and deposition. The default aerosol model supported in the CAM6 and CAM7 physics is the Modal Aerosol Model (MAM). CAM7 uses the 5-mode MAM (MAM5), in which three modes (Aitken, accumulation, and coarse modes) by default contain dust.

2.31.5.2. Conventional Bin Partition

In the early CESM versions, CAM employed the Bulk Aerosol Model (BAM) as the default aerosol model. Thus, the emission fluxes in CTSM is by default partitioned into 4 bins before passing to the coupler:

(2.31.31)\[F_{j} = F_d \sum _{i=1}^{I}M_{i,j}\]

where \(F_{j}\) is the mass emission flux from the \(j\) th aerosol bin. The current way of paritioning the emission fluxes before passing to the coupler is still being used, but the partition into the different bins is only required by BAM, not the current default MAM in CAM6 and CAM7. Therefore, for MAM, the four \(F_{j}\) are summed up to become one total flux \(F_d\) again inside MAM. It is then redistributed inside MAM to different individual MAM modes following the BFT. So, the following details are inherited from previous versions of the code and still works fine, but we may clean up the code in the future.

Equation (2.31.31) comes from Zender et al. (2003) and used by Mahowald et al. (2006). It sums \(M_{i,\, j}\) over \(I=3\) source modes \(i\) where \(M_{i,\, j}\) is the mass fraction of each source mode \(i\) carried in each of \(J=4\) transport bins \(j\)

(2.31.32)\[M_{i,j} =\frac{m_{i} }{2} \left[{\rm erf}\left(\frac{\ln {\textstyle\frac{D_{j,\max } }{\tilde{D}_{v,i} }} }{\sqrt{2} \ln \sigma _{g,i} } \right)-{\rm erf}\left(\frac{\ln {\textstyle\frac{D_{j,\min } }{\tilde{D}_{v,i} }} }{\sqrt{2} \ln \sigma _{g,i} } \right)\right]\]

where \(m_{i}\), \(\tilde{D}_{v,\, i}\), and \(\sigma _{g,\, i}\) are the mass fraction, mass median diameter, and geometric standard deviation assigned to each particle source mode \(i\) (Table 2.31.1), while \(D_{j,\, \min }\) and \(D_{j,\, \max }\) are the minimum and maximum diameters (m) in each transport bin \(j\) (Table 2.31.2).

Note that in CAM, dust emission flux will be scaled by another global dust tuning factor for matching the observed atmospheric dust constraints.

Table 2.31.1 Mass fraction \(m_{i}\) , mass median diameter \(\tilde{D}_{v,\, i}\) , and geometric standard deviation \(\sigma _{g,\, i}\) , per dust source mode \(i\)
\(i\)	\(m_{i}\) (fraction)	\(\tilde{D}_{v,\, i}\) (m)	\(\sigma _{g,\, i}\)
1	0.036	0.832 x 10\({}^{-6}\)	2.1
2	0.957	4.820 x 10\({}^{-6}\)	1.9
3	0.007	19.38 x 10\({}^{-6}\)	1.6

Table 2.31.2 Minimum and maximum particle diameters in each dust transport bin \(j\)
\(j\)	\(D_{j,\, \min }\) (m)	\(D_{j,\, \max }\) (m)
1	0.1 x 10\({}^{-6}\)	1.0 x 10\({}^{-6}\)
2	1.0 x 10\({}^{-6}\)	2.5 x 10\({}^{-6}\)
3	2.5 x 10\({}^{-6}\)	5.0 x 10\({}^{-6}\)
4	5.0 x 10\({}^{-6}\)	10.0 x 10\({}^{-6}\)