An alternative method¶
dparcel.parcel.Parcel
functions well and produces reasonable results, but
the method it is built on has some inherent inefficiencies that would make it
unsuitable for applications requiring large amounts of calculation with low
computational cost, such as convection parametrisation in global climate
models.
dparcel.parcel.FastParcel
implements the far more elegant and efficient method
described by [sherwood_et_al_2013] in their Section 4 for homogeneous
parcels. There is very good agreement between dparcel.parcel.Parcel
and
dparcel.parcel.FastParcel
, and dparcel.parcel.FastParcel.motion
is
approximately twice as fast as dparcel.parcel.Parcel.motion
. The agreement is
best for small entrainment rates.
The major advantage of the new method is that the parcel’s properties at any height \(z\) can be found without needing to calculate them at every intermediate level between \(z_0\) and \(z\) in a stepwise fashion. This eliminates any accumulation of errors and allows us to achieve the same accuracy with fewer iterations. It is also not necessary to constantly check whether evaporation or condensation must occur in the parcel at each entrainment step.
The approach consists of the following steps.
Equivalent potential temperature¶
[sherwood_et_al_2013] assumed that entrainment mixes equivalent potential temperature between the parcel and its environment in the same way that it mixes temperature and moisture, as expressed in their Equation (6):
where \(\theta_e\) is the parcel equivalent potential temperature and \(\bar{\theta_e}\) is the corresponding environmental value at the same height. The approach then seeks to use the known \(\theta_e\) to determine the temperature at any height.
This works for updrafts, but we must modify it for downdrafts so that the \(\theta_e\) approaches \(\bar{\theta_e}\) regardless of the direction of motion from the initial height \(z_0\):
This is implemented in
dparcel.parcel.FastParcel.equivalent_potential_temperature
.
At this stage, we have not explicitly accounted for the effect of introducing
environmental air on the rate of evaporation in the parcel, but since
\(\theta_e\) is conserved both for adiabatic descent and evaporation at
constant pressure, any additional evaporation should not invalidate the result.
Total water content¶
Pressure and equivalent potential temperature alone do not define a unique temperature; we must also determine the specific humidity. As long as the parcel descends moist adiabatically, we will have \(q = q^*(p,T)\), enabling the calculation of \(T\). However, there may come a point where the saturation specific humidity becomes equal to the total amount of water (liquid and vapour) in the parcel (i.e., no liquid is left). Beyond this point the descent will be dry adiabatic.
In order to find the transition point between moist and dry descent, we reason that since the total amount of water \(Q\) in the parcel is conserved by adiabatic descent and evaporation without entrainment, it is mixed by entrainment in the same manner as \(\theta_e\):
Once we have found \(Q(z)\), we may find the point at which \(Q(z) = q^*(p,T)\) – this is the transition point. Beyond this point we will have \(q = Q(z) < q^*(p,T)\).
Temperature for moist descent¶
We begin by finding temperature as a function of height, assuming moist descent, by solving
where \(\theta_e^{\mathrm{sol}}(z)\) is the previously obtained parcel value. We use Newton’s method, with
We will first need to evaluate the saturation equivalent potential temperature \(\theta_e^*(p,T) = \theta_e(p, T, q^*(p,T))\) and its partial derivative with respect to temperature.
We will use the approximation presented by [bolton_1980]:
where \(r\) is the mixing ratio and
is the potential temperature at the LCL, where \(T_K\) is the absolute temperature and
is the temperature at the LCL.
For a saturated parcel, we can see that \(T_K = T_D = T_L\). Thus
and
We first note that \(\partial \theta_E / \partial T' = \theta_E (\partial \log \theta_E / \partial T_K)\), and
and
with
and, lastly,
where \(a = 17.67\), \(b = 243.5\) K, \(e_0 = 6.112\) mbar and \(C = 273.15\) K.
We implement the calculation in
dparcel.thermo.saturation_equivalent_potential_temperature
.
We then continue with the Newton’s method solution, using the approximation of
[davies-jones_2008] for the case of non-entraining moist pseudoadiabatic
descent as a first guess. This is implemented in
dparcel.parcel.FastParcel.properties_moist
.
Temperature for dry descent¶
We may perform a similar computation, assuming dry descent, now solving
with Newton’s method giving
We must first find \(\frac{\partial}{\partial T} [\theta_e(p, T, Q(z))]\). Returning to the approximation of [bolton_1980], now with specific humidity independent of temperature, we have
Now,
and
Bolton’s (10) for the saturation vapour pressure \(e_s\) implies that
with \(U\) being the relative humidity. Then
Bolton’s (10) in its original form gives
which implies
Finally, we can also substitute the exact relation
where \(\epsilon\) is the ratio of the molar mass of dry air to the molar mass of water vapour.
We use the non-entraining dry adiabatic value as a first guess.
This step is implemented in
dparcel.parcel.FastParcel.properties_dry
.
Transition between moist and dry descent¶
We then find the transition point where \(Q(z) = q^*(p,T)\)
(or equivalently \(l(z) = Q(z) - q^*(p,T) = 0\) where \(l\) is the
liquid water mass ratio). This is implemented in
dparcel.parcel.FastParcel.transition_point
.
Temperature for general descent¶
At last, we may combine moist and dry descent to find the final temperature,
specific humidity and liquid ratio as a function of height. This is implemented
in dparcel.parcel.FastParcel.properties
.
Buoyancy¶
It is then a simple matter to find the buoyancy:
where \(T_v\) and \(\bar{T_v}\) are the parcel and environmental
virtual temperatures. This is implemented in
dparcel.parcel.FastParcel.buoyancy
.
Motion¶
Knowing the buoyancy as a function of height, we simply substitute the new
buoyancy function into the existing code from
dparcel.parcel.Parcel.motion
to simulating the parcel’s motion.
This is implemented in dparcel.parcel.FastParcel.motion
.