A new parameterisation for homogeneous ice nucleation driven by highly variable dynamical forcings

The present work aims to extend the parameterisation of homogeneous ice nucleation introduced in Dolaptchiev et al. (2023) by incorporating variable ice mean mass and generalising the approach under different conditions. The proposed method involves introducing an empirically derived correction to the parameterisation based on a large data set of parcel model simulations to account for the mean ice particle mass variations. The method is validated against ensemble simulations using time-resolved double-moment ice microphysics, showing a mean relative error of less than 16 %, with robust performance across a range of conditions. The uncertainty of the proposed parameterisation is evaluated for increasing integration time steps. The method remains computationally efficient and produces sufficiently accurate results, even with larger time steps, making it suitable for integration into numerical weather prediction and climate models. It is shown that the generalised approach not only provides a good representation of individual nucleation events but also effectively captures the statistics across the ensemble data. The prediction of the ice mixing ratio is also assessed against the reference double-moment system results. It is demonstrated that, following the nucleation event and re-initialisation of the ice mixing ratio, the system equilibrates toward the reference solution within a few time steps. This refined parameterisation offers a more accurate prediction of ice number concentration and ice mixing ratio, is not limited to perturbations induced by gravity waves (GWs), and can be supplemented by other relevant dynamical effects, such as large-scale motions or even turbulence.

1 Introduction

Cirrus clouds correspond to up to 30 % of the total cloud cover; however, they are still not well understood, and their overall impact on climate in particular on the radiation forcing remains highly uncertain (e.g. Gasparini et al., 2018; Joos et al., 2014; Boucher et al., 2013). The albedo effect and the greenhouse effect of these clouds are similar in magnitude, meaning that the microphysical details, such as the number and shape of the ice crystals, can significantly influence their net radiative effect (Krämer et al., 2020; Matus and L'Ecuyer, 2017; Wang et al., 2020; Zhang et al., 1999). Some studies also indicate the indirect effects of cirrus clouds on future climate (e.g. Gettelman et al., 2012), which occur mostly due to increases in ice crystal concentrations due to homogeneous nucleation from anthropogenic sulfur emissions.

These microphysical characteristics of cirrus are affected by a complex interaction of the small-scale cloud processes and the dynamics of the surrounding atmosphere. Observations suggest that gravity wave (GW) dynamics significantly impact the properties and life cycle of cirrus clouds (Kärcher and Ström, 2003; Atlas and Bretherton, 2023).

Vertical velocity perturbations, associated with gravity waves, turbulence, or other processes, lead to the homogeneous freezing of aqueous solution droplets (Koop et al., 2000; Baumgartner et al., 2022; Kärcher et al., 2024) at low-temperature conditions in the tropopause region. The effect of gravity waves and other local dynamical processes on the nucleation has been widely investigated utilising different approaches to describe the updrafts. Those interactions are studied, for example, in Dinh et al. (2016) and Kaercher and Podglajen (2019) utilising a stochastic approach, in Joos et al. (2008) using a GW linear theory, and in Zhou et al. (2016) using the large-scale grid-resolved updraft velocity alongside subgrid-scale turbulence-induced updrafts.

Based on the results from Baumgartner and Spichtinger (2019) studying the homogeneous nucleation due to constant updraft velocities, Dolaptchiev et al. (2023) has extended the asymptotic approach to incorporate GW dynamics. The present work complements Dolaptchiev et al. (2023) and has been designed for further usage in connection with GW parameterisations for a more realistic representation of GW–cirrus interactions in coarse-grid models. Recently, transient and lateral propagation effects have been introduced into GW parameterisations (Kim et al., 2021; Achatz et al., 2023; Voelker et al., 2024), offering further potential for a more realistic simulation of these interactions.

The mean mass variability in the ice particles has a significant effect on the cloud ice number concentration prediction, which is also observed in numerical parcel model simulations in Dolaptchiev et al. (2023). The influence of a cloud ice particle size uncertainty in a climate model was assessed in Wang et al. (2020) and illustrates the sensitivity of the surface precipitation to the change in particle effective radius due to the larger energy imbalance and due to spatial variability in short- and longwave radiation. The effects of variable mass become more prominent in conditions of a higher number of pre-existing ice particles. Finally, it was previously shown (Gierens et al., 2003) that the deposition (or accommodation) coefficient, controlling the growth of ice crystals in the kinetic regime, largely affects the nucleation process, while at the same time being dependent on the radius or mass of ice crystals.

The current work focuses on generalising the approach introduced in Dolaptchiev et al. (2023) to include mean mass variability for a wide range of initial conditions, e.g. dynamical forcing and environmental conditions. The asymptotic approach is supplemented by an empirical correction that can be used for various forcing representations. The proposed method is not restricted to the GW-induced dynamics and can be used for any kind of air motion, e.g. as large-scale updrafts or even turbulent fluctuations. This work is also relevant for a further implementation into a global numerical weather prediction (NWP) model, where it can be coupled to the dynamical parameterisations (e.g. GW and turbulence). However, due to the generalised formulation of the approach and robustness for larger time steps, it can also be applied to climate models with large time steps.

The article is organised as follows: Sect. 2.1 provides a brief overview of the approach used and the motivation for its generalisation. Section 2.2 and 2.3 describe the chosen approach for constructing the parameterisation and technical details of the data set used. In Sect. 3, the proposed parameterisation with correction is verified against reference calculations, and the range of applicability is tested based on ensemble calculations. Concluding remarks can be found in Sect. 4.

2 Theory and methodology

2.1 Double-moment approach to ice physics representation coupled to GW

In the present work, cirrus clouds are modelled using a double-moment bulk microphysics scheme (based on Spichtinger et al., 2023) that assumes a unimodal ice mass distribution function, where the zeroth and first moments of the distribution function are ice number concentration n_i and ice mixing ratio q_i, respectively. Ice crystals are assumed to have a spherical shape, which simplifies the description of the properties. The assumption mentioned is reasonable for very small, newly nucleated ice particles formed through homogeneous nucleation. However, as ice particles grow through deposition, their shape becomes more complex, introducing some uncertainty in the calculations. For consistency, we maintain the same shape assumption in deposition calculations, acknowledging that this introduces a certain error. Although a more sophisticated shape approximation would be beneficial for studying ice growth dynamics, this does not affect the developed correction and goes beyond the topic of the current study; former nucleation parameterisations also assume spherical ice particles at their early stage of life (see e.g. Kärcher and Lohmann, 2002; Ren and Mackenzie, 2005; Barahona and Nenes, 2008). The ice physics scheme used here is described in Spichtinger and Gierens (2009) and Spreitzer et al. (2017) with the addition of sedimentation sinks. The influence of Lagrangian pressure and temperature variations, e.g. due to GWs, on homogeneous nucleation and on mean mass variability is incorporated in Dolaptchiev et al. (2023). The resulting system of equations for n_i,q_i, and the saturation ratio over ice S reads as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{D}{n}_{\text{i}}}{\mathrm{D}t}}=J\exp \left(B\right(S-S_{\text{c}}\left)\right)+{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\partial \left(\mathit{\rho }n{v}_{n_{\text{i}}}\right)}{\partial z}},\text{(2)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{D}{q}_{\text{i}}}{\mathrm{D}t}}& ={\displaystyle \frac{p_{\text{si}}}{p}}D\left({\displaystyle \frac{q_{\text{i}}}{n_{\text{i}}}}{)}^{\mathrm{1}/\mathrm{3}}\right(S-\mathrm{1})T{n}_{\text{i}}\\ & +\widehat{m}J\exp \left(B\right(S-S_{\text{c}}\left)\right)+{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\partial \left(\mathit{\rho }q{v}_{q_{\text{i}}}\right)}{\partial z}},\end{array}\text{(3)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{D}S}{\mathrm{D}t}}& =-D\left({\displaystyle \frac{q_{\text{i}}}{n_{\text{i}}}}{)}^{\mathrm{1}/\mathrm{3}}\right(S-\mathrm{1})T{n}_{\text{i}}-{\displaystyle \frac{p}{p_{\text{si}}}}J\exp (B(S-S_{\text{c}}))\\ & -S{\displaystyle \frac{\mathrm{1}}{\mathit{\pi }}}{\displaystyle \frac{\mathrm{D}\mathit{\pi }}{\mathrm{D}t}}({\displaystyle \frac{L_{\text{i}}}{T{R}_{\text{v}}}}-{\displaystyle \frac{c_{\text{p}}}{R}}).\end{array}\end{array}$$

