Adding time integrator version of non-equilibrium formulation

Project.toml

@@ -1,7 +1,7 @@
 name = "CloudMicrophysics"
 uuid = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b"
 authors = ["Climate Modeling Alliance"]
-version = "0.22.1"
+version = "0.22.3"
 
 [deps]
 ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c"

docs/bibliography.bib

@@ -87,6 +87,17 @@ @article{Bigg1953
   doi = {10.1088/0370-1301/66/8/309}
 }
 
+@article{Blahak2010,
+  title = {Efficient approximation of the incomplete gamma function for use in cloud model applications},
+  author = {Blahak, U.},
+  journal = {Geoscientific Model Development},
+  volume = {3},
+  year = {2010},
+  number = {2},
+  pages = {329--336},
+  doi = {10.5194/gmd-3-329-2010}
+}
+
 @article{BrownFrancis1995,
   title = {Improved Measurements of the Ice Water Content in Cirrus Using a Total-Water Probe},
   author = {Philip R. A. Brown and Peter N. Francis},
@@ -120,6 +131,16 @@ @article{China2017
   year = {2017}
 }
 
+@article{Choletteetal2019,
+  author = {Mélissa Cholette et al},
+  title = {Parameterization of the Bulk Liquid Fraction on Mixed-Phase Particles in the Predicted Particle Properties (P3) Scheme: Description and Idealized Simulations},
+  journal = {Journal of the Atmospheric Sciences},
+  year = {2019},
+  volume = {76},
+  doi = {10.1175/JAS-D-18-0278.1},
+  pages= {561--582}
+}
+
 @article{Desai2019,
   title = {Aerosol-Mediated Glaciation of Mixed-Phase Clouds: Steady-State Laboratory Measurements},
   author = {Desai, Neel and Chandrakar, KK and Kinney, G and Cantrell, W and Shaw, RA},
@@ -616,6 +637,17 @@ @article{SeifertBeheng2006
   publisher = {Springer}
 }
 
+@article{Simmeletal2002,
+  title = {Numerical solution of the stochastic collection equation—comparison of the Linear Discrete Method with other methods},
+  author = {Simmel, Martin and Trautmann, Thomas and Tetzlaff, Gerd},
+  journal = {Atmoshperic Research},
+  volume = {61},
+  number = {2},
+  year = {2002},
+  doi = {10.1016/S0169-8095(01)00131-4},
+  publisher = {Elsevier}
+}
+
 @article{Spichtinger2023,
   AUTHOR = {Spichtinger, P. and Marschalik, P. and Baumgartner, M.},
   TITLE = {Impact of formulations of the homogeneous nucleation rate on ice nucleation events in cirrus},

docs/src/API.md

@@ -10,6 +10,7 @@ MicrophysicsNonEq
 MicrophysicsNonEq.τ_relax
 MicrophysicsNonEq.conv_q_vap_to_q_liq_ice
 MicrophysicsNonEq.conv_q_vap_to_q_liq_ice_MM2015
+MicrophysicsNonEq.conv_q_vap_to_q_liq_ice_MM2015_timeintegrator
 ```
 
 # 0-moment precipitation microphysics
@@ -69,6 +70,7 @@ P3Scheme.thresholds
 P3Scheme.distribution_parameter_solver
 P3Scheme.ice_terminal_velocity
 P3Scheme.het_ice_nucleation
+P3Scheme.ice_melt
 ```
 
 # Aerosol model

docs/src/MicrophysicsNonEq.md

@@ -93,25 +93,145 @@ where:
 - ``R_v`` is the gas constant of water vapor,
 - ``L_v`` and ``L_s`` is the latent heat of vaporization and sublimation.
 
-    Note that these forms of condensation/sublimation and deposition/sublimation are equivalent to those described in the adiabatic parcel model with some rearrangements and assumptions. It is just necessary to use the definitions of ```\tau```, ```q_{sl}```, and the thermal diffusivity ```D_v```:
-
-    ```math
-    \begin{equation}
-      \tau = 4 \pi N_{tot} \bar{r}, \;\;\;\;\;\;\;\;
-      q_{sl} = \frac{e_{sl}}{\rho R_v T}, \;\;\;\;\;\;\;\;
-      D_v = \frac{K}{\rho c_p}
-    \end{equation}
-    ```
-    and if we assume that the supersaturation S can be approximated by specific humidities (this is only exactly true for mass mixing ratios):
-    ```math
-    \begin{equation}
-        S_l = \frac{q_{vap}}{q_{sl}}
-    \end{equation}
-    ```
-    then we can write:
-    ```math
-    \begin{equation}
-      q_{vap} - q_{sl} = q_{sl}(S_l - 1)
-    \end{equation}
-    ```
-    and ```Gamma``` is equivalent to the ```G``` function used in our parcel and parameterizations.
+Note that these forms of condensation/sublimation and deposition/sublimation are equivalent to those described in the adiabatic parcel model with some rearrangements and assumptions. It is just necessary to use the definitions of ```\tau```, ```q_{sat}```, and the thermal diffusivity ```D_v```:
+
+```math
+\begin{equation}
+  \tau = 4 \pi N_{tot} \bar{r}, \;\;\;\;\;\;\;\;
+  q_{sat} = \frac{e_{sat}}{\rho R_v T}, \;\;\;\;\;\;\;\;
+  D_v = \frac{K}{\rho c_p}
+\end{equation}
+```
+and if we assume that the supersaturation S can be approximated by specific humidities (this is only exactly true for mass mixing ratios):
+```math
+\begin{equation}
+    S = \frac{q_{vap}}{q_{sat}}
+\end{equation}
+```
+then we can write:
+```math
+\begin{equation}
+  q_{vap} - q_{sat} = q_{sat}(S - 1)
+\end{equation}
+```
+and ```Gamma``` is equivalent to the ```G``` function used in our parcel and parameterizations.
+
+
+### Time integrated version
+
+Finally, we can also use the time integrated formulation of condensation/evaporation and deposition/sublimation described by [MorrisonMilbrandt2015](@cite). We assume that the condensation/evaporation and deposition/sublimation over the course of the timestep can be calculated using the average difference between the $q_v$ value and those of saturation. Ie, they write:
+
+```math
+\begin{equation}
+   \left( \frac{d q_l}{dt} \right) _{cond} = \frac{\bar{\delta_l}}{\tau_l \Gamma_l}
+\end{equation}
+```
+where $\delta$ is defined as
+```math
+\begin{equation}
+    \delta = q_v - q_{sl}
+\end{equation}
+```
+
+The same can be done for deposition/sublimation, with an addition of the correction for the difference between the specific humidities of supersaturation for liquid and for ice:
+
+```math
+\begin{equation}
+   \left( \frac{d q_i}{dt} \right) _{cond} = \frac{\bar{\delta_l}}{\tau_i \Gamma_i} + \frac{q_{sl} - q_{si}}{\tau_i \Gamma_i}
+\end{equation}
+```
+
+In order to calculate $\bar{\delta}$, we consider the derivative of delta:
+
+```math
+\begin{equation}
+    \frac{d \delta}{dt} = \frac{dq_v}{dt} - \frac{dq_{sl}}{dt}
+\end{equation}
+```
+
+We can write the derivation for the saturated specific humidity:
+
+```math
+\begin{equation}
+    \frac{d q_{s}}{dt} = \frac{dq_{s}}{dT} \frac{dT}{dt} + \frac{dq_{s}}{dp} \frac{dp}{dt} = \frac{dq_{s}}{dT} \frac{dT}{dt} + \frac{q_{s} \rho_a g w}{p - e}
+\end{equation}
+```
+where
+- $\rho_a$ is the air density
+- $p$ is the air pressure
+- $e$ is the saturation vapor pressure
+- This makes sense given the assumptions that $\frac{dp}{dt} = -\rho_a gw$ from hydrostatic balance, and $q_s = \frac{e_s}{p-e_s}$  for water vapor, where $e_s$ is the saturation water vapor pressure. Then $\frac{d q_s}{dp} = - \frac{e_s}{(p-e_s)^2} = - \frac{q_s}{p-e_s}$ .
+
+and they write the change in temperature with time as:
+
+```math
+\begin{equation}
+    \frac{dT}{dt} = - \frac{g w}{c_p} + \frac{L_v}{c_p} \frac{\delta}{\tau_l \Gamma_l} + \frac{L_s}{c_p} \left(\frac{\delta}{\tau_i \Gamma_i} + \frac{q_{sl} - q_{si}}{\tau_i \Gamma_i} \right)
+\end{equation}
+```
+
+Then
+
+```math
+\begin{equation}
+       \frac{d q_{sl}}{dt} = \frac{dq_{sl}}{dT} \left[- \frac{g w}{c_p} + \frac{L_v}{c_p} \frac{\delta}{\tau_l \Gamma_l} + \frac{L_s}{c_p} \left(\frac{\delta}{\tau_i \Gamma_i} + \frac{q_{sl} - q_{si}}{\tau_i \Gamma_i} \right) \right] + \frac{q_{sl} \rho_a g w}{p - e}
+\end{equation}
+```
+
+!!! note
+    we neglect terms due to radiation and mixing (included in the [MorrisonGrabowski2008_supersat](@cite) and [MorrisonMilbrandt2015](@cite) versions) in the parcel case.
+
+Putting these back into the $\delta$ equations and rearranging based on the definitions of $\Gamma_l$ and $\Gamma_i$, we get:
+
+```math
+\begin{equation}
+    \frac{d \delta}{dt} = \frac{- q_{sl} \rho_a g w}{p - e_{sl}} + \frac{gw}{c_p} \frac{dq_{sl}}{dT} - \frac{\delta}{\tau}
+\end{equation}
+```
+
+where $\tau$ is
+
+```math
+\begin{equation}
+    \tau = \left( \tau_{l}^{-1} + \left( 1 + \frac{L_{s}}{c_p} \frac{dq_{sl}}{dT} \right) \frac{\tau_i^{-1}}{\Gamma_i} \right)^{-1}
+\end{equation}
+```
+
+where
+
+- $\tau_l$ is the relaxation timescale for liquid droplets
+- $\tau_i$ is the relaxation timescale for ice droplets
+- $L_s$ is the latent heat of sublimation
+- $c_p$ is the specific heat of air at constant pressure
+- $T$ is temperature
+
+If we assume that changes in $\frac{dq_s}{dT}$ and $\tau$ are small in comparison to their magnitudes, then we can approximate them as constant and both $\delta$ equations are linear differential equations with a solution of the form:
+
+```math
+\begin{equation}
+    \delta (t) = A_c \tau + (\delta_{t=0} - A_c \tau) e^{-t/\tau}
+\end{equation}
+```
+
+Then the $A_c$ values can be calculated based on the previous derivatives.
+
+Condensation Ac:
+```math
+\begin{equation}
+        A_c = - \frac{q_{sl} \rho g w}{p - e_s} + \frac{dq_{sl}}{dT} \frac{wg}{c_p} - \frac{(q_{sl} - q_{si})}{\tau_i \Gamma_i} \left( 1 + \frac{L_s}{c_p} \frac{dq_{sl}}{dT} \right)
+\end{equation}
+```
+
+Finally, to get the values for condensation/evaporation and deposition/sublimation over the timestep, we calculate the time average of $\delta$ by dividing by $\Delta t$ (the timestep) and then by the relaxation time as described in the above equations.
+
+```math
+\begin{equation}
+    \left( \frac{d q_l}{dt} \right)_{cond} = \frac{A_c \tau}{\tau_i \Gamma_l} + (\delta_{t=0} - A_c \tau) \frac{\tau}{\Delta t \tau_i \Gamma_l} (1-e^{-\Delta t / \tau} )
+\end{equation}
+```
+
+```math
+\begin{equation}
+    \left( \frac{d q_{si}}{dt} \right)_{dep} = A_c \frac{\tau}{\tau_i \Gamma_i} + (\delta_{t=0} - A_c \tau) \frac{\tau}{\Delta t \tau_i \Gamma_i} (1-e^{-\Delta t / \tau} ) + \frac{(q_{sl} - q_{si})}{\tau_i \Gamma_i}
+\end{equation}
+```