|
|
|
Documentation:
|
|
|
|
Below-cloud Interaction Model (BCIM)
|
|
|
|
A model to simulate the isotopic evolution of a single hydrometeor falling from the cloud through a vertical air column to the ground
|
|
|
|
|
|
|
|
Developed and tested by Pascal Graf, 2016 - 2018
|
|
|
|
pascal.graf@env.ethz.ch
|
|
|
|
|
|
|
|
Associated MATLAB functions
|
|
|
|
Main functions:
|
|
|
|
|
|
|
|
Function name
|
|
|
|
Purpose
|
|
|
|
below_cloud_model.m
|
|
|
|
Main model. Calculates the phyiscal and isotopic properties of a single falling hydrometeor from the given inputs (see “input”) and returns an output struct with fields listed under “output”.
|
|
|
|
calc_isotopes.m
|
|
|
|
Diagnoses the isotopic composition of a hydrometeor along its fall trajectory. The mass of the hydrometeor and the isotopic composition of the ambient air at every height step have to be known.
|
|
|
|
find_input_diameter.m
|
|
|
|
Iteratively determines the initial diameter of a hydrometer, which reaches a target terminal diameter when it arrives on the ground. This function is needed when the terminal instead of the default initial diameter is given as model input.
|
|
|
|
get_profile_value.m
|
|
|
|
Get the interpolated value of a profile variable at a specified height.
|
|
|
|
rayleigh_ascent.m
|
|
|
|
Calculates vertical profiles of T, h, δ2H and d-excess of a (moist-)adiabatically ascending air parcel, which undergoes Rayleigh fractionation.
|
|
|
|
|
|
|
|
Helper functions:
|
|
|
|
|
|
|
|
Function name
|
|
|
|
|
|
|
|
Purpose
|
|
|
|
air_density.m
|
|
|
|
|
|
|
|
Calculates the air density ρ.
|
|
|
|
alphaH.m; alphaO.m
|
|
|
|
|
|
|
|
Calculates the liquid-vapour equilibrium fractionation factor (for 2H or 18O)
|
|
|
|
alphaHk.m; alphaOk.m
|
|
|
|
|
|
|
|
Calculates the kinetic fractionation factor (for 2H or 18O).
|
|
|
|
alphaHs.m; alphaOs.m
|
|
|
|
|
|
|
|
Calculates the solid-vapour equilibrium fractionation factor (for 2H or 18O).
|
|
|
|
diffusion_coefficient.m
|
|
|
|
|
|
|
|
Calculates the diffusion coefficient.
|
|
|
|
drop_mass_tendency.m
|
|
|
|
|
|
|
|
Calculates the hydrometeor mass change due to evaporation and condensation.
|
|
|
|
drop_T_tendency.m
|
|
|
|
|
|
|
|
Calculates the hydrometeor temperature change due to evaporation, condensation and heat diffusion.
|
|
|
|
dynamic_viscosity.m
|
|
|
|
|
|
|
|
Calculates the dynamic viscosity.
|
|
|
|
heat_vent_coeff.m
|
|
|
|
|
|
|
|
Calculates the heat ventilation coefficient.
|
|
|
|
latheat_evap.m
|
|
|
|
|
|
|
|
Calculates the latent heat of evaporation.
|
|
|
|
latheat_sublim.m
|
|
|
|
|
|
|
|
Calculates the latent heat of sublimation.
|
|
|
|
mass_mix_ratio.m
|
|
|
|
|
|
|
|
Calculates the mass mixing ration of vapour in dry air.
|
|
|
|
mass_vent_coeff.m
|
|
|
|
|
|
|
|
Calculates the mass ventilation coefficient (possible for each isotope species)
|
|
|
|
satvappress_ice.m
|
|
|
|
|
|
|
|
Calculates the saturation vapour pressure over ice.
|
|
|
|
satvappress.m
|
|
|
|
|
|
|
|
Calculates the saturation vapour pressure over liquid.
|
|
|
|
terminal_velocity_solid.m
|
|
|
|
|
|
|
|
Calculates the terminal velocity of a solid hydrometeor.
|
|
|
|
terminal_velocity.m
|
|
|
|
|
|
|
|
Calculates the terminal velocity of a liquid hydrometeor.
|
|
|
|
therm_cond.m
|
|
|
|
|
|
|
|
Calculates the thermal conductivity of moist air.
|
|
|
|
|
|
|
|
Inputs:
|
|
|
|
Input variable
|
|
|
|
Description
|
|
|
|
diam_input
|
|
|
|
Cell array (1x2), containing the equivalent diameter of the hydrometeor in [m] and an indication whether this is the diameter at the starting height 'start' or the diameter when arriving at the surface 'end'. The diameter at the starting height is iteratively determined in the latter case.
|
|
|
|
Example: diam_input={0.001,'end'} → Hydrometeor arrives on the ground with a diameter of 1 mm.
|
|
|
|
prof_deltaH
|
|
|
|
(1xn) array, containing the vertical delta2H profile of environmental water vapour. If n=2, the profile is linearly interpolated between the lowest and highest model level (defined with the input variable z). If n=length(z), the entire profile of vapour delta2H is defined. Values of prof_deltaH should be given in [permil]. Dimension has to be (1x2) if namelist.profilecalc=1. Then, only the first value will be used as starting value for the Rayleigh ascent.
|
|
|
|
Example: prof_deltaH=[-50 -300] → Linear delta2H profile with -50 permil at the surface (z(1)) and -300 permil at the highest model level (z(end)).
|
|
|
|
prof_dexc
|
|
|
|
Same as prof_deltaH, but for deuterium-excess in [permil]
|
|
|
|
prof_rh
|
|
|
|
Same as prof_deltaH, but for relative humidity in [] (e.g., 50% = 0.5)
|
|
|
|
prof_T
|
|
|
|
Same as prof_deltaH, but for air temperature in [°C.]
|
|
|
|
z
|
|
|
|
(1xn) array, containing the vertical height levels in [m], on which the above arrays are defined and on which the output (profiles) will be defined.
|
|
|
|
Example: z=linspace(500,3500,3001) → Linear array of height levels from 500 m to 3500 m with steps of 1 m: [500, 501, 502,... ,3500]
|
|
|
|
namelist
|
|
|
|
Structure of different model input parameters in the format "namelist.fieldname". See below for possible fieldnames.
|
|
|
|
|
|
|
|
The following list shows possible field names for namelist. If they are not defined, the default (highlighted in grey) values are used. The namelist must contain at least one field in order to run below_cloud_model.m.
|
|
|
|
|
|
|
|
Namelist field
|
|
|
|
Value
|
|
|
|
Description
|
|
|
|
drop_T
|
|
|
|
1
|
|
|
|
Allow drop temperature to be different from air temperature.
|
|
|
|
0
|
|
|
|
Drop temperature equals air temperature
|
|
|
|
p0
|
|
|
|
100'000
|
|
|
|
Pressure at the ground (z(1)) in [Pa]
Here: 100'000 [Pa] = 1000 [hPa]
|
|
|
|
massincrease
|
|
|
|
1
|
|
|
|
Allow drop mass increase.
|
|
|
|
0
|
|
|
|
Only allow drop mass decrease.
|
|
|
|
dt
|
|
|
|
0.1
|
|
|
|
Time integration time step [s].
|
|
|
|
0
|
|
|
|
Time integration disabled. Integrate after dz. dz is defined in its own namelist field.
|
|
|
|
profilecalc
|
|
|
|
1
|
|
|
|
Calculate profiles of T, rh, delta2H and dexc with a Rayleight ascent. Only the first (lowest) value of prof_T, prof_rh, prof_delta2H and prof_dexc is taken to initiate the ascent of the air parcel from the ground. The rest of the profiles are calculated with the MATLAB function "rayleigh_ascent.m", by (moist-) adiabatically lifting the air parcel. All condensed moisture is rained out immediately. prof_delta2H needs to be an array with size (1x2) if namelist.profilecalc=1.
|
|
|
|
0
|
|
|
|
Use input profiles of T, rh, delta2H and dexc. If input profiles are of size (1x2), they are linearly interpolated between z(1) and z(end).
|
|
|
|
solidterminalvelocity
|
|
|
|
0
|
|
|
|
Use terminal velocity of a liquid drop for both liquid and solid hydrometeors.
|
|
|
|
1
|
|
|
|
Use a separate terminal velocity for solid hydrometeors.
|
|
|
|
fract_T
|
|
|
|
'drop'
|
|
|
|
Drop temperature is used for the fractionation equation.
|
|
|
|
'env'
|
|
|
|
Environmental temperature is used for the fractionation equation.
|
|
|
|
'mix'
|
|
|
|
Average between drop and environmental temperature is used for the fractionation equation.
|
|
|
|
z_form
|
|
|
|
0
|
|
|
|
Formation and starting height of the hydrometeor is at the highest model level (z(end)).
|
|
|
|
3000
|
|
|
|
starting height is at 3000 m above sea level (values between z(1) and z(end) are possible).
|
|
|
|
-1
|
|
|
|
starting height is at the lowermost level where the air is supersaturated with respect to ice (exm > 0).
|
|
|
|
formmech
|
|
|
|
'wbf'
|
|
|
|
Hydrometeor is formed by the Wegener-Bergeron-Findeisen process.
|
|
|
|
'depo'
|
|
|
|
Hydrometeor is formed by vapour deposition on ice.
|
|
|
|
'riming'
|
|
|
|
Hydrometeor is formed by riming of liquid droplets on ice. Corresponds to condensation of vapour on liquid. The degree of riming can be defined with namelist.rimingfrac (see below).
|
|
|
|
rimingfrac
|
|
|
|
0.5
|
|
|
|
Defines the degree (or fraction) of riming if namelist.formmech='riming'. Here: 50% of the hydrometeor mass is formed by riming, the remaining 50% are formed by the WBF process. Values between 0 and 1 are possible.
|
|
|
|
dz
|
|
|
|
1
|
|
|
|
Defines step size of the vertical integration and of the output in [m] if namelist.dt=0.
|
|
|
|
Here: Vertical resolution of integration and output is 1 m.
|
|
|
|
isocalc
|
|
|
|
'explicit'
|
|
|
|
Isotopic composition of the hydrometeor is calculated using separate mass equations for all isotopic species.
|
|
|
|
'stewart'
|
|
|
|
Isotopic composition of the hydrometeor is calculated according to Stewart (1975) (Eq. 2).
|
|
|
|
temptendswitch
|
|
|
|
'both'
|
|
|
|
Drop temperature changes due to evaporative cooling of adaptation to the environmental air temperature are both switched on.
|
|
|
|
'evap'
|
|
|
|
Only evaporative cooling is enabled.
|
|
|
|
'adap'
|
|
|
|
Only adaptation to environmental air temperature is enabled.
|
|
|
|
constisotopeprofile
|
|
|
|
0
|
|
|
|
Take the isotope profile of vapour from the input (namelist.profilecalc=0) or calculate it with a Rayleigh ascent (namelist.profilecalc=1).
|
|
|
|
1
|
|
|
|
Define the isotope profile of vapour by a Rayleigh ascent with starting values T=12 and rh=0.75.
|
|
|
|
This can be useful to test the sensitivity to temperature or humidity profiles, but keeping the isotope profile of vapour constant.
|
|
|
|
fulloutput
|
|
|
|
0
|
|
|
|
Subfunction "calc_isotopes" returns profiles of deltaHp, deltaOp, deltaHp_eq, deltaOp_eq, f and ftot.
|
|
|
|
1
|
|
|
|
Subfunction "calc_isotopes.m" additionally returns "mvaptot" (total evaporated vapour mass at each level), "RHvap" and "ROvap" (isotope ratios of the evaporated vapour).
|
|
|
|
soundingprofile
|
|
|
|
0
|
|
|
|
Input profiles of T and rh are not from an atmospheric sounding or are not be treated like one. Rayleigh ascent (if namelist.profilecalc=1) uses the values of the lowest index (e.g. T(1)) as starting values.
|
|
|
|
1
|
|
|
|
Take the average of the lowermost 100 m of the profiles of T and rh as starting values for the Rayleigh ascent. This should be switched on when taking profiles of T and rh from atmospheric soundings, where surface effects can cause starting values which are very different from the air layer above and thus potentially cause unrealistic profiles.
|
|
|
|
|
|
|
|
|
|
|
|
Outputs:
|
|
|
|
below_cloud_model.m returns a struct "profiles" with the arrays listed below. The array sizes depends on the model inputs.
|
|
|
|
|
|
|
|
Variables that describe the profile of ambient air are of size (1xn), where n is the length of the input array "z". "z" defines the grid on which these profiles are defined (background grid).
|
|
|
|
For example:
|
|
|
|
z: (1x3001) array; p: (1x3001) array
|
|
|
|
z(305) = 805 m
|
|
|
|
p(305) = 93000 Pa
|
|
|
|
→ The pressure at model level 305 is 930 hPa, which corresponds to 805 m above sea level.
|
|
|
|
Variables that describe properties of the hydrometeor are arrays of size (1xm), where m is the number of model levels, which are at and below the starting altitude (z_form) of the hydrometeor. The corresponding grid (hydrometeor grid) is defined in "z_out" (and "height")
|
|
|
|
For example:
|
|
|
|
z = 0:1:3000 -> (1x3001) array
|
|
|
|
z_form = 2000
|
|
|
|
z_out contains all z<=z_form:
|
|
|
|
z_out = 0:1:2000 -> (1x2001) array
|
|
|
|
m=n if z_form=z(end) (if namelist.z_form=0)
|
|
|
|
|
|
|
|
Array name
|
|
|
|
Size
|
|
|
|
Description
|
|
|
|
z
|
|
|
|
(1xn)
|
|
|
|
Height above sea level of all model levels [m]
|
|
|
|
z_out
|
|
|
|
(1xm)
|
|
|
|
Height above sea level of all model levels at and below the formation height of the hydrometeor [m]
|
|
|
|
height
|
|
|
|
(1xm)
|
|
|
|
Same as z_out.
|
|
|
|
h
|
|
|
|
(1xn)
|
|
|
|
Relative humidity of ambient air []
|
|
|
|
T
|
|
|
|
(1xn)
|
|
|
|
Temperature of ambient air [°C]
|
|
|
|
p
|
|
|
|
(1xn)
|
|
|
|
Pressure of ambient air [Pa]
|
|
|
|
exm
|
|
|
|
(1xn)
|
|
|
|
Excess moisture of ambient air (specific humidity divided by the saturation specific humidity over ice at ambient temperature; exm>=0) [kg/kg]
|
|
|
|
m
|
|
|
|
(1xm)
|
|
|
|
Hydrometeor mass [kg]
|
|
|
|
diam
|
|
|
|
(1xm)
|
|
|
|
Hydrometeor diameter [m]
|
|
|
|
Td
|
|
|
|
(1xm)
|
|
|
|
Hydrometeor temperature [°C]
|
|
|
|
Tdelta
|
|
|
|
(1xm)
|
|
|
|
Temperature difference between hydrometeor and surrounding air [°C]
|
|
|
|
equih
|
|
|
|
(1xm)
|
|
|
|
Relative humidity of the surrounding air with respect to Td, instead of T (equivalent relative humidity) []
|
|
|
|
masstend
|
|
|
|
(1xm)
|
|
|
|
Hydrometeor mass tendency (rate of change of mass dm/dt) [kg/s]
|
|
|
|
temptend
|
|
|
|
(1xm)
|
|
|
|
Hydrometeor temperature tendency (rate of change of temperature dTd/dt) [°C/s]
|
|
|
|
f
|
|
|
|
(1xm)
|
|
|
|
Remaining fraction of mass of the hydrometeor relative to the next higher model level []
|
|
|
|
ftot
|
|
|
|
(1xm)
|
|
|
|
Remaining fraction of mass of the hydrometeor relative to the initial mass (at z_form) []
|
|
|
|
evapfrac
|
|
|
|
(1xm)
|
|
|
|
Evaporated fraction of mass (current mass relative to the initial mass of the hydrometeor; evapfrac=1-ftot) []
|
|
|
|
deltaHv
|
|
|
|
(1xn)
|
|
|
|
delta2H value of the ambient vapour [permil]
|
|
|
|
deltaOv
|
|
|
|
(1xn)
|
|
|
|
delta18O value of the ambient vapour [permil]
|
|
|
|
dexcv
|
|
|
|
(1xn)
|
|
|
|
d-excess value of the ambient vapour [permil]
|
|
|
|
deltaHp
|
|
|
|
(1xm)
|
|
|
|
delta2H value of the hydrometeor [permil]
|
|
|
|
dexcp
|
|
|
|
(1xm)
|
|
|
|
d-excess value of the hydrometeor [permil]
|
|
|
|
deltaHp_eq
|
|
|
|
(1xm)
|
|
|
|
delta2H value of equilibrium vapour from the hydrometeor at ambient temperature [permil]
|
|
|
|
dexcp_eq
|
|
|
|
(1xm)
|
|
|
|
d-excess value of equilibrium vapour from the hydrometeor at ambient temperature [permil]
|
|
|
|
z_intsteps
|
|
|
|
(1xk)
|
|
|
|
Height above sea level of the individual integration steps [m] (z_intsteps defines the integration grid and is different from z if time integration was chosen). k = number of integration steps.
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Examples:
|
|
|
|
A hydrometeor falling from 3500 meters above sea level to the ground at 500 m a.s.l. (with surface pressure of 950 hPa). The hydrometeor falls through background air, which is defined by a (moist-)adiabatically ascending air parcel from the surface with starting values of T=12°C., h=75%, delta2H=-150 permil and d-excess=10 permil. Default settings are used if not defined by the namelist (see above), i.e., integration after time, mass increase allowed, explicit treatment of all isotope mass equations, wegener-bergeron-findeisen formation mechanism, liquid terminal velocity, drop temperature different from that of the environment, which can be modified by both evaporation and heat diffusion and which is used to calculate the fractionation factor.
|
|
|
|
|
|
|
|
namelist.p0=95000 % Surface pressure 950hPa
|
|
|
|
namelist.z_form=3500; % Formation height at 3500 m above sea level
|
|
|
|
z=linspace(500,3500,3000) % 3000 height steps from 500 to 3500 m
|
|
|
|
[T,h,~,~,~]=rayleigh_ascent(12,0.75,namelist.p0,-150,10,z);
|
|
|
|
% Profiles of T and h are defined by an air parcel that ascends % (moist-)adiabatically from 950hPa with starting temperature of 12°C.
|
|
|
|
% and starting relative humidity of 75%. Isotopic values (-150 and 10) % have to be given to rayleigh_ascent() as input, but are not of
|
|
|
|
% relevance here.
|
|
|
|
diam=1; % Hydrometeor diameter is 1mm
|
|
|
|
[profiles]=below_cloud_model({diam,’end’},[-150 0],[10 0],h,T,z,namelist) |
