Science:Saltation model

From SwissExperiment

Saltation model

Background

A saltation model was written by Judith Doorschot (SLF) and Michael Lehning (SLF) in 2001-2004 as part of J. Doorschot's Ph.D. <ref name=doorschot2001a/> The original model was written in C, and validated using data collected at the Wiessfluhjoch Versuchsfeld in Davos.

The model is a combination of a lagrangian particle motion model and a logarithmic velocity boundary layer. As part of the modelling process, the boundary layer profile is modified to account for reduced shear stress near the surface. This reduction of shear stress corresponds to momentum transfer to the surface, required to entrain particles.

The model was ported into Matlab by Andy Clifton (SLF) as part of his Ph.D. work in 2006-2008, with the assistance of Katherine Leonard (then LDEO at Columbia, later SLF). Several aspects of the model physics - such as threshold velocity and the parameters describing the entrainment process - were updated with new experiment data, resulting in a code that can be applied to a wider range of situations than those used by Doorschot.

While not a product of SwissEx:Science, this code may be needed as a part of the project. The wiki also serves as a way to distribute and document the code.

Model Physics

The model calculates the rate of particle entrainment from the surface and the trajectories followed by particles from the surface and through the boundary layer. To do this, it requires models for the entrainment process, particle flight and air boundary layer, and how the feedback between saltation and the boundary layer profile occurs. These are shown in the diagrams below.

In the model, particles are assumed to be uniformly distributed on a surface over which a fluid is flowing. The velocity boundary layer in the fluid (typically air) is defined using the roughness length and friction velocity as observed above the saltion layer. Particles are entrained from the surface and move through the boundary layer following generally parabolic paths, which are modified by lift and drag.

In order to improve the accuracy of the model, the surface grain size, ejection angle and ejection velocity are assumed to follow log-normal distributions, and every combination of 'classes' is simulated before being combined to give a total horizontal mass flux. The horizontal mass fluxes occuring at different sizes and ejection parameters are combined according to the frequency of occurence, so that the mass flux of a case corresponding to 10% of particle diameters, 12% of ejection velocities, 8% of ejection angles and 5% of ejection speeds is weighted to be 0.1 x 0.12 x 0.08 x 0.05 = 0.000048 or 0.048% of the 'total' mass flux.

DCLL saltation model overview	DCLL saltation model mass flux loop

There are slight differences in the physics included in the model used by Doorschot and Lehning (the 'DL' model), and that developed by Clifton, Lehning and Leonard and from the original version (the 'DCLL' model). These differences include;

Parameter	Symbol	DL	DCLL
Threshold velocity	$u * t$	user specified (from measurements)	$u_{*t} = A \sqrt{\left(\frac{\rho_{ice}-\rho_{air}}{\rho_{air}}\right)gd_p},$ where $A = 0.18$ for alpine (cohesive) snow.
Mean ejection velocity	$v e j$	$2.5 u *$	$u_\left[3.5-\left(1-e^{-5\left(\frac{u_}{u_{*t}}-1\right)}\right)\right]$
Standard deviation of ejection velocity	$\sigma\left(v_{ej}/ u_*\right)$	$0.75$	$0.5$
Mean ejection angle	$α e j$	$25^\circ$	$75-55\left(1+e^{-d/175}\right)$
Standard deviation of ejection angle	$\sigma\left(\alpha_{ej}\right)$	$15^\circ$	$\frac{1}{2}\overline{\alpha_{ej}}$
Roughness length	$z 0 m$	user specified (from profile measurements) or $d / 10$	$0.021\frac{u_{*}^2}{2g}$
Entrainment rate	$ζ a e$	$\frac{1}{8\pi {d_p}^{2}}$	$\frac{\rho_{snow}}{\rho_{ice}}\times\frac{1}{\pi {d_p}^{2}/4}$
Surface coefficient of restitution	$r$	0.35	$\frac{\rho_{snow}}{\rho_{ice}}$

As well as this, the matlab code developed for the DCLL model allows the vertical flux profile to be resolved and the mass flux and particle size distribution through a particular 'slice' of the profile - equivalent to deploying a particle trap or snow particle counter ('SPC') in the flow. This required considerable alterations to the code used by DL, although the two codes deliver the same final mass flux, given the same inputs.

