Dust Emission and Transport#
Warning
This site is Work in Progress.
Theoretical Concept#
Introduction#
With the Dust Emission and Transport (DET) module, PALM is able to simulate saltationinduced dust release and transport for different dustsized particles. The implementation and physical background is based on the publications of Klamt et al. (2024) and the AFWA scheme of LeGrand et al. (2019).
Attention: In contrast to Klamt et al. (2024), where the molecular viscosity of air has been calculated at each grid point and for each time step, a constant viscosity for an air temperature defined via namelist parameter pt_surface is used. This way, the performance of the DET is improved by one order of magnitude, while the effect on the results is only marginal. Furthermore, the calculated settling velocity is also based on pt_surface, e.g. the height dependence of temperature is ignored.
All particles with a diameter less than \(63\,\mu\)m are interpreted as dust, i.e., clay (\(\leq 4\,\mu\)m) and silt (\(4\,\mu\)m \(< D \leq 63\,\mu\)m) particles belong to this group (Klose, 2014). In general, the implementation follows an Eulerian approach, i.e., the simulation of dust transport is realized using a Cartesian 3D grid, and all calculations are performed at the grid points of this grid. The following dustrelated processes are considered in the implementation: dust release from the (flat/smooth) ground due to saltation of particles, gravitational settling, dry deposition, largescale subsidence, passive advection with the resolved turbulent flow, and subgridscale turbulent transport. Prognostic variables in DET are the mass concentrations (in kg m^{3}) of different dust size bins that can be defined by the user.
Dust Emission#
The parameterization of dust emission is based on the Air Force Weather Agency (AFWA) dust emission scheme of the GOCART aerosol model LeGrand et al. (2019), which is part of WRFChem (the Weather Research and Forecasting (WRF) model coupled with Chemistry). In the AFWA scheme, dust emission is handled as a twopart process, wherein large particle saltation from coarser dust and sandsized particles is triggered by wind shear and leads to a fineparticle (dustsized) bulk emission flux by saltation bombardment and aggregate disintegration. The bulk dust emission flux is then further distributed among different size bins.
For estimating the vertical bulk dust emission flux, particlesize dependent horizontal saltation fluxes are calculated for each given saltation size bin. These saltation fluxes are then weighted depending on the soil configuration to form the socalled total verticallyintegrated streamwise (horizontal) saltation flux. The default configuration of the saltation size bins is given in Table 1. The soil category 1 ("sand") of the STATSGOFAO Database (Pérez et al., 2011) is used here. This category assumes 92% sand particles, 5% silt particles, and 3% clay particles. The mass weighting factors are the result of the product of these percentage values and the mass fractions given in Table 1.
Saltation size bin  1  2  3  4  5  6  7  8  9  10 

Effective diameter (\(\mu\)m)  1.42  8  20  32  44  70  130  200  620  1500 
Soil separate class  Clay  Silt  Silt  Silt  Silt  Sand  Sand  Sand  Sand  Sand 
Particle density (kg m^{3})  2500  2650  2650  2650  2650  2650  2650  2650  2650  2650 
Mass fraction  1  0.25  0.25  0.25  0.25  0.0205  0.0410  0.0359  0.3897  0.5128 
Mass weighting factor  0.03  0.0125  0.0125  0.0125  0.0125  0.0189  0.0377  0.0330  0.3585  0.4718 
Table 1: Configuration of Saltation Size Bins and Associated Attributes for the AFWA Scheme Assuming the Soil Category 1 ("sand") of the STATSGOFAO Database (Pérez et al., 2011; LeGrand et al. (2019)).
Once the vertical bulk dust emission flux is calculated, it is distributed among different dust size bins by using the socalled brittle fragmentation theory (Kok, 2011). This theory enables the calculation of dust distribution weighting factors that are multiplied with the vertical bulk dust emission flux to form the binspecific emission flux. This procedure requires information about the configuration of dust size bins, similar to the configuration of saltation size bins above. Table 2 shows the default configuration derived from LeGrand et al. (2019).
Dust size bin  1  2  3  4  5 