Here the saturation ratio over ice is described as $S=\frac{p_{\text{v}}}{p_{\text{si}}}$; p_v is the water vapour pressure; p_si is the saturation pressure over ice; p is the pressure; $\widehat{m}$ is the mean mass of newly nucleated ice particles; J is a homogeneous nucleation rate, also including the number concentration of solution droplets; B is a nucleation parameter; S_c is a critical saturation ratio; D is a growth coefficient, including the diffusion coefficient and the Howell factor; v${}_{{\text{n}}_{\text{i}}}$, and v${}_{{\text{q}}_{\text{i}}}$ are sedimentation velocities; T is the temperature; ρ is the density; π is the Exner pressure; L_i denotes the latent heat of sublimation; c_p is the specific heat capacity of air at constant pressure; R_v is the ideal gas constant for water vapour; and R is the ideal gas constant for dry air. The values of the parameters used are presented in Table A1. In general formulation, the critical value of the saturation ratio over ice is temperature-dependent; in the current work, the value is kept constant, while other studies (such as Spichtinger et al., 2023) show an insignificant influence of S_c on nucleation processes. Equations (1)–(3) are written in terms of material derivative $\frac{\mathrm{D}}{\mathrm{D}t}$; i.e. all quantities are described following the motion of air parcels. The right-hand parts of the system of equations describe the ice crystal growth due to the deposition of water vapour, the formation of new ice crystals through homogeneous nucleation, and the sedimentation of the ice crystals due to gravity. The latter is not further included in the consideration because of the much larger sedimentation timescale (see Dolaptchiev et al., 2023).

In Baumgartner and Spichtinger (2019), detailed air parcel simulations are conducted using a box model and a bulk microphysics scheme (based on Spichtinger and Gierens, 2009), shown for various environmental conditions, assuming a constant mean mass is generally valid over the duration of the nucleation. However, in case of pre-existing ice particles just before nucleation, when the saturation ratio exceeds 1, the ice crystals' mean mass increases, leading to a higher deposition rate. This increased deposition affects the saturation ratio over ice (as a sink), which in turn influences the number of nucleated ice crystals, reducing the actual nucleation rate. Given that, it is crucial to include these effects in the current model, where the mean mass of ice crystals can be determined from the relation $q_{\text{i}}/n_{\text{i}}$. This is achieved by incorporating the evolution equation for q_i and substituting $m\left(t\right)=q_{\text{i}}/n_{\text{i}}$ in the diffusional growth term.

Timescale separation is employed in Dolaptchiev et al. (2023) to assess the different processes and link the ice physics equations with the effects of GWs. The details of the coupling between the processes' timescales and assumptions of the GW theory are presented in Appendix A.

Further study is done considering the ice physics of a single air parcel influenced by GW dynamics, adopting a Lagrangian framework. In this case, material derivatives in the system of equations are replaced by time derivatives $\frac{\mathrm{D}}{\mathrm{D}t}=\frac{\mathrm{d}}{\mathrm{d}t}$.