MATLAB Code

The following section describes the code deposited in the wiki as File:DCLL saltation standalone.zip. For details of the contents, see Matlab File Package (below).

Concept

The code is split into several parts. These prepare the simulation, set global variables (flags) and call the main simulation function. This is done through 3 main files;

run_DCLL_saltation_standalone.m; this is a script that is used to specify the current case to simulate, call the function to set the global variables, run the simulations and plot the results. This could be extended to include saving the input variables and output results to file, if required. Values for $z 0 m$ and $u * t$ are calculated in this script from other inputs, before being given as inputs to the main function. DCLL_saltation_standalone_GLOBALS.m; this function sets the global values. These values are used by many different parts of the saltation model, and care should be exercised when changing some values. There are also options such as the vertical resolution to use for the mass flux profile and the height of the top and bottom of the 'trap' to extract particle size distribution information, at. DCLL_saltation_standalone_main.m; this is the actual main function of the saltation model. Any changes made here may alter the model function and results. Changes (corresponding to a particular modelling case) should be limited to run_DCLL_saltation_standalone.m and DCLL_saltation_standalone_GLOBALS.m	Matlab DCLL saltation model concept

Setting up a simulation

A simulation is defined by environmental variables that describe things like the air temperature and pressure and surface density and grain size. These variables are defined in run_DCLL_saltation_standalone.m and are used to define the input arguments to the saltation model function, DCLL_saltation_standalone_main.

The syntax to call the saltation model from within matlab is flux_out = DCLL_saltation_standalone_main( tempS, pressS, rhoS, angleS, u_star, z0m, grain_size, grain_sdev, n_tau, n_dg, n_a, n_v) , where values in brackets on the right are input arguments. The output is on the left, flux_out.

The input arguments to the saltation model DCLL_saltation_standalone_main are;

Variable	Description	Units	Note
tempS	surface layer temperature	[K]
pressS	surface pressure	[Pa]
rhoS	surface density	[Kg/m3]
angleS	slope angle	[radians]
u_star	friction velocity	[m/s]
u_star_th	threshold friction velocity	[m/s]	calculated in run_DCLL_saltation_standalone.m using $u_{*t} = A \sqrt{\dfrac{\rho_p-\rho_f}{\rho_f}gd_p}$
z0m	roughness length	[m]	calculated in run_DCLL_saltation_standalone.m using $z_{0m} = 0.021\dfrac{{u_*}^2}{2g}$
grain_size	mean grain size	[m]
grain_sdev	standard deviation of grain size	[m]
n_tau	number of friction velocity classes to consider (integer)	[-]	strongly recommend n_tau=1
n_dg	number of grain classes to consider (integer)	[-]	strongly recommend n_dg>=11
n_a	number of angle classes to consider (integer)	[-]	strongly recommend n_a>=11
n_v	number of velocity classes to consider (integer)	[-]	strongly recommend n_v>=11

The following global variables are hard-coded within DCLL_saltation_standalone_GLOBALS.m. The use of global variables allows them to be manipulated both within the code and by scripts that are used to call the saltation model.

Influence	Variable	Description	Default
Material Properties	RHO_ICE	Density of ice	917 kg/m3
	RHO_P	Particle density	RHO_ICE
	BAGNOLD_A	$A$ in the threshold friction velocity equation	0.18 for cohesive snow 0.1 for non-cohesive snow
Local Environment	ANGLE_SLOPE	Slope angle [radians]. Inherited from input variable angleS.	0
	GRAV	acceleration due to gravity	9.81
	KAPPA	von karman constant	0.41
Numerical Scheme	DT_SALT	time step to use in the simulations	0.0005
	DZ	Vertical resolution to use in the simulations	0.00002
	CONV_TEST	Maximum allowed change in mass flux between iterations.	0.001 (i.e. 0.1%)
Particle tracking	ON_TRAP	Enables the trap. 0 = trap off, 1 = trap on.	0
	ON_PROFILE	Generates a mass-flux profile	0 (automatically enabled if ON_TRAP = 1)
	SPCZMAX	Upper edge of the trap [m].	N/A
	SPCZMIN	Lower edge of the trap [m].	N/A
	SPCETA	Trap efficiency [0-1.00].	1