Lower bound diameter (\(\mu\)m)  0.2  2  3.6  6  12 
Upper bound diameter (\(\mu\)m)  2  3.6  6  12  20 
Effective diameter (\(\mu\)m)  1.46  2.8  4.8  9  16 
Particle density (kg m^{3})  2500  2650  2650  2650  2650 
Table 2: Configuration of Dust Size Bins and Associated Attributes for the AFWA Scheme (LeGrand et al. (2019)).
Gravitational Settling#
The vertical divergence of the gravitational settling flux causes a change in the dust mass concentration which is why this process is considered in DET. The gravitational settling flux at a given position \((x,y,z)\) and time \(t\) is given in the Stoke's regime (laminar regime) by
with the dust mass concentration \(c\) and the gravitational settling velocity \(v_\text{g}\), which is also known as the terminal velocity. This velocity is calculated based on Jacobson (2005) (Eq. (20.4)) via:
where \(\rho_{\text{a}}\) is the air density, \(D_{\text{d},j}\) is the dust particle effective diameter, \(\rho_{\text{d},j}\) is the dust particle density, \(\eta_\text{a}\) is the molecular (dynamic) viscosity of air, \(g\) is the gravitational acceleration, \(C_j\) is the Cunningham slipflow correction for small particles, and the index \(j\) marks a specific dust bin. According to Farrell and Sherman (2015), the assumption of the Stoke's regime (now commonly referred to as Stoke's law) is valid for particles with diameters less than \(100\,\mu\)m, i.e., the whole range of dustsized particles. Thus, it is not recommended to use DET for particles other than dust.
Deposition#
Instead of only considering the binspecific dust emission flux, a net surface flux is eventually calculated which is the sum of \(F_{\text{e},j}\), the emission flux, and the deposition flux \(F_{\text{d},j}\):
By definition, \(F_{\text{d},j}\) is negative at the surface (directed downwards). So far, wet deposition is not considered. The deposition flux is directly proportional to the local concentration of the dust size bin \(j\) (\(c_j\)) at a reference height (\(z_\text{r}\)) at which the binspecific dry deposition velocity \(v_{\text{d},j}\) is evaluated:
The dry deposition velocity tries to capture all physical processes governing dry deposition. The reference height at which \(c_j\) and \(v_{\text{d},j}\) are evaluated is \(zu(1)=dz/2\) for default surfaces at the ground with \(dz\) being the vertical grid spacing. It is equal to the height of the constant flux layer where Monin–Obukhov similarity theory is assumed.
To calculate \(v_{\text{d},j}\) (being downward positive), we follow the dry deposition scheme Z01 of WRFChem (Zeng et al., 2020), which is based on electrical resistance models of Slinn (1982) and Zhang et al. (2001). The dry deposition velocity is calculated for each dust size bin as
The overall transfer resistance is divided into an aerodynamic resistance above the canopy \(R_\text{a}\), accounting for the particle turbulent diffusion in the constant flux layer, and a surface resistance \(R_\text{s}\). The latter summarizes the quasilaminar sublayer (or boundary) resistance and the canopy resistance (or bulk surface resistance). Brownian diffusion, impaction, and the rebound effect are considered here. Interception is neglected.
Usage#
PALM can be started using the following example setup:
!
! INITIALIZATION PARAMETER NAMELIST
!
&initialization_parameters
!
! grid parameters
!
nx = 399, ! Number of gridboxes in xdirection (nx+1)
ny = 399, ! Number of gridboxes in ydirection (ny+1)
nz = 144, ! Number of gridboxes in zdirection (nz)
dx = 10.0, ! Size of single gridbox in xdirection
dy = 10.0, ! Size of single gridbox in ydirection
dz = 10.0, ! Size of single gridbox in zdirection
dz_stretch_level = 1200.0, ! Height (in m) where stretching starts
dz_stretch_factor = 1.08, ! dz(k+1) = dz(k) * dz_stretch_factor
!
! mode
!
large_scale_subsidence = .TRUE., ! Parameter to enable largescale subsidence
subs_vertical_gradient = 0.0023, 0.0, ! Gradient(s) of the profile for the
! largescale subsidence/ascent velocity
! (in (m/s) / 100 m)
subs_vertical_gradient_level = 0.0, 1000.0, ! Height level from which on the gradient
! for the subsidence velocity is effective
! (in m, uvgrid)
!
! initialization
!
initializing_actions = 'set_constant_profiles', ! initial conditions
ug_surface = 0.0, ! ucomp of geostrophic wind at surface
vg_surface = 0.0, ! vcomp of geostrophic wind at surface
pt_surface = 300.0, ! initial surface potential temp
pt_vertical_gradient = 0.0,
2.0, ! piecewise temp gradients
pt_vertical_gradient_level = 0.0,
1000.0, ! height level of temp gradients
reference_state = 'horizontal_average', ! The instantaneous horizontally
! averaged potential temperature
! profile will be used as reference
! state in the buoyancy term. Doing
! so the pressure pertubation
! will fluctuate around zero in a
! xyplane instead of an increasing
! p* value of 25Pa/50Pa/...
!
! boundary conditions
!
bc_lr = 'cyclic',
bc_ns = 'cyclic',
bc_uv_b = 'dirichlet', ! (noslip condition)
bc_uv_t = 'neumann', ! (freeslip condition)
surface_heatflux = 0.24, ! sensible heat flux at the bottom surface
bc_pt_b = 'neumann', ! required with surface_heatflux
roughness_length = 0.1, ! Roughness length (in m)
constant_flux_layer = .TRUE., ! Parameter to switch on a constant flux layer
! at the bottom boundary.
!
! numerics
!
rayleigh_damping_height = 1300.0, ! Well above inversion height
rayleigh_damping_factor = 0.1, ! 0.0 <= rayleigh_damping_factor <= 1.0
fft_method = 'fftw', ! http://www.fftw.org/
momentum_advec = 'wsscheme', ! 5th order upwind scheme
scalar_advec = 'wsscheme', ! 5th order upwind scheme
!
! physics
!
omega = 7.29212E5, ! Angular velocity of the rotating system (Earth, in rad/s).
latitude = 52.37, ! Geographical latitude (in degrees north) of Hannover
/ ! end of initialization parameter namelist
!
! RUNTIME PARAMETER NAMELIST
!
&runtime_parameters
!
! run steering
!
end_time = 3600.0, ! simulation time of the 3D model
npex = 16, npey = 16, ! Virtual processor topology.
create_disturbances = .TRUE., ! randomly perturbate horiz. velocity
dt_disturb = 150.0, ! interval for random perturbations
disturbance_energy_limit = 0.01, ! upper limit for perturbation energy
data_output_2d_on_each_pe = .FALSE., ! don't do 2D output on each MPI rank
!
! data output
!
netcdf_data_format = 2, ! netCDF 64bitoffset format
dt_run_control = 0.0, ! output interval for run control
dt_data_output = 900.0, ! output interval for general data
dt_data_output_av = 900.0, ! output interval for averaged data
dt_dopr = 900.0, ! output interval for profile data
nz_do3d = 60, ! 3d output up to 600m
data_output = 'u_xy', 'u_xz', 'u_xy_av', 'u_xz_av',
'v_xy', 'v_xz', 'v_xy_av', 'v_xz_av',
'w_xy', 'w_xz', 'w_xy_av', 'w_xz_av',
'p_xy', 'p_xz', 'p_xy_av', 'p_xz_av',
'theta_xy', 'theta_xz',
'wspeed', 'wdir',
'dust_mc_bin1', 'dust_mc_bin2', 'dust_mc_bin3',
'dust_mc_bin4', 'dust_mc_bin5',
'dust_mc_bin2_xy', 'dust_mc_bin2_xz', 'dust_mc_bin2_yz',
'dust_mc_bin2_av',
'dust_mc_bin3_xy', 'dust_mc_bin3_xz', 'dust_mc_bin3_yz',
'dust_mc_bin3_av',
'dust_emis_flux*_bin1_xy', 'dust_emis_flux*_bin2_xy',
'dust_emis_flux*_bin3_xy', 'dust_emis_flux*_bin4_xy',
'dust_emis_flux*_bin5_xy',
'dust_emis_flux*_bin1_xy_av',
'dust_depo_flux*_bin1_xy', 'dust_depo_flux*_bin2_xy',
'dust_depo_flux*_bin3_xy', 'dust_depo_flux*_bin4_xy',
'dust_depo_flux*_bin5_xy',
'dust_depo_flux*_bin1_xy_av',
'clay', 'clay_av', 'clay_xy', 'clay_xy_av',
'dust', 'dust_av', 'dust_xz', 'dust_xz_av',
'silt', 'silt_av', 'silt_yz', 'silt_yz_av',
data_output_pr = '#u', '#v', 'w', 'w_subs',
'wtheta', 'w"theta"', 'w*theta*',
'e', 'e*',
'u*2', 'v*2', 'w*2', 'theta*2',
'#theta', '#km', '#kh',
'hyp', 'p', 'rho',
'clay', 'silt', 'dust',
section_xy = 1, 5, 10, ! grid index for 2D XY cross sections
section_xz = 200, 1, ! grid index for 2D XZ cross sections
section_yz = 200, 1, ! grid index for 2D YZ cross sections
averaging_interval = 900.0, ! averaging interval general data
dt_averaging_input = 10.0, ! averaging general data sampling rate
averaging_interval_pr = 900.0, ! averaging interval profile data
dt_averaging_input_pr = 10.0, ! averaging profile data sampling rate
/ ! end of runtime parameter namelist
!
! DUST EMISSION AND TRANSPORT PARAMETER NAMELIST
!
&det_parameters
!
! general parameters
!
deposition_scheme = 'Z01', ! scheme for calculating deposition
det_start_time = 0.0, ! simulation time after det is active
switch_off_module = .FALSE., ! switch of (de)activating det
!
! parameters for dry deposition
!
alpha_imp = 50.0, ! land use dependent parameter for impaction efficiency calculation
brownian_diffusion_coefficient = 0.54, ! parameter for Brownian diffusion
!
! boundary conditions
!
bc_dm_b = 'neumann', ! bottom boundary condition for prognostic det variables
bc_dm_l = 'cyclic', ! west/left boundary condition for prognostic det variables
bc_dm_n = 'cyclic', ! north boundary condition for prognostic det variables
bc_dm_r = 'cyclic', ! east/right boundary condition for prognostic det variables
bc_dm_s = 'cyclic', ! south boundary condition for prognostic det variables
bc_dm_t = 'neumann', ! top boundary condition for prognostic det variables
!
! saltation bin parameters
!
bin_mass_fraction_ssc = 1.0, 0.25, 0.25, 0.25, 0.25, 0.0205, 0.0410,
0.0359, 0.3897, 0.5128, ! binspecific mass fraction of corresponding soil separate class
diameter_saltation = 1.42E6, 8.0E6, 20.0E6, 32.0E6, 44.0E6,
70.0E6, 130.0E6, 200.0E6, 620.0E6,
1500.0E6, ! effective diameter of a saltation size bin
mass_fraction_ssc = 0.03, 0.05, 0.05, 0.05, 0.05, 0.92, 0.92, 0.92,
0.92, 0.92, ! mass fraction of soil separate class
n_saltation_bins = 10, ! number of saltation size bins
particle_density_saltation = 2500.0, 2650.0, 2650.0, 2650.0, 2650.0, 2650.0,
2650.0, 2650.0, 2650.0, 2650.0, ! particle density of a saltation size bin
!
! dust bin parameters
!
diameter_dust = 1.46E6, 2.8E6, 4.8E6, 9.0E6, 16.0E6, ! effective diameter of a dust size bin
lower_bound_diameter = 0.2E6, 2.0E6, 3.6E6, 6.0E6, 12.0E6, ! minimum effective diameters represented by the dust size bin
n_dust_bins = 5, ! number of dust size bins
particle_density_dust = 2500.0, 2650.0, 2650.0, 2650.0, 2650.0, ! particle density of a dust size bin
upper_bound_diameter = 2.0E6, 3.6E6, 6.0E6, 12.0E6, 20.0E6, ! maximum effective diameters represented by the dust size bin
/ ! end of det parameter namelist
Note: DET can also be used with an empty namelist because the default configuration sets every namelist parameter appropriately. Output of time series of horizontally averaged bulk deposition flux and bulk emission flux (where bulk means the integral over all size bins), as well as maximum values of the bulk emission flux and minimum values of the bulk deposition flux, are added by default to the standard time series output.
Notes, shortcommings and open issues#