As shown in Dolaptchiev et al. (2023), considering only dominant processes, the system of Eqs. (1)–(3) can be reduced to the following form:

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\mathrm{d}{n}_{\text{i}}}{\mathrm{d}t}}=J\exp \left(B\right(S-S_{\text{c}}\left)\right),\text{(5)}& {\displaystyle \frac{\mathrm{d}{q}_{\text{i}}}{\mathrm{d}t}}={\displaystyle \frac{p_{\text{si}}}{p}}D\left({\displaystyle \frac{q_{\text{i}}}{n_{\text{i}}}}{)}^{\mathrm{1}/\mathrm{3}}\right(S-\mathrm{1})T{n}_{\text{i}},\text{(6)}& {\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=-D\left({\displaystyle \frac{q_{\text{i}}}{n_{\text{i}}}}{)}^{\mathrm{1}/\mathrm{3}}\right(S-\mathrm{1})T{n}_{\text{i}}+SF(t).\end{array}$$

A forcing term F(t) incorporating constant background updraft and GWs dynamics is written in the following form:

$$F\left(t\right)={\displaystyle \frac{g{L}_{\text{i}}}{c_{\text{p}}{R}_{\text{v}}{T}^{\mathrm{2}}}}(w_{\mathrm{00}}+\sum _j{\widehat{w}}_j\cos ({\mathit{\omega }}_{j}t+{\mathit{\varphi }}_j\left)\right),$$

where g is the acceleration due to gravity; w₀₀ is a constant background updraft; and ${\widehat{w}}_j$, ω_j, and ϕ_j are the vertical wind amplitude, frequency, and phase of GW number j. It is worth mentioning that the application of the full system of Eqs. (1)–(3) is not a standard approach, whereas the reduced system (Eqs. 4–6) is close to the reduced model in Baumgartner and Spichtinger (2019), and variations of it are commonly used in NWP and climate applications but including the ice mass concentration.

Referring to the previous results from Dolaptchiev et al. (2023), in the case of constant mass, the asymptotic technique can be applied to obtain the prototype parameterisation of the homogeneous nucleation process. By employing asymptotic expansions, the authors systematically decompose the governing equations into different scale regimes corresponding to the key microphysical processes: nucleation and deposition. In the nucleation regime, where supersaturation rapidly increases due to gravity-wave-induced cooling, they derive an analytical solution for the evolution of ice crystal number concentration. In the growth (or deposition) regime, where ice crystals grow primarily by vapour deposition, the governing equations are solved separately under the assumption of gradual saturation adjustment. The solutions from both regimes are then matched using asymptotic matching techniques. This approach results in a reduced analytical expression for the number of nucleated ice crystals per nucleation event, explicitly linking it to wave parameters (such as frequency and amplitude) and background atmospheric conditions (such as temperature and humidity). The nucleation event is roughly determined where the saturation ratio reaches the critical value S_c, and the predicted value for the ice number concentration is found as follows:

$$\begin{array}{}\text{(7)}& n_{\text{post}}=\left\{\begin{array}{ll}\mathrm{2}\frac{S_{\text{c}}F\left(t_{\mathrm{0}}\right)}{D{m}_{\text{c}}^{\mathrm{1}/\mathrm{3}}T(S_{\text{c}}-\mathrm{1})}-n_{\text{pre}},& \text{if }{n}_{\text{pre}}<\frac{S_{\text{c}}F\left(t_{\mathrm{0}}\right)}{D{m}_{\text{c}}^{\mathrm{1}/\mathrm{3}}T(S_{\text{c}}-\mathrm{1})};\\ n_{\text{pre}},& \text{otherwise}.\end{array}\right.\end{array}$$

Here $m_{\text{c}}={\mathrm{10}}^{-\mathrm{12}}$ kg is the reference mass, n_pre and n_post are the ice number concentrations before and after the nucleation event, and t₀ is the time of a nucleation event when S=S_c. Note that nucleation already occurs for values of S close to S_c; see asymptotic analysis in Spichtinger et al. (2023).

In Dolaptchiev et al. (2023), the asymptotic approach was extended to account for the variable mass in the depositional growth coefficient in Eqs. (5) and (6) under the regime of slowly varying the ice mixing ratio q_i during the nucleation event. In such a case, the different processes can be decoupled, and only n_i undergoes a quick change during the nucleation regime.

Numerical studies of the reduced system (4)–(6) based on a simplified parcel model showed several possible regimes relevant to the nucleation process (see Fig. 1). One regime at high values of the pre-existing ice crystal number concentration (shown in Fig. 1 in black) follows the same tendency as a constant-mass case. The evolution of parameters can be separated into three phases (see Baumgartner and Spichtinger, 2019), where the nucleation part describes only the phase of rapid change in ice number concentration. In contrast, the saturation ratio and ice mixing ratio change slowly during this short time and can therefore be assumed not to vary during this phase. However, conditions of low to zero pre-existing ice (shown in red) lead to a different behaviour due to the growth coefficient change. Assumptions of the time separation for the nucleation regime do not hold, since the ice mixing ratio is changing on the same timescale as the ice number concentration. The extension to be proposed here aims to cover the whole range of initial conditions, allowing the use of the parameterisation accounting for mass variability independently from the regime.

Figure 1 Evolution of the ice number concentration n_i, ice mixing ratio q_i, saturation ratio over ice S, and mean mass of ice particles $m=q_{\text{i}}/n_{\text{i}}$ calculated based on system (4)–(6). The two cases of initial conditions are low pre-existing ice n_init=10² kg⁻¹ and higher values of pre-existing ice n_init=10⁶ kg⁻¹. The blue line corresponds to the critical saturation ratio, and dashed lines mark the time of the nucleation event, where S=S_c is in red for the case n_init=10² kg⁻¹ and in black for the case n_init=10⁶ kg⁻¹.

2.2 Correction of the growth coefficient for variable mean mass

In the present study, we utilise the previously obtained result for the constant-mass case (7). This approach would allow a clear physical relation between n_post, the initial number concentration n_init, and the forcing term F(t₀) at the point t₀ of the nucleation event. It would also account for different background conditions for temperature by design. Since mass variability changes the growth coefficient directly, as one can see from Eq. (6), the straightforward solution is to correct the growth coefficient in the constant-mass parameterisation to a factor representing the mass at the point of the nucleation event. The corrected parameterisation has the following form:

$$\begin{array}{}\text{(8)}& n_{\text{post}}=\left\{\begin{array}{ll}\mathrm{2}\frac{S_{\text{c}}F\left(t_{\mathrm{0}}\right)}{D\left(m_{\text{eff}}{)}^{\mathrm{1}/\mathrm{3}}T\right(S_{\text{c}}-\mathrm{1})}-n_{\text{pre}},& \text{if }{n}_{\text{pre}}<\frac{S_{\text{c}}F\left(t_{\mathrm{0}}\right)}{D\left(m_{\text{eff}}{)}^{\mathrm{1}/\mathrm{3}}T\right(S_{\text{c}}-\mathrm{1})};\\ n_{\text{pre}},& \text{otherwise},\end{array}\right.\end{array}$$

where $(m_{\text{eff}}{)}^{\mathrm{1}/\mathrm{3}}$ is a growth coefficient correction factor representing the effective mean mass of the ice crystals at the start of the nucleation event.

In order to find the desired growth coefficient, one needs to define the parameters that can be used for the correction. Figure 2 shows the evolution of the parameters calculated from the reduced system (4)–(6) and known pre-nucleation regime parameters, which can be used for deriving the correction and parameterisation. Here a pre-nucleation regime is defined as the solution of the Eqs. (4)–(6) system in the regime, where initial ice number concentration stays constant, which leads to the solution of Eqs. (5) and (6) with $n=n_{\text{init}}=\mathrm{const}$ decoupled from Eq. (4). Sharp changes in the ice mixing ratio and mass in the case of a reduced system make it impossible to approximate the necessary quantities q_i(t₀) or $\frac{q_{\text{i}}\left(t_{\mathrm{0}}\right)}{n_{\text{i}}\left(t_{\mathrm{0}}\right)}$ from pre-nucleation regime values alone. Therefore, the introduced effective mass is approximated based on a large data set. Variables that can be used in the approximation include the forcing term at the point of the nucleation event F(t₀), the initial number concentration n_init, and the mass m^* predicted from the pre-nucleation regime, assuming no change in ice number concentration before the nucleation event.

Figure 2 Evolution of $n_{\text{i}},q_{\text{i}},S,m$ in the reduced system (4)–(6) case and in the pre-nucleation regime, solving Eqs. (5) and (6) with assumption of $n_{\text{i}}=n_{\text{init}}=\mathrm{const}$. The blue line corresponds to the critical saturation ratio, and the dashed lines mark the time of the nucleation event, where S=S_c.

The exact effective mass m_eff is determined from the results of the full system, such that the constant-mass approach leads to the observed n_post. Assuming a known evolution of the parameters during the nucleation event, Eq. (8), under the conditions allowing further nucleation, can be solved for m_eff:

$$\begin{array}{}\text{(9)}& m_{\text{eff}}=({\displaystyle \frac{\mathrm{2}{S}_{\text{c}}F\left(t_{\mathrm{0}}\right)}{(n_{\text{post}}+n_{\text{pre}})DT(S_{\text{c}}-\mathrm{1})}}{)}^{\mathrm{3}},\end{array}$$

where t₀, n_pre, and n_post are known from the resolved full system (1)–(3). A similar approach was used in Baumgartner and Spichtinger (2019), where the increase in the ice mass during the nucleation was found to be linear with a slope proportional to some constant mean mass.

Next, we discuss the construction of the data set (Kosareva, 2025b) from ensemble simulations based on Eqs. (1)–(3), disregarding the sedimentation process. This simplification can be justified by the small sizes of ice particles during and shortly after the nucleation event. The variability during the nucleation event depends on the initial conditions, such as ice number concentration, ice mixing ratio, or mean mass of ice particles, and on the GW forcing. Since we want to represent the nucleation process itself, the initial saturation ratio is set to a value of $S(t=\mathrm{0})=\mathrm{1.4}$, close to the critical value. Other initial conditions are varied randomly using a uniform distribution in logarithmic space, i.e. log (n_init) in physically meaningful ranges motivated by observations (Krämer et al., 2016): we choose 10⁻⁴ kg⁻¹ $<n_{\text{init}}(t=\mathrm{0})<{\mathrm{10}}^{\mathrm{7}}$ kg⁻¹ and 10⁻¹⁶ kg $<m(t=\mathrm{0})<{\mathrm{10}}^{-\mathrm{12}}$ kg. The lower bound for the initial number concentration is chosen to represent the conditions of low to no pre-existing ice, which can occur in applications to the global NWP models.

Construction of the forcing is done based on the output from the global ICOsahedral Non-hydrostatic (ICON) model (Zängl et al., 2015) coupled with the GW parameterisation Multi-Scale Gravity Wave Model (MS-GWaM) (Bölöni et al., 2021; Kim et al., 2021; Achatz et al., 2023; Kim et al., 2024; Voelker et al., 2024). This parameterisation is based on Wentzel–Kramers–Brillouin (WKB) theory and is implemented using Lagrangian ray volumes, which are regarded as carriers of the GW fields' wave action density (Bölöni et al., 2021). The current version of MS-GWaM allows us to account for convectively generated GWs and GWs generated from processes other than convection or flow over mountains. The influence of the GWs generated by the mountains is incorporated by using the original scheme (Lott and Miller, 1997) included in the ICON model. For constructing the forcing term, we utilise ICON version 2.6.5-nwp1, with 120 vertical levels and at horizontal resolution R2B5 (∼ 80 km). The initial conditions are taken for 1 May 2010, and the model is integrated forward for 3 weeks, accounting for 2 weeks of spin-up. Information on subgrid-scale gravity-wave-related vertical wind perturbations is retrieved from MS-GWaM using the ray volume with the maximum vertical wind amplitude per cell and the corresponding frequency. In order to cover the variety of possible forcing-term configurations, the following cases are included in the data set construction:

1. Large-scale vertical velocity from ICON run w₀₀. In such a case, the forcing term reads $F\left(t\right)=\frac{g{L}_{\text{i}}}{c_{\text{p}}{R}_{\text{v}}{T}^{\mathrm{2}}}{w}_{\mathrm{00}}$.

2. A single GW from MS-GWaM with maximum wind amplitude $\widehat{w}$ within the cell, the corresponding frequency ω, and initial phase ϕ which is picked randomly: $F\left(t\right)=\frac{g{L}_{\text{i}}\widehat{w}}{c_{\text{p}}{R}_{\text{v}}{T}^{\mathrm{2}}}\cos (\mathit{\omega }t+\mathit{\varphi })$.

3. A superposition of GWs, where the main GW is taken from the MS-GWaM and other GWs have randomly varying amplitudes A_j, frequencies ω_j, and phases ϕ_j independent from the main GW and rescaled in a way that the total momentum flux added is 1mPa, similar to approach used in Dolaptchiev et al. (2023): $F\left(t\right)=\frac{g{L}_{\text{i}}\widehat{w}}{c_{\text{p}}{R}_{\text{v}}{T}^{\mathrm{2}}}\cos (\mathit{\omega }t+\mathit{\varphi })+{\sum }_j\frac{g{L}_{\text{i}}{A}_j}{c_{\text{p}}{R}_{\text{v}}{T}^{\mathrm{2}}}\cos ({\mathit{\omega }}_{j}t+{\mathit{\varphi }}_j)$.

4. As in case 2 in the presence of (large-scale) background updraft w₀₀,
   $F\left(t\right)=\frac{g{L}_{\text{i}}}{c_{\text{p}}{R}_{\text{v}}{T}^{\mathrm{2}}}(w_{\mathrm{00}}+\widehat{w}\cos (\mathit{\omega }t+\mathit{\varphi }\left)\right)$.

5. As in case 3 in the presence of (large-scale) background updraft w₀₀,
   $F\left(t\right)=\frac{g{L}_{\text{i}}}{c_{\text{p}}{R}_{\text{v}}{T}^{\mathrm{2}}}(w_{\mathrm{00}}+\widehat{w}\cos (\mathit{\omega }t+\mathit{\varphi })+{\sum }_j\frac{g{L}_{\text{i}}{A}_j}{c_{\text{p}}{R}_{\text{v}}{T}^{\mathrm{2}}}\cos ({\mathit{\omega }}_{j}t+{\mathit{\varphi }}_j\left)\right)$.

It is worth mentioning that the GW diagnosed from the MS-GWaM model dominates in cases 3 and 5 due to the larger amplitudes of the vertical velocity perturbations.

Distributions of the large-scale vertical winds and vertical winds from the GWs used for constructing the forcing are presented in Fig. 3a, b. The probabilities for the corresponding frequency ω are also presented, along with those for the calculated GW momentum flux ρuw for the considered GW (see Fig. 3c, d). The data are taken from the instantaneous output of the global ICON simulation at altitudes from 8 to 14 km, most relevant for cirrus formation. While this range does not explicitly cover the full tropical tropopause layer (TTL), including altitudes up to 18 km, we note that the selected range captures the majority of ice formation conditions observed in the upper troposphere/lower stratosphere (UTLS). Moreover, the vertical velocity distributions in this range are comparable to those at higher altitudes, and the proposed parameterisation is additionally tested under conditions representative of stronger forcing (see Appendix C). These tests suggest that the approach remains robust and physically consistent even beyond the original range of sampled data.

It is seen in Fig. 3 that the large-scale updraft varies in range from −0.2 up to 0.4 m s⁻¹, where, in the majority of the cells, vertical velocity is close to zero. Therefore, we include cases (2) and (3) without the large-scale updraft when creating the data set. Vertical wind perturbations $\widehat{w}$ obtained from the GW parameterisation show stronger vertical velocity, with a rare maximum above 2 m s⁻¹. Corresponding frequencies ω are sampled in a range between a Coriolis parameter of $f={\mathrm{10}}^{-\mathrm{4}}$ s⁻¹ and a buoyancy frequency of $N=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{2}}$ s⁻¹ in the tropopause region, and the form of the ω spectrum resembles the Desaubies form (VanZandt, 1982). It is worth noting that the resulting momentum flux based on the highest-amplitude ray volume is aligned with the observational statistics obtained from the balloon measurements (Hertzog et al., 2012) and previously calculated momentum fluxes from the MS-GWaM parameterisation (Kim et al., 2021). Therefore, the total momentum flux considered lies in the physical range, and the forcing term applied in the parcel simulations is realistic to some degree. Further support for the realism of the applied conditions comes from Fig. 3b, where the probability density function (PDF) of gravity-wave-induced vertical velocity perturbations is $\widehat{w}$. The spread and frequency of occurrence show a good agreement with available observations from balloon measurements (Podglajen et al., 2016; Köhler et al., 2023).

Figure 3 Frequencies of occurrences of (a) the large-scale vertical velocity w, (b) the vertical velocity amplitude $\widehat{w}$, (c) the frequency of the GWs ω defined from the MS-GWaM parameterisation, and (d) the momentum flux corresponding to the detected GW.

The various vertical velocities used in the construction of the forcing are sampled from the wide range of conditions, including the large updrafts up to 10 m s⁻¹. Because of the general model formulation, one can add the vertical velocity contributions to the forcing term based on other information provided by the model, for example, by the turbulence scheme. However, such an extension of the forcings is not tested in the current work.

The exact value of the effective mass m_eff is found using Eq. (9) based on the ensemble calculations for $\sim \mathrm{2}\times {\mathrm{10}}^{\mathrm{5}}$ nucleation events; the results are shown in Fig. 4. The data set is created using random initial conditions for n_i, q_i, and various GW-forcing configurations as discussed above. The internal variability in the data is quite high, especially depending on the forcing term at the point of the nucleation event. However, by fixing certain parameters (blue and black markers), one can see a clear dependence of m_eff on the chosen parameters. The variability with the forcing term is much larger, whereas, under the conditions of smaller initial number concentrations, m_eff stays nearly constant. The dependence on the forcing term is decaying with large variability for smaller values of the forcing at the point of the nucleation event. The increase in the mass with an increase in the initial number concentration is connected to the fact that pre-existing ice suppresses the nucleation process; hence, the number of newly nucleated ice crystals is smaller. Because the initial ice mean mass is quite high, and the new particles are 4 orders of magnitude lighter, the recalculated mean mass increases if more pre-existing ice is present. The observed behaviour of the mass in dependence on the forcing term is the opposite: it decreases with increasing F(t₀). A larger forcing term at the point of the nucleation is associated with the larger updrafts leading to a faster change in the saturation ratio over ice. The depletion of the water vapour via deposition in such conditions is ineffective, leading to smaller mass values.

Figure 4 The exact m_eff correction depending on the initial number concentration (a) and on the forcing-term value at the nucleation event point (b). Blue markers correspond to calculations with constant n_init, and black markers correspond to calculations made with approximately constant values of F(t₀).

2.3 Mean mass fitting for parameterisation

In the present section, we consider options for approximating the dependence of m_eff based on the parameters known from the pre-nucleation regime. Since the evolution of the ice mixing ratio q_i is not captured in this regime (due to fixed growth coefficient and n_i), we propose to fit a function depending on the known parameters, such as initial ice number concentration n_init and forcing term at the point of the nucleation event F(t₀), to the ensemble simulation results. These parameters are selected based on the definition of the correction given in Eq. (9), and their relevance is further supported by the clear dependence observed in Fig. 4, where only a slight spread is visible in cases with fixed n_init and F(t₀).

Based on the created data set and the relative behaviour of m_eff in dependence on n_init and F(t₀), an approximation to the data is derived. The non-linear least-squares method is applied to fit an assumed function to the known data from the ensemble simulations. The method requires us to suggest a certain functional dependency on the parameters. Given the data set behaviour, the chosen simplified function is

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{m}_{\text{eff}}\left(F\right(t_{\mathrm{0}}),n_{\text{init}})& =\exp (a_{\mathrm{1}}+a_{\mathrm{2}}F(t_{\mathrm{0}}{)}^{\mathrm{1}/\mathrm{3}}+a_{\mathrm{3}}{n}_{\text{init}}^{\mathrm{1}/\mathrm{3}}\\ & +a_{\mathrm{4}}{n}_{\text{init}}F\left(t_{\mathrm{0}}{)}^{\mathrm{1}/\mathrm{3}}\right).\end{array}\end{array}$$

The mean error for m_eff during the target data set reproduction is around 30 %, which translates into a mean error of ∼ 10 % in radius, for the set of ∼ 2 ×10⁵ nucleation events with different GW-forcing representations.

As an extension, a more complex dependency on n_init and F(t₀) and involving one more parameter m^*, defined from the pre-nucleation regime, is introduced and tested. The results show an improvement compared to the previously introduced fit (10); see the corresponding figures and discussion in Appendix B. While the use of an advanced fit improves the prediction of the effective mass (m_eff), the ice number concentration after nucleation (n_post), which is the primary target of the parameterisation, is already well captured by the simple fit. The more complex approach yields only a marginal improvement of about 2 % in the mean error. Besides, the advanced fit adds more complexity to the calculation; requires the use of an additional parameter; and, because of the additional terms used, can lead to diverging results in cases out of sample range with regard to the reference data used for the parameter optimisation. Therefore, given the efficiency and robustness of the original simple fit, we recommend its use and present all further analysis based on this version.

3 Results: evaluation of the model

Further results are based on the parcel model code (Kosareva, 2025a) and are shown for an independent data set consisting of 490 nucleation events in order to validate the chosen fit. Initial conditions for n_i and m_i were created randomly, and the forcing term based on ICON simulation is different from the one used for the construction of the fit. In the following, we will first assess the proposed fit. This will be followed by stability and accuracy tests, depending on the time step chosen for the integration of Eqs. (5)–(6) in the parameterised case.

We acknowledge that developing a parameterisation based on a limited data set is inherently challenging and introduces certain limitations. The considered conditions are representative of the UTLS region and are specific to the processes as captured by the ICON model. As a result, extreme or anomalous cases may be underrepresented. Nevertheless, we argue that the physically grounded formulation of the proposed parameterisation, along with its robustness to variations in forcing, allows it to provide reasonable predictions even outside the original fitting range. To support this claim, we validated the parameterisation on an independent data set and included results for selected test cases with unusually strong vertical velocities in Appendix C.

3.1 Extended parameterisation

A comparison of the application of the chosen fit (10) and the exact effective mass calculated from the nucleation event simulations is shown in Fig. 5. A direct comparison of predictions with the exact values for m_eff shows a good agreement for the major part of the data points. Nonetheless, the fit is not able to capture the higher possible m_eff values at higher values of n_init. Figure 6 shows the exact post-nucleation ice number concentration and the one obtained with the use of the extended parameterisation (8). The underestimation of the m_eff values at higher values of pre-existing ice number concentrations leads to an overestimation of the ice number concentration as shown in Fig. 6.

Figure 5 Exact m_eff correction depending on the initial number concentration and m_eff prediction from the constructed fit (10) for independent ensemble calculations.

Referring to Fig. 5, the most frequently occurring cases of m_eff in the range between 10⁻¹³ and 10⁻¹² kg and the corresponding final n_post in Fig. 6 in the range between 10⁵ and 10⁶ are successfully represented by the chosen fit (10) and the corrected parameterisation (8). However, the outlying extreme cases are poorly captured due to the simplicity of the fit.

Figure 6 Post-nucleation ice number concentrations defined from the full system and prediction from the constructed parameterisation with correction (10).

The overall behaviour of the fit stays similar to the exact data; therefore, with decreasing n_init, the prediction for the m_eff stays nearly independent from n_init (see Fig. 5a). The cases with an increased forcing term correspond to the smaller values of m_eff, increasing the final prediction for n_post in the exact data, and the same behaviour is captured by the fit.

To assess how well the chosen approach represents individual cases, the relative errors of m_eff, effective ice particle radius R_eff, and n_post with respect to their exact values are calculated. The effective radius is calculated based on the effective mass as follows: $R_{\text{eff}}={\left(\frac{\mathrm{3}{m}_{\text{eff}}}{\mathrm{4}\mathit{\pi }{\mathit{\rho }}_{\text{i}}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}}$. The probability density function (PDF) of the relative errors is presented in Fig. 7. Here the relative errors are calculated as $\mathrm{\Delta }M=(M_{\text{predicted}}-M_{\text{exact}})/M_{\text{exact}}\cdot \mathrm{100}$ %, with $M=m_{\text{eff}},R_{\text{eff}},n_{\text{post}}$. The mean errors for m_eff and n_post predictions are 24 % and 16 %, respectively. Despite the substantial deviation in effective mass, the prediction for the effective radius shows a good agreement, and the mean error is only about 10 %. More importantly, 90 % of the error for the prediction of the ice number concentration lies below 37 %. The majority of the individual cases are represented sufficiently well, whereas outliers and poorly captured events with a deviation in n_post above 50 % correspond to less than 10 % of the distribution.

Figure 7 Probability density functions of the relative errors in the prediction of m_eff, R_eff, and n_post with respect to the exact solutions. Vertical lines correspond to the 90th, 50th, and 10th percentiles of the distribution.

For a better understanding of the prediction deviations, we also consider the PDF of the exact parameters and the approximated ones in order to assess the ability of our fit to capture the statistical properties of the data set. Previous ensemble simulations using constant-mass parameterisation by Dolaptchiev et al. (2023) showed that, while individual events can vary significantly or even be missed, the statistical properties for the ensemble still align well.

Comparisons of the PDFs for the effective mass m_eff and the predicted ice number concentration are shown in Fig. 8.

Figure 8 Probability density functions for m_eff and n_post exact ensemble solutions and predictions. Vertical lines correspond to the mean of the distribution.

The PDF of the predicted mass from Fig. 8 at the point of the nucleation exhibits a slight underestimation of the mean value, due to the inability to correctly predict larger m_eff (also shown in Fig. 5). For the same reason, the tail of the distribution is shorter, and the maximum value has a mismatch of approximately 1 order of magnitude. The lower values of m_eff are represented better, and the minimum values are in the same range. A comparison of the PDFs of the ice number concentration shows a good correspondence of major statistical indicators, such as mean, maximum, and minimum values. However, the extended parameterisation gives occasional overestimation of n_post values for larger values of the forcing term and underestimation of the n_post in cases of small forcing (see the right panel in Fig. 6).

Apart from the values of n_i, the prediction for the ice mixing ratio q_i is made based on the estimated mass m_eff and the ice number concentration n_post, calculated using the corrected parameterisation, via q_{i, post}=m_eff n_post. Such an approach contains some physical inconsistencies because the approximated ice mean mass m_eff differs from the mean mass after the nucleation event, which is supposed to be used. However, the estimation gives a suitable approximation for the further simulation of q_i. Furthermore, the evolution of the parameters for two isolated cases is presented in Figs. 13 and 14 to show the q_i prediction relevance.

3.2 Influence of the time step

The results shown in Sect. 3.1 are conducted using the same time step for the solution of the full system (1)–(3) and coupled system (5)–(6) employed in the parameterisation. Such a decision is made to keep the results comparable without introducing further uncertainties. However, the major goal of the parameterisation is to increase the numerical performance, especially by using larger time steps as in coarse-resolution models, e.g. NWP and climate models. Therefore, the following section is devoted to the depiction of the uncertainties and evaluation of the results under different time steps.

The investigation of the stability of numerical integration and applicability of the defined approach is done based on a data set with six GWs and a constant background updraft w₀₀. One of the sources of uncertainty is the accuracy of the detection of the forcing term at the time of the nucleation where S=S_c. The nucleation event can be misplaced, changing the defined time of the nucleation event t₀, depending on the time step used for the calculation of the saturation ratio evolution (6). The calculation of the forcing term is sensitive to the changes in t₀, which leads to further errors in the estimation of the F(t₀) and other parameters.

Furthermore, the relative errors of the predictions from the exact values ΔM are calculated as mentioned before in Sect. 3.1 for the considered parameters $M=F\left(t_{\mathrm{0}}\right),m_{\text{eff}},n_{\text{post}}$. Firstly, in Fig. 9, we compare relative errors in the forcing F(t₀) when solving slow-scale processes for S,q_i based on system (5) and (6) with five different time step choices. In Figs. 9–11, horizontal bars represent percentiles of the distribution to illustrate data variability. Note that, in several cases, the 10th percentile lies near zero and may not be easily distinguishable in the plots due to scale limitations. The smallest chosen time step of 1 s gives the exact solution, whereas longer time steps lead to deviations from the exact prediction. Such deviations may influence the estimated m_eff via errors introduced in t₀ and F(t₀) determination; however, this would have a crucial impact on the n_post prediction by directly changing F(t₀) in the formula (8).

Figure 9 Probability density function of the relative error of the parameterised forcing term F(t₀) compared to the value detected from the exact full-system solution in %. Horizontal bars correspond to the 90th, 50th, and 10th percentiles of the distribution.

Maximum differences above 50 % occur rarely and are associated with outliers, for instance, some cases where the nucleation event is missed with the usage of the larger time step. One such case is illustrated further in this section in Fig. 13. The mean relative error in the forcing term corresponding to the case of 30 s time steps is 2.6 %, and the 90th percentile is 3.5 %. The mean and the 90th percentile of the PDF are doubled when the time step is doubled and rise up to 15 % and 31 % for the largest considered time step of 2 min. Considering the exact differences instead of relative errors, it is found that larger time steps mostly lead to an underestimation of the forcing. Hence, because of the direct relation to the ice number concentration, one would expect underestimation of the n_post with increasing time step.

Uncertainties from the detection of the forcing term then lead to a variation in m_eff, considered in Fig. 10. The overall error we introduce by using the fit is illustrated by the first PDF when the same time step as that in the full system is used. The additional deviation introduced by differences in F(t₀) leads to a variation of about 2 %–3 % in the mean error of m_eff, which is negligible compared to the overall error introduced by the fit.

Figure 10 Probability density function of the relative error of m_eff predicted by fit (10) with respect to the exact mass calculated from the full system using (9) in %. Horizontal bars correspond to the 90th, 50th, and 10th percentiles of the distribution.

Comparing the error of the ice number concentration prediction from the exact solution in Fig. 11, one observes the increase in the errors when going to larger time steps. The mean and 90th percentile of the distribution stay nearly the same when the time step is increased from 1 to 60 s. The increase in mean error becomes larger for a time step dt=120 s, and, at a time step dt=240 s, it becomes closer to 50 %. We conclude that one can use a time step of 60 s without a significant error increase. To estimate the quality of the new parameterisation, one can compare the mean errors to the relative accuracy of measurements (see e.g. Krämer et al., 2016, 2020). The resulting mean error using a time step of 60 s is of the same order as the accuracy of measurements, which is about 15 % (cf. Luebke et al., 2016). Therefore, the proposed approach and chosen time step are well suited for the representation of the process.

Figure 11 Probability density function of the relative error of n_post predicted using fit (10) by (8) compared to the exact mass calculated from the full system in %. Horizontal bars correspond to the 90th, 50th, and 10th percentiles of the distribution.

In order to check the dependence of the statistics on the time step used in integrating Eqs. (5) and (6), the comparison of the PDFs of m_eff and n_post is shown in Fig. 12. The overall shape of the approximated PDF of mass m_eff agrees with the exact solution. However, the case with the largest time step of 2 min shows an underestimated m_eff minimum for several orders of magnitude. This leads to poor representation of the cases with higher n_post and consequent overestimation of the n_post maximum (see Fig. 12 lower panel). Predictions of the ice number concentrations are not far from the exact full system results if one compares the centre of the PDF. The peak of the predicted n_post PDFs stays in the range of 2 % difference compared to the full system for all the parameterisation cases. However, the prediction made using the largest time step substantially overestimates the maximum n_post value. It is worth mentioning that the ice number concentrations obtained with the parameterisation and full system are in the range of the observed quantities for the homogeneous nucleation (Krämer et al., 2020).

Figure 12 Probability density functions of m_eff,n_post, and q_i from the full model and predicted by the parameterisation (8) with different used time steps.

Furthermore, we present two nucleation event cases in order to directly compare the evolution and parameterisation prediction with different time steps (see Figs. 13 and 14). By means of such comparison, one can draw an additional conclusion on the correctness and possibility of use of q_i as mentioned earlier and assess the numerical stability of the method under different time steps. Initial conditions for the presented cases are picked from the ensemble calculation with six superposed GWs with the presence of a background updraft.

Figure 13 Comparison of the parameter evolution in full system, reduced system, and parameterisation with proposed correction. Initial conditions: S=1.4, n_i=50 kg⁻¹, $q_{\text{i}}=\mathrm{2.6}\times {\mathrm{10}}^{-\mathrm{13}}$ kg kg⁻¹.

Figure 14 Comparison of the parameter evolution in full system, reduced system, and parameterisation with proposed correction. Initial conditions: S=1.4, n_i=0.018 kg⁻¹, $q_{\text{i}}=\mathrm{6}\times {\mathrm{10}}^{-\mathrm{16}}$ kg kg⁻¹.

Neither case shows any difference between the reduced system (4)–(6) and the full system (1)–(3). The corrected parameterisation (8) is tested for larger time steps of 30, 120, and 240 s, omitting the case with the smallest time step. The first case in Fig. 13 illustrates that the parameterisation gives a good agreement with the full system, even with larger time steps. The evolution of the saturation ratio over ice, ice mixing ratio, and mass is aligned with the full system results after some time following the nucleation event. The prediction of the ice number concentration is not that exact for cases with larger time steps, but the difference is still in the range of 30 %. After the prediction of q_i using the m_eff and n_post from the parameterisation, one can see that integration is done without numerical problems. The evolution of the ice mixing ratio q_i shows a jump due to the nucleation process, followed by a short relaxation period, after which the q_i prediction is adjusted and overlaps with the full system result.

The second test case presented in Fig. 14 has been picked on purpose to show several possible outcomes of using the parameterisation and the consequences of usage of larger time steps. The evolution of the slowly varying parameters. such as $S,q_{\text{i}},m$. is captured relatively well by the parameterisation with time steps of 30 and 120 s. The largest considered time step, however, leads to smaller variations in the saturation ratio; hence the nucleation event is not captured by the parameterisation.

The parameterisation using a time step of 30 s gives a sufficiently good agreement with a relative error compared to the exact solution not larger than 22 %. However, even though the parameterisation with a 2 min time step predicts the nucleation event, the difference in the final ice number concentration is increased to 46 %. Other variables and the evolution of the slowly varying part are described well in cases where the nucleation is captured as such. Considering different cases and referring to the representation of the PDFs, a time step of 60 s can be suggested for parameterisation for the optimal result. Further increase in time step leads to a decrease in accuracy and, in cases of the largest time step of 4 min, even to a wrong representation of the process, where some events are not captured. Such behaviour is observed due to the fact that the pre-existing ice grows too fast, depleting the available water vapour and thus preventing the saturation ratio from increasing to values triggering the nucleation.

4 Conclusions

This study aimed to extend the applicability of the nucleation parameterisation of Dolaptchiev et al. (2023) to a broader range of atmospheric conditions. Specifically, it addressed the regime of low pre-existing ice, where the rapid evolution of the ice mixing ratio limits the validity of the asymptotic approach used in earlier work. To overcome this limitation, a generalisation of the parameterisation derived in Dolaptchiev et al. (2023), including the effect of variable ice mean mass, is proposed in the article. Compared to the constant mean mass assumption, the proposed approach has the advantage of incorporating a physically grounded correction to the growth coefficient, reflecting the actual variation in ice particle growth under different conditions. The considered method allows the inclusion of the mean mass variation, which leads to a better representation of the predicted ice number concentration distribution.

The new approach is validated against ensemble simulations based on the full double-moment scheme, showing a good agreement with the exact solution. The representation of individual nucleation events is done with reasonable accuracy, and the mean relative error with respect to the full system is less than 16 % for the ice number concentration, which is similar to the usual accuracy of measurements. A vast number of cases are represented accurately, but the substantial differences remain for the extreme cases of the high forcing-term values. The probability density function of the occurring n_i determined from independent ensemble simulations captures the major statistical characteristics, such as the mean and the centre of the PDF.

Because of the computational efficiency demand for NWPs and climate applications, the effect of uncertainties resulting from large ice microphysics time steps in the proposed extended parameterisation is assessed as well. The optimal representation of individual nucleation events combined with the sufficiently well captured statistics of n_post occurrence is achieved with dt=60 s. This recommendation is comparable to the fast physics time step of about 4 min, typically used in NWP models. Evaluation of the numerical effectiveness shows that parameterisation is several times faster than the solution of the full system with a time step of 1 s.

Given the strong agreement of the proposed method with resolved physics and its robustness and sufficient accuracy, the approach is now ready for further application in NWP models, such as ICON. The generalisation of the method allows direct coupling with forcings from dynamical parameterisations, thereby enhancing its integration into operational systems. For instance, the detailed representation of local subgrid-scale GW dynamical fields impacting cirrus clouds can be incorporated using the Multi-Scale Gravity Wave Model (MS-GWaM) parameterisation (Bölöni et al., 2021; Voelker et al., 2024) coupled to the described ice nucleation parameterisation. Additionally, this coupling approach can be extended to other relevant dynamical effects, e.g. turbulence.

Appendix A

Table A1 Values or ranges for the physical quantities within the relevant ranges for the UTLS region; see also Baumgartner and Spichtinger (2019) and Dolaptchiev et al. (2023).

Coupling of the ice physics processes to the GW theory is done under the assumptions of mid-frequency and high-frequency GWs, where the GW timescale T_w is tied to the Brunt–Väisälä frequency $N=\sqrt{\frac{g}{\mathit{\theta }}\frac{\mathrm{d}\stackrel{\mathrm{‾}}{\mathit{\theta }}}{\mathrm{d}z}}$ by the following relation:

$$T_{\text{w}}={\displaystyle \frac{\mathrm{1}}{{\mathit{\epsilon }}^{\mathit{\beta }}N}},\mathit{\beta }=\mathrm{0},\mathrm{1};\mathit{\epsilon }=\mathrm{0.1},$$

where β=0 characterises the high-frequency GW timescale, β=1 characterises the mid-frequency GW timescale, and ε is an introduced small parameter used for defining the distinguished limit in the asymptotic analysis. For instance, the characteristic timescales in the conditions of the troposphere, where $N={\mathrm{10}}^{-\mathrm{2}}$ s⁻¹, are T_w=100 and T_w=1000 s for the high- and mid-frequency waves, respectively. The diffusional growth timescale is estimated from $T_{\text{d}}=(D{m}_{\text{c}}^{\mathrm{1}/\mathrm{3}}\frac{p_{\text{si, c}}}{p_{\mathrm{00}}{q}_{\mathrm{v},\text{c}}}{T}_{\mathrm{00}}{n}_{\text{c}}{)}^{-\mathrm{1}}$ and has a value of about ∼ 340 s. Here, n_c is the reference ice number concentration, and m_c is the reference mean mass of ice particles. Given the characteristic timescales and expanding the Exner pressure π in terms of a small parameter ε, one can write the following expression for the derivative of π (in non-dimensional form; see Dolaptchiev et al., 2023):

$$\begin{array}{}\text{(A1)}& {\displaystyle \frac{\mathrm{D}\mathit{\pi }}{\mathrm{D}t}}={\mathit{\epsilon }}^{\mathrm{2}+\mathit{\alpha }}{\displaystyle \frac{\mathrm{D}{\mathit{\pi }}^{\prime }}{\mathrm{D}t}}+\mathit{\epsilon }w{\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{‾}}{\mathit{\pi }}}{\mathrm{d}{z}_s}},\end{array}$$

where $\mathit{\pi }={\left(\frac{p}{p_{\mathrm{00}}}\right)}^{\frac{R}{c_{\text{p}}}}$, α=0 corresponds to the strong stratification case, α=1 corresponds to the weak stratification case (Achatz et al., 2017, 2023), $\stackrel{\mathrm{‾}}{\mathit{\pi }}$ is the Exner pressure of the hydrostatically balanced reference atmosphere, R is the gas constant for dry air, and π^′ describes the fluctuations due to the wave field.

As shown in Dolaptchiev et al. (2023), applying the scaling analysis, one can neglect the first term on the right-hand side of Eq. (A1). In the case where the leading-order vertical velocity is solely due to a single GW, w in the second term of Eq. (A1) can be written as

$$w=\left|\widehat{w}{|}^{\left(\mathrm{0}\right)}\cos \right(\mathit{\omega }t+\mathit{\varphi }).$$

For a superposition of j gravity waves, the following expression applies:

$$w=\sum _j\left|{\widehat{w}}_j\right|\cos ({\mathit{\omega }}_{j}t+{\mathit{\varphi }}_j).$$

Here ${\widehat{w}}_j$, ω_j, and ϕ_j are the vertical wind amplitude, frequency, and phase of GW number j.

Substituting the following expressions in the Eq. (6) and also taking into account constant background updraft w₀₀, the last term can be written in the following form:

$$\begin{array}{}\text{(A2)}& F\left(t\right)={\displaystyle \frac{g{L}_{\text{i}}}{c_{\text{p}}{R}_{\text{v}}{T}^{\mathrm{2}}}}\left(\sum _j{\widehat{w}}_j\cos \right({\mathit{\omega }}_{j}t+{\mathit{\varphi }}_j)+w_{\mathrm{00}}),\end{array}$$

where L_i is the latent heat of sublimation, g is the acceleration due to gravity, c_p is the heat capacity, R_v is the gas constant for dry air, and T is the background reference temperature.

Appendix B

The parameter selection and limitations of the proposed fit are further evaluated using additional results obtained with a more advanced fit formulation for m_eff. This version incorporates a third parameter, m^*, and introduces a more complex functional dependency, expanding the simple fit to a third-order polynomial and increasing the number of fitting coefficients to 17. This extension is designed to better capture the spread observed in the data set shown before in Fig. 4.

The PDFs of the relative errors of the m_eff (mean radius of ice particle R₀) and n_post (predicted with the use of the more complex fit) with respect to the reference are shown in Fig. B1. As mentioned in Sect. 2.3, the improvement in the target parameter n_post is minor. The PDFs of predicted m_eff and n_post (see Fig. B2) indicate a better representation of the maximum and right tail of the PDF of ice number concentration. Despite those improvements, the use of a more comprehensive fit is not considered further, due to efficiency reasons and divergent behaviour of the polynomials in the extreme cases.

Figure B1 Probability density functions of the relative errors in the predictions of m_eff, R_eff, and n_post made with the advanced fit, with respect to the exact solutions. Vertical lines correspond to the 90th, 50th, and 10th percentiles of the distribution.

Figure B2 Probability density functions for m_eff, n_post, and q_i exact ensemble solutions and predictions made with the advanced fit. Vertical lines correspond to the mean of the distribution.

Appendix C

Figure C1 The ice number concentration calculated based on the proposed approach with the conditions of T= 220 K, a time step of 10 s, and the prescribed constant updraft (param.) and single GW forcing (param. GW forcing) is shown depending on the w or $\widehat{w}$. For comparison, the results based on the reference system (1)–(3) without sedimentation are shown for the GW forcing (ref. GW forcing). As a reference for the constant updraft, the values from detailed microphysics calculations from Kärcher and Lohmann (2002) are shown, where the upper/middle/lower symbols correspond to an initialisation of narrower/baseline/broader aerosol size distributions.

In order to further quantify the proposed parameterisation and test it for the larger vertical velocity perturbations, a set of additional calculations was performed using the constant updraft and extreme amplitudes of GW-induced vertical velocity. The results are compared to the commonly used homogeneous nucleation parameterisations, such as Kärcher and Lohmann (2002), and a reference system (1)–(3) with resolved ice physics. The ice number concentration prediction is presented depending on constant vertical updraft and $\widehat{w}$ vertical updraft created by the GW (see Fig. C1). The proposed parameterisation shows a good agreement with the traditional approaches and captures the larger n_i for the larger vertical updraft, which was also found by Kärcher and Lohmann (2002) using a bin model. Cases where vertical updraft was forced by the GWs also show the same behaviour as the constant updraft; however, the validation is done against the resolved version. The extended parameterisation overestimates n_i compared to the resolved physics, but the final predictions stay in the same order of magnitude and follow the expected physical trend. Notably, the larger deviations occur primarily in scenarios with extreme vertical velocity perturbations, up to 15 m s⁻¹, which represent rare and highly energetic events in the UTLS region. Despite this, the parameterisation successfully captures the qualitative and quantitative response of the system under such conditions. This demonstrates both the robustness and the flexibility of the approach, making it a reliable and efficient choice for representing a wide range of atmospheric scenarios, including those with strong dynamic forcing.