Suspension	ON_SUSPLINK	Enables the suspension monitoring. 0 = no monitoring, 1 = monitoring enabled. This part of the code has not been validated and is not recommended for use	0
Debugging	DEBUG_DIST	Display information about the particle size probability distribution function and compares it to data from Budd (1966)	0
	DEBUG_TRAJECT	Display information about the particle trajectory.	0
	DEBUG_PROFILE	Display information about the mass flux profile	0
	DEBUG_TRAP	Display information about the 'trap'	0
	DEBUG_SUSPLINK		0
	DEBUG_AEROENTRAIN		0
	DEBUG_VW	Display information about particle motion	0

N.B. re ON_SUSPLINK and DEBUG_SUSPLINK. A potential link to a suspension code was implemented by assuming that any particle exceeding a certain height could become suspendend. This approach has not been validated and is not recommended for use. However, the function susplink.m could be modified to use any approach.

Outputs

The function DCLL_saltation_standalone outputs a single structure flux_out; flux_out = DCLL_saltation_standalone(...

This structure includes several elements. They can be accessed through flux_out.<element>, where <element> could be one of the following list;

z_profile; vertical position (array)
m_tot_profile; horizontal mass flux associated with the vertical position z_profile
THIS LIST NOT COMPLETE

Matlab file package

The following matlab files and functions are in the DCLL_saltation_standalone directory included in File:DCLL saltation standalone.zip;

File	Description
DCLL_saltation_standalone_GLOBALS.m	describes and defines the globals used in the model
DCLL_saltation_standalone_main
get_rho_air	calculates the density of air from the gas equation, considering relative humidity
get_bagnold_threshold.m	calculates the drift threshold friction velocity from the particle diameter, ice density and air density
run_DCLL_saltation_standalone.m	sets variables used in the simulation, calls the model and demonstrates post-processing

The following functions are in the private directory;

File	Description
get_gamma_dis.m	Generates a gamma distribution of particle sizes from an input mean grain diameter and standard deviation.
aero_entrain.m	Calculates the mean velocity, rise height and impact and ejection angle of a saltating particle
budd_snowdist.m	returns the surface snow size distribution observed by Budd (1966, table 5)
dynvisc
get_dvisc_air
mass_flux	calculates the mass leaving the ground
get_lognorm_dis	generates a log-normal distribution using the same approach as in the original DL model
prob_profile	calculates the probability of a particle passing through a given height
saltation_DCLL_standalone_mflux_loop	calculates the horizontal mass flux for one particular combination of parameters.
SQR	calculates the square root of a number. Legacy code to match a macro in the DL model.
susplink
test_saltation	checks for the presence of 'strong' (path lengths increasing over time) or 'weak' (path lengths decreasing) saltation
traject	calculates a particle trajectory, including lift and drag
vw	calculates the wind velocity profile from the friction velocity and roughness length
vw2	calculates a modified wind velocity profile using the saltation height

Potential modifications

General

Matlab code

The matlab code is designed to allow complete debugging and data display. Consequently it is slow and needs large amounts of memory. A faster, leaner 'operational' version could be written
The function vw.m requires an estimate of the saltation layer height; this is calculated using a not-entirely clear expression.

References

Doorschot J, Lehning M. 2002. Equilibrium saltation: mass fluxes, aerodynamic entrainment, and dependence on grain properties. Boundary-Layer Meteorology 104(1): 111-130. DOI: 10.1023/A:1015516420286.
Doorschot J, Lehning M, Vrouwe A. 2004. Field measurements of snow drift and mass fluxes and comparison with model simulations. Boundary-Layer Meteorology 113: 347-368. DOI: 10.1007/s10546-004-8659-z.
Clifton A, Lehning M. 2008. Improvement and validation of a snow saltation model using wind tunnel measurements, 33(14): 2156-2173. DOI: 10.1002/esp.1673.

Swiss Experiment

Project Information

SwissEx Science

Work Packages

Products

Data

Data Input/Output