In contrast to Klamt et al. (2024), where the molecular viscosity of air has been calculated at each grid point and for each time step, a constant viscosity for an air temperature defined via namelist parameter pt_surface is used. Furthermore, the calculated settling velocity is also based on pt_surface, e.g. the height dependence of temperature is ignored.

It is not recommended to simulate sand particles (\(63\,\mu\)m \(< D \leq 2\,\)mm) with the provided implementation, especially not the larger fraction (\(\gg 100\,\mu\)m), because, on the one hand, the implemented emission algorithm is specially designed and tested for dustsized particles and on the other hand, the gravitational settling is calculated from Stoke's law, which can only be regarded as valid for particle diameters \(< 100\,\mu\)m.

DET is realized for horizontal upwardfacing default surfaces only. Land and urban surfaces are not considered. All soil parameters given in DET are assumed to be horizontally homogeneous.

You may use the LSM together with the DET, but other LSM surface types than flat sandy surfaces should be simulated with caution, especially rougher surfaces with plants. Also, heterogeneous LSM surfaces may not work appropriately. If you plan such simulations, check every setting of the parameters carefully. The equations of the deposition scheme might also required to be adjusted in such cases.

DET can not be used with topography.

The effects of soil moisture and wet deposition are not considered.

The configuration of the saltation size bins in Table 1 must start with at least one saltation size bin that is assigned to the soil separate class clay. Otherwise, the sandblasting efficiency, which connects the total horizontal saltation flux with the vertical bulk dust flux, can not be calculated (Marticorena and Bergametti, 1995).

Other, probably more accurate deposition schemes, should be implemented. Zhang and Shao (2014) and Zhang and He (2014) are two promising options. They were evaluated by Bergametti et al., 2018 and Khan and Perlinger, 2017.

So far not implemented:
 horizontal boundary conditions for the prognostic quantities of DET that are not cyclic,
 masked output.