A new parameterization of ice heterogeneous nucleation coupled to aerosol chemistry in WRF-Chem model version 3.5.1: evaluation through ISDAC measurements

In the Arctic, during polar night and early spring, ice clouds are separated into two leading types of ice clouds (TICs): (1) TIC1 clouds characterized by a large concentration of very small crystals and TIC2 clouds characterized by a low concentration of large ice crystals. Using a suitable parameterization of heterogeneous ice nucleation is essential for properly representing ice clouds in meteorological and climate models and subsequently understanding their interactions with aerosols and radiation. Here, we describe a new parameterization for ice crystal formation by heterogeneous nucleation in water-subsaturated conditions coupled to aerosol chemistry in the Weather Research and Forecasting model coupled with chemistry (WRF-Chem). The parameterization is implemented in the Milbrandt and Yau (2005a, b) two-moment cloud microphysics scheme, and we assess how the WRF-Chem model responds to the run-time interaction between chemistry and the new parameterization. Well-documented reference cases provided us with in situ data from the spring 2008 Indirect and Semi-Direct Aerosol Campaign (ISDAC) over Alaska. Our analysis reveals that the new parameterization clearly improves the representation of the ice water content (IWC) in polluted or unpolluted air masses and shows the poor performance of the reference parameterization in representing ice clouds with low IWC. The new parameterization is able to represent TIC1 and TIC2 microphysical characteristics at the top of the clouds, where heterogenous ice nucleation is most likely occurring, even with the known bias of simulated aerosols by WRF-Chem over the Arctic.

1 Introduction

The Arctic is warming faster than the global mean, and projections for the future suggest that this tendency will continue (IPCC, 2013). The contribution of aerosols to the changing climate of the Arctic is poorly known. Aerosols perturb the radiative balance directly by absorbing radiation and indirectly due to aerosol effects on cloud properties. This leads to increases in shortwave scattering efficiency and infrared radiation (IR) emissivity alterations of Arctic clouds (Zhao and Garrett, 2015; Shindell and Faluvegi, 2009). The radiative properties and lifetime of clouds are particularly sensitive to aerosol concentration, composition and size. While the uncertainties associated with the indirect effects of aerosols on liquid clouds are still large, the effect of ice nucleation is even less well understood. Ice particle formation in tropospheric clouds significantly changes cloud microphysical properties, radiation balance and precipitation efficiency. At the core of the problem, ice nucleation causes multiple changes to cloud behavior, which at present are difficult to quantify. In its latest report, the IPCC (Intergovernmental Panel on Climate Change) was unable to estimate the radiative forcing of aerosols on clouds through ice nucleation (Boucher et al., 2013).

The detailed process of ice nucleation in cold clouds is complex and remains a major challenge for parameterization in atmospheric models. This is especially the case for polar ice clouds, for which the paucity of observations is a serious limitation (Curry et al., 1996; Kanji et al., 2017; McFarquhar et al., 2017). For instance, instead of assuming that cloud particles are distributed homogeneously, to investigate model response and climate sensitivity, some models have based their parameterization on in situ observations (Kay et al., 2016; Cirisan et al., 2020). However, the strong coupling between clouds and state variables, particularly temperature and moisture or relative humidity, requires a dynamic coupling of the cloud microphysics interactively with the atmospheric state variables. Among these coupling processes, the efficiency of ice-nucleating particles (INPs) to activate cloud formation is critical given the rarity of INPs in the pristine atmosphere. Two approaches are used to treat the INP efficiency: a singular and deterministic method, or a stochastic method (Pruppacher and Klett, 1997). While the singular approach assumes nucleation to occur at specific relative humidity and temperature (e.g., Wheeler and Bertram, 2012; Murray et al., 2012), the stochastic method allows for time-dependent state variables following the classical nucleation theory (CNT) (Pruppacher and Klett, 1997; Cirisan et al., 2020). It is also our approach in this study, whereby we assume that freezing occurs at any location on the INP surface with equal probability. This is one attempt to best represent in situ observations, yet it is still not fully physically comprehensive and one exploration step. The ultimate general method is still a matter of intense research (Vali, 2014; Wright and Petters, 2013).

Most atmospheric models use simple time-independent parameterizations of ice nucleation, predicting ice crystal number concentration either as a function of temperature (Fletcher, 1962; Cooper, 1986) or ice supersaturation (e.g., Meyers et al., 1992). These parameterizations do not include a limitation of ice crystal number concentration by the number of available ice nuclei particles and can lead to very poor estimation of ice crystal number concentration, in particular if they are applied outside the range of measurements used to constrain them (Prenni et al., 2007). This is particularly true for ice clouds in Arctic conditions (Keita and Girard, 2016). In the CNT model, a crucial fitting parameter is the contact angle (θ), quantifying the wettability of a solid particle surface by ice via the Young–Dupré equation. It is generally described as a single contact angle for an entire aerosol population, which does not work well for predicting the fractions of INPs on dust aerosol or on particles that have heterogeneous surfaces (Hoose and Möhler, 2012).

In recent years, with increasing data on ice nucleation from field and laboratory studies, new time-independent parameterizations have been developed, often based on empirical fits to atmospheric INP measurements as a function of temperature and aerosol particle size distributions (e.g., Connolly et al., 2013; Welti et al., 2012; Phillips et al., 2013; DeMott et al., 2010, 2015; Cirisan et al., 2020). Despite significant advances, they are of limited use in large-scale models operating over a wide range of temperatures. More complex CNT parameterizations than those using contact angle (θ probability density function or PDF) come at high computational costs (Welti et al., 2012; Murray et al., 2012; Niedermeier et al., 2014). In the particular context of climate simulations in Arctic atmospheric and chemical conditions, there is a need for efficient parameterizations of heterogeneous ice nucleation using simplified approaches to limit computational time.

In Keita et al. (2019), the parameterization of Girard et al. (2013) for water-subsaturated conditions based upon the CNT approach was implemented in the online Weather Research and Forecasting model coupled with chemistry (WRF-Chem) (Grell et al., 2005). This parameterization is suitable to represent the formation of ice clouds in the Arctic. It assumes that INPs are mainly mineral dust particles, which is consistent with recent results from the NETCARE (Network on Climate and Aerosols: Addressing Key Uncertainties in Remote Canadian Environments) project (Abbatt et al., 2019). This parameterization considered physicochemical properties of INPs, which are important in Arctic conditions, especially during winter and early spring (Eastwood et al., 2009; Keita and Girard, 2016) when sulfuric acid is often a dominant component of the aerosol, known as Arctic haze. Two types of ice clouds (TICs) were characterized Grenier et al. (2009). A TIC1 is an ice cloud seen by lidar but unseen by radar and is composed of a relatively large number of nonprecipitating small ice crystals; its ice crystal number concentration is higher than 10 L⁻¹. This cloud can have an upper part composed of low concentrated precipitating ice crystals. The second type, TIC2, is an ice cloud seen by radar and lidar characterized by a low concentration of larger precipitating ice crystals with an ice crystal number concentration lower than 10 L⁻¹. After spatial and temporal evaluation of the model, Keita et al. (2019) showed the ability of the parameterization to discriminate TIC1 and TIC2 clouds observed during the Indirect and Semi-Direct Aerosol Campaign (ISDAC) (McFarquhar et al., 2011). However, the study of Keita et al. (2019) was constrained by a prescribed concentration of aerosols with a fixed acid concentration.

In this paper, we investigate ice heterogeneous nucleation for the first time in a fully coupled aerosol and chemistry parameterization. We evaluate the response of the WRF-Chem model to the realistic time-dependent interaction between aerosols predicted by the chemistry module and the contact angle approach proposed by Girard et al. (2013). The new parameterization significantly improves the treatment of ice nucleation by discriminating TIC1 and TIC2 cloud formation as a function of the aerosol chemical composition. Each cloud is closely analyzed against observational data from three detailed flights conducted during ISDAC (2008). This study is part of the NETCARE project addressing key uncertainties in remote Canadian environments with the objective of assessing the impact of aerosols on Arctic ice clouds.

The paper is organized as follows. Section 2 briefly describes the Milbrandt and Yau (2005a, b) scheme for cloud microphysics and the presentation of ice heterogeneous nucleation parameterization coupled with aerosol chemistry. Section 3 presents the test cases from ISDAC and Sect. 4 the evaluation of the new parameterization against ISDAC. Section 5 is dedicated to the conclusion.

2 Description of the new scheme for ice heterogeneous nucleation in WRF-Chem

The new scheme for ice crystal formation by heterogeneous nucleation in the deposition mode is implemented in WRF-Chem version 3.5.1. WRF-Chem is a regional, fully coupled “online” model (Grell et al., 2005), for which all prognostic meteorological, chemical and aerosol variables are fully integrated within WRF-ARW, a mesoscale meteorological model; it uses the same grid, time step, advection scheme and physics schemes as WRF-ARW. Several schemes are available in WRF-Chem for cloud microphysics. We choose the Milbrandt and Yau (2005a, b), MY05, for its ability to simulate Arctic clouds in previous works (Keita et al., 2019; Keita and Girard, 2016).

2.1 Overview of the two-moment version of the cloud microphysical scheme MY05

MY05 (Milbrandt and Yau, 2005a, b) is a bulk cloud microphysics parameterization with one-, two- and three-moment versions. We use the two-moment version available in WRF-Chem. It includes the following prognostic variables: the mass mixing ratio q_x and the number concentration N_x, with $x\in (c,r,i,s,h,g)$ respectively representing cloud liquid water (c), cloud ice water (i), rain (r), snow (s), hail (h) and graupel (g). The time evolutions of the hydrometeor mass mixing ratio and number concentration are respectively governed by the following prognostic equations:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\displaystyle \frac{\partial q_x}{\partial t}}=& -{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}\mathrm{\nabla }\cdot \left(\mathit{\rho }{q}_x\mathit{U}\right)+\mathrm{\nabla }\cdot \left(\mathbf{K}\mathrm{\nabla }{q}_x\right)\\ & +{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\rho }{q}_x{V}_{Qx}\right)+{\displaystyle \frac{\mathrm{d}{q}_x}{\mathrm{d}t}}{|}_{\mathrm{s}},\end{array}\end{array}$$

and

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\displaystyle \frac{\partial N_{\mathrm{T},x}}{\partial t}}=& -\mathrm{\nabla }\cdot \left(N_{\mathrm{T},x}\mathit{U}\right)+\mathrm{\nabla }\cdot \left(\mathbf{K}\mathrm{\nabla }{N}_{\mathrm{T},x}\right)\\ & +{\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\rho }{N}_{\mathrm{T},x}{V}_{Nx}\right)+{\displaystyle \frac{\mathrm{d}{N}_{\mathrm{T},x}}{\mathrm{d}t}}{|}_{\mathrm{s}},\end{array}\end{array}$$

where ρ is the density of air, U the 3D velocity vector, V_Qx the mass-weighted fall speed, N_T,x the total number concentration per unit volume, V_Nx the number-weighted fall speed and K the turbulent diffusion matrix. The right-hand side terms of both equations respectively represent advection and divergence, turbulent mixing, sedimentation, and microphysical tendencies (marked by the “s” subscript).

The mass of a single hydrometeor for the x category is parameterized as a power law of the form

$$\begin{array}{}\text{(3)}& m_x\left(D\right)=c_x{D}^{\mathrm{3}},\end{array}$$

where $c_x={\mathit{\rho }}_x\frac{\mathit{\pi }}{\mathrm{6}}$, with ρ_x the bulk density (Table 1) for spherical particles x (cloud liquid water, rain, snow, graupel and hail). Cloud ice crystals are assumed to be bullet rosettes (Schoenberg Ferrier, 1994) with c_i=440 kg m⁻³. The size spectrum of each category is described by a common generalized gamma distribution function (Schoenberg Ferrier, 1994) of the form

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{N}_x\left(D\right)=& N_{\mathrm{T},x}{\displaystyle \frac{{\mathit{\nu }}_x}{\mathrm{\Gamma }(\mathrm{1}+{\mathit{\alpha }}_x)}}{\mathit{\lambda }}_x^{{\mathit{\nu }}_x(\mathrm{1}+{\mathit{\alpha }}_x)}\\ & D^{{\mathit{\nu }}_x(\mathrm{1}+{\mathit{\alpha }}_x)-\mathrm{1}}\exp \left(-({\mathit{\lambda }}_{x}D{)}^{{\mathit{\nu }}_x}\right),\end{array}\end{array}$$

where N_x(D) is the number concentration of hydrometeor x per unit volume per unit diameter D, α_x is the shape parameter controlling the size dispersion, λ_x is the slope and ν_x is a second size dispersion parameter. The size distribution of cloud droplets is represented in MY05 by α_x=1 and λ_x=3. For all other hydrometeors ν_x=1, leading to the form

$$\begin{array}{}\text{(5)}& N_x\left(D\right)=N_{\mathrm{0}x}{D}^{{\mathit{\alpha }}_x}\exp (-{\mathit{\lambda }}_{x}D),\end{array}$$

where N_0x is the intercept parameter given by

$$\begin{array}{}\text{(6)}& N_{\mathrm{0}x}=N_{\mathrm{T},x}{\displaystyle \frac{\mathrm{1}}{\mathrm{\Gamma }(\mathrm{1}+{\mathit{\alpha }}_x)}}{\mathit{\lambda }}_x^{(\mathrm{1}+{\mathit{\alpha }}_x)}.\end{array}$$

Table 1 Bulk density for each hydrometeor category.

The four ice-phase hydrometeors follow the size distribution above. The cloud ice water category represents pristine ice crystals. The snow category includes crystals with radii greater than 100 µm and aggregates. The graupel category includes moderate-density graupel formed from heavily rimed ice or snow. The hail category corresponds to high-density hail and frozen raindrops. For each ice-phase hydrometeor x, the total number concentration N_T,x (kg⁻¹) and the mass mixing ratio q_T,x (kg kg⁻¹) are respectively given by

$$\begin{array}{}\text{(7)}& N_{\mathrm{T},x}=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{N}_{\mathrm{0}x}{D}^{{\mathit{\alpha }}_x}\exp (-{\mathit{\lambda }}_{x}D)\mathrm{d}D,\end{array}$$

and

$$\begin{array}{}\text{(8)}& q_{\mathrm{T},x}=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{m}_x\left(D\right)N_{\mathrm{0}x}{D}^{{\mathit{\alpha }}_x}\exp (-{\mathit{\lambda }}_{x}D)\mathrm{d}D,\end{array}$$

where m_x(D) is obtained from Eq. (3).

Table 2 Source and sink terms listed according to the hydrometeor category, which gains mass and number, except for self-collections or when the loss is to water vapor.

Microphysical processes represented in MY05 are summarized in Table 2, where processes are listed according to the hydrometeor category. The source and sink terms for the two-moment (mass content) scheme are from previous studies (Kong and Yau, 1997; Schoenberg Ferrier, 1994) and depend on the size distribution function. The primary sources of ice crystals in the atmosphere are heterogeneous and homogeneous ice nucleation. Homogeneous freezing is the spontaneous freezing of a water (or haze) droplet. According to Pruppacher and Klett (1997), the homogeneous freezing rate of cloud droplets is dominant at temperatures below $\sim -\mathrm{32}$ ^∘C. In the range −30 to −50 ^∘C, MY05 follows DeMott et al. (1994) with

$$\begin{array}{}\text{(9)}& \mathrm{\Delta }{N}_{\mathrm{freeze}}=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}(\mathrm{1}-\exp (-JV\mathrm{\Delta }t\left)\right)N_{T_{\mathrm{c}}}\left(D\right)\mathrm{d}D.\end{array}$$

In a given time step Δt, ΔN_freeze is the number of droplets freezing homogeneously and J is the nucleation rate for pure water. For homogeneous nucleation,

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{\log }_{\mathrm{10}}\left(J\right)& =-\mathrm{606.3952}-\mathrm{52.6611}{T}_{\mathrm{c}}-\mathrm{1.7439}{T}_{\mathrm{c}}^{\mathrm{2}}\\ & -\mathrm{2.65}\times {\mathrm{10}}^{-\mathrm{2}}{T}_{\mathrm{c}}^{\mathrm{3}}-\mathrm{1.536}\times {\mathrm{10}}^{-\mathrm{4}}{T}_{\mathrm{c}}^{\mathrm{4}},\end{array}\end{array}$$

with the volume V approximated by the mean droplet diameter in centimeters. Therefore, the fraction of cloud droplets freezing in one time step may be written as

$$\begin{array}{}\text{(11)}& F_{\mathrm{freeze}}={\displaystyle \frac{\mathrm{\Delta }{N}_{\mathrm{freeze}}}{N_{T_{\mathrm{c}}}}}\left(\mathrm{1}-\exp (-J{\displaystyle \frac{\mathit{\pi }}{\mathrm{6}}}{D}_{\mathrm{mc}}^{\mathrm{3}}\mathrm{\Delta }t)\right),\end{array}$$

where D_mc is the mean volume diameter of cloud droplets. Heterogeneous ice nucleation needs INPs, a minor fraction of the tropospheric aerosol, which exhibits micro surface structures to facilitate the formation of ice crystals. In the presence of INPs, if thermodynamic conditions are favorable, ice crystals can form by heterogeneous nucleation through four different modes. Deposition nucleation and condensation freezing can occur without the presence of supercooled droplets. For clouds below 0 ^∘C primarily composed of supercooled liquid droplets, ice crystal can form by immersion and contact freezing. This conceptual definition of heterogeneous ice nucleation (Pruppacher and Klett, 1997) is used in MY05. Contact freezing follows Young (1974) wherein the number concentration of contact INPs is a function of temperature according to Meyers et al. (1992). In the contact-freezing formation mode, ice nucleation occurs on a solid particle colliding with a supercooled liquid droplet. Immersion freezing of raindrops and cloud water droplets follows the parameterization of Bigg (1953). The deposition mode involves the growth of ice directly from the vapor phase, whereas condensation freezing occurs if the ice phase is formed immediately after condensation of water vapor on a solid particle as liquid intermediate. In the original version of MY05, deposition and condensation freezing are functions of water vapor supersaturation with respect to ice, S_i, following Meyers et al. (1992):

$$\begin{array}{}\text{(12)}& N_{m,\mathrm{i}}\left(S_{\mathrm{i}}\right)=\mathrm{1000}\exp \left[\mathrm{12.96}(S_{\mathrm{i}}-\mathrm{1})-\mathrm{0.639}\right],\end{array}$$

where N_m,i is the number of ice crystals predicted per unit volume due to deposition and condensation freezing. The Meyers et al. (1992) parameterization for deposition and condensation freezing depends only on supersaturation. It was derived from ground-based measurements. These approximations may lead to an overestimation of N_m,i when the number concentration of particles acting as INPs is low, such as in Arctic conditions (Eidhammer et al., 2009). Moreover, the immersion-freezing mode from Pruppacher and Klett (1997) has been extended to include freezing of immersed INPs inside an aqueous solution or wet aerosol (Vali et al., 2015), which is a significant process of Arctic ice cloud formation (Eastwood et al., 2008).

2.2 A new parameterization of ice heterogeneous nucleation coupled with chemistry for MY05 in WRF-Chem

The new parameterization focuses on heterogeneous ice nucleation for uncoated INPs and for sulfuric-acid-coated INPs in the deposition mode, i.e., in water-subsaturated conditions. In this approach, INPs are assumed to be mineral dust particles following Girard et al. (2013). For contact freezing and immersion freezing from supercooled cloud droplets, the parameterizations remain unchanged. As condensation freezing is uncertain (Vali et al., 2015), this process is no longer included in the model. The modified version of MY05 including our new parameterization described below is hereafter referred to as MYKE.

The parameterization is based on the CNT, a stochastic approach in which the nucleation rate J_d depends on the contact angle between an ice embryo and its INPs. Following CNT, in each time step Δt the number concentration of nucleated ice crystals N_f is given by

$$\begin{array}{}\text{(13)}& N_f\left(\mathrm{\Delta }t\right)=N_{\mathrm{t}}\exp \left[\mathrm{1}-J_d{A}_{\mathrm{d}}\mathrm{\Delta }t\right],\end{array}$$

where A_d is the total surface area of dust particles and N_t is the total number concentration of available INPs. In previous studies using this approach (Keita and Girard, 2016; Keita et al., 2019; Girard et al., 2013; Khvorostyanov and Curry, 2009; Morrison et al., 2005b; Liu et al., 2007; Hoose et al., 2010; Chen et al., 2008), A_d and N_t were prescribed and constant over time, although the concentration of atmospheric INPs varied tremendously in time and space, as well as in their composition and origins. The new MYKE parameterization within WRF-chem now considers the temporal and spatial variation of A_d and N_t. J_d, the nucleation rate of embryos per unit surface of particles (Pruppacher and Klett, 1997; Martin, 2000; Hung et al., 2003; Pant et al., 2004; Parsons et al., 2004b; Archuleta et al., 2005; Pant et al., 2006), is defined as

$$\begin{array}{}\text{(14)}& J_d=B\exp \left({\displaystyle \frac{-\mathrm{\Delta }{G}^{\ast }}{kT}}\right),\end{array}$$

where B=10²⁶ ${\mathrm{cm}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ is the kinetic coefficient (Pruppacher and Klett, 1997), k is the Boltzmann constant (J K⁻¹), T is the temperature in Kelvin, and ΔG^∗ is the critical Gibbs free energy for the formation of an ice embryo (J) and is defined as

$$\begin{array}{}\text{(15)}& \mathrm{\Delta }{G}^{\ast }={\displaystyle \frac{\mathrm{16}\mathit{\pi }{\mathit{\sigma }}_{\mathrm{iv}}^{\mathrm{3}}f\left(\cos \mathit{\theta }\right)}{\mathrm{3}{\mathit{\rho }}_{\mathrm{i}}^{\mathrm{2}}{R}_{\mathrm{v}}^{\mathrm{2}}{T}^{\mathrm{2}}{\ln }^{\mathrm{2}}{S}_{\mathrm{i}}}},\end{array}$$

where ${\mathit{\sigma }}_{\mathrm{iv}}=\mathrm{106.5}\times {\mathrm{10}}^{-\mathrm{3}}$ J m⁻² is the surface tension between ice and water vapor, ρ_i=0.5 g cm⁻³ is the bulk ice density and R_v=461.5 $\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ is the gas constant for water vapor. The function f(cos θ) is a monotonic decreasing function of the cosine of the contact angle θ as defined by Pruppacher and Klett (1997) for a curved substrate:

$$\begin{array}{}\text{(16)}& \begin{array}{rl}f\left(\cos \mathit{\theta }\right)=& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left\{\mathrm{1}+{\left({\displaystyle \frac{\mathrm{1}-q\cos \mathit{\theta }}{\mathit{\varphi }}}\right)}^{\mathrm{3}}+q^{\mathrm{3}}\left[\mathrm{2}-\mathrm{3}\left({\displaystyle \frac{q-\cos \mathit{\theta }}{\mathit{\varphi }}}\right)\right.\right.\\ & \left.\left.+{\left({\displaystyle \frac{q-\cos \mathit{\theta }}{\mathit{\varphi }}}\right)}^{\mathrm{3}}+\mathrm{3}{q}^{\mathrm{2}}\cos \mathit{\theta }\left({\displaystyle \frac{q-\cos \mathit{\theta }}{\mathit{\varphi }}}-\mathrm{1}\right)\right]\right\},\end{array}\end{array}$$

where $\mathit{\varphi }=\sqrt{\mathrm{1}-\mathrm{2}q\cos \mathit{\theta }+q^{\mathrm{2}}}$ and $q=\frac{r_n}{r_{\mathrm{g}}}$, with r_g being the critical germ size expressed as

$$\begin{array}{}\text{(17)}& r_{\mathrm{g}}={\displaystyle \frac{\mathrm{2}{\mathit{\nu }}_{\mathrm{w}}{\mathit{\sigma }}_{\mathrm{iv}}}{kT\ln \left(S_{\mathrm{i}}\right)}},\end{array}$$

where ν_w is the volume of a water molecule.

In the CNT, the contact angle θ is a very important variable because it represents the ability of an INP to form ice. The lower the contact angle, the better an INP the aerosol is. Numerous laboratory studies have found realistic values of θ based on the physicochemical composition of aerosols (e.g., Marcolli et al., 2007; Eastwood et al., 2008; Fornea et al., 2009; Welti et al., 2009; Kanji and Abbatt, 2010; Welti et al., 2009). The CNT approach using these values was subsequently applied successfully in climate and forecast models at different scales (Khvorostyanov and Curry, 2009; Morrison et al., 2005a; Liu et al., 2007; Chen et al., 2008). For example, using the parameterization of Girard et al. (2013) based on laboratory studies from Eastwood et al. (2008, 2009), Keita et al. (2019) were able to simulate Arctic clouds forming in polluted and clean air masses with a prescribed contact angle of 26 and 12^∘, respectively. These studies were, however, limited because the contact angles represent extreme cases that must be prescribed arbitrarily before the simulation, and they assumed homogeneity of the degree of acidity of clouds in space and time throughout the whole domain.

For the first time, a real-time variable contact angle is used here in the CNT approach by coupling MY05 with the chemical module in WRF-Chem. This coupling is between MY05 and the MOSAIC (Model for Simulating Aerosol Interactions and Chemistry) aerosol module (Zaveri et al., 2008). MOSAIC simulates a wide variety of aerosol species: sulfates, methanesulfonate, nitrate, chloride, carbonate, ammonium, sodium, calcium, black carbon (BC), primary organic mass (OC), liquid water and other inorganic mass (OIN). OIN represents unspecified inorganic species such as silica (SiO₂), other inert minerals and trace metals lumped together assimilated to mineral dust. MOSAIC uses a sectional approach to represent aerosol size distributions by dividing the size distribution for each species into several size bins (four or eight available in WRF-Chem) and assumes that the aerosols are internally mixed in each bin. MOSAIC considers the major aerosol processes of inorganic aerosol thermodynamic equilibrium, binary aerosol nucleation, coagulation and condensation but does not include secondary organic aerosol (SOA) formation in the version used in this study. MOSAIC is a good compromise between accuracy and computing performance. It is used in WRF-Chem with four chemical mechanisms.

The coupling is done by expressing θ as a function of the aerosol neutralization fraction f_n in dust particles internally mixed with sulfate, nitrate and ammonium (Zhang et al., 2007; Fisher et al., 2011), which is between 0 and 1 and is defined as

$$\begin{array}{}\text{(18)}& f_{\mathrm{n}}={\displaystyle \frac{\left[{\mathrm{NH}}_{\mathrm{4}}^{+}\right]}{\mathrm{2}\left[{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}\right]+\left[{\mathrm{NO}}_{\mathrm{3}}^{-}\right]}}.\end{array}$$

This was motivated by several previous studies (Jouan et al., 2012; Grenier and Blanchet, 2010; DeMott et al., 2010; Blanchet and Girard, 1994; Keita et al., 2019; Keita and Girard, 2016) suggesting that the acidification of ice nuclei by the oxidation of sulfur dioxide forming sulfuric acid in the Arctic greatly alters the microphysical response of ice clouds. Such ice clouds tend to have bigger and fewer ice crystals than ice clouds formed in pristine environments. For instance, Kulkarni et al. (2014) showed that, except for quartz, acid-coated dust makes less effective INPs in the deposition mode but has similar effectiveness in the immersion-freezing mode, i.e., in water-supersaturated regime. Based on X-ray diffraction analyses, they argued that acid treatment caused structural deformations of the surface dust, and the lack of structured order reduced the ice nucleation properties of coated particles in the deposition mode. Moreover, they suggested that, at water-supersaturated conditions, surface chemical reactions might not change the original ice-nucleating properties permanently because coating material could be removed by dissolution. Panda et al. (2010) concluded that sulfuric-acid-treated kaolinite particles could result in the formation of aluminum sulfate that can easily be dissolved in water. Considering these recent findings and our objective to develop a simplified parameterization to limit computational time, we choose to use the CNT formula for deposition mode but with a specific factor, the neutralization fraction f_n, indicating the degree of acidity of the coating of dust particles.

Moreover, θ has been derived by Eastwood et al. (2008, 2009) from heterogeneous nucleation rates on kaolinite particles obtained in laboratory measurements. As a best fit, they found limiting values of $\mathit{\theta }=\mathrm{26}{}^{\circ }$ in polluted air and $\mathit{\theta }=\mathrm{12}{}^{\circ }$ in clean air. Kaolinite represents a significant component of mineral dust (Glaccum and Prospero, 1980). It is also found to be an efficient ice nucleus in the deposition mode, requiring relative humidity with respect to ice (RH_i) below 112 % in order to initiate ice crystal formation (Eastwood et al., 2009). This is a typical microphysical condition found in Arctic ice clouds. Recent studies from Kumar et al. (2018, 2019a, b) showed that

1. the relevance of quartz particles as atmospheric INPs is uncertain,

2. INP activity of dust particles not only depends on their composition but also on their chemical exposure history, and

3. the exposition of dust particles to acidic air masses decreases their INP activity.

Thus, using kaolinite as a proxy for dust particles in our parameterization is reasonable in the current state of knowledge on dust particle composition in the atmosphere, in particular in the Arctic atmosphere where our parameterization applies. Keita and Girard (2016), after analyzing the slope between the nucleation rate and the saturation over ice for TIC1 and TIC2 clouds (see Fig. 16 in Keita and Girard, 2016), observed for a given S_i that

1. the slope is the largest for the smallest accessible contact angle, and

2. the decrease in the slope with the increasing contact angle is very nonlinear.

These results are consistent with laboratory experiments (Sullivan et al., 2010) showing a rapid increase in the contact angle with acidity on coated INPs. These results motivated us to parameterize the contact angle θ as a function of the aerosol neutralization fraction f_n under a concave form. Simple concave functions follow the power law:

$$\begin{array}{}\text{(19)}& \mathit{\theta }=\mathrm{26}-\mathrm{14}{f}_{\mathrm{n}}^p,\end{array}$$

with p larger than 1. We have chosen a quadratic (p=2) form for simplicity:

$$\begin{array}{}\text{(20)}& \mathit{\theta }=\mathrm{26}-\mathrm{14}{f}_{\mathrm{n}}^{\mathrm{2}}.\end{array}$$

We have also added a sensitivity simulation under a biquadratic form (p=4) to test the influence of the exponent p on the concave form of the contact angle with the neutralization fraction:

$$\begin{array}{}\text{(21)}& \mathit{\theta }=\mathrm{26}-\mathrm{14}{f}_{\mathrm{n}}^{\mathrm{4}}.\end{array}$$

Both formulations, referred to as MYKE2 (Eq. 20) and MYKE4 (Eq. 21), are implemented in MY05 and tested hereafter. They imply that θ is close to 26^∘ for $\mathrm{0}<f_{\mathrm{n}}<\mathrm{0.5}$ with a more (Eq. 21) or less (Eq. 20) rapid decrease between 0.5 and 1 as shown in Fig. 1. The coupling between MY05 and MOSAIC is done by taking information from MOSAIC for A_d and N_t as needed to compute Eq. (13) and for f_n to compute Eqs. (20) and (21). These parameters are computed assuming the same aerosol size bin definition as in MOSAIC.

Figure 1 Variation of f_n with θ for MYKE2 (blue line) and MYKE4 (green line).

3 Configuration of the model for typical TIC1 and TIC2 clouds observed during ISDAC

ISDAC took place during April 2008 at the North Slope of Alaska. The objective was to study the role of Arctic aerosols in cloud microphysical properties and in the surface energy budget. Numerous studies have been based upon data from ISDAC (McFarquhar et al., 2011; Lawson et al., 2019). Among them, several studies investigated detailed parameters of ice clouds by analyzing the ISDAC database (Jouan et al., 2012; DeMott et al., 2010) or by running atmospheric models on case studies highlighted during the campaign (Keita and Girard, 2016; Matrosov et al., 2019; Keita et al., 2019). For instance, Keita et al. (2019) analyzed microphysical properties of TICs for ISDAC flights in nonpolluted and polluted environment using WRF simulations. Flights F13, F21 and F29 studied by Keita et al. (2019) were typical of a TIC1 cloud (F13) formed in a pristine air mass and of two TIC2 representative cloud cases (F21 and F29) formed in a polluted air mass. Here, our goal is to show the potential of the new ice nucleation parameterization to discriminate TIC1 and TIC2 cloud formation as a function of the aerosol chemical composition. Each cloud type is closely investigated using detailed observations from three flights conducted during ISDAC.

The simulations with WRF-Chem including MYKE are done over the whole period of ISDAC (McFarquhar et al., 2011), from 1 to 30 April 2008, on the domain shown in Fig. 2, which is identical to that described by Keita et al. (2019). The three test cases (F13, F21 and F29) are included in this period. The domain is based on a Lambert projection centered on Barrow, Alaska (now Utqiagvik), over 160 × 100 grid cells with a horizontal resolution of 10 km and 55 vertical levels between the surface and 50 hPa. The first 4 d of the simulation (1 to 4 April included) are used for model spin-up. Three simulations are performed: the first one uses the original MY05 scheme (the REF simulation), the second one uses the new parameterization given in Eq. (20) (the MYKE2 simulation) and the third one uses the new parameterization described by Eq. (21) (the MYKE4 simulation). WRF-Chem options and parameterizations used in these simulations are summarized in Table 3. As in Keita et al. (2019), meteorological initial and boundary conditions use NCEP (National Centers for Environmental Prediction) Global Forecast System (GFS) Final Analysis (FNL) data ($\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$), and the simulations are nudged to GFS FNL updated every 6 h above the planetary boundary layer (PBL).

Figure 2 Model domain (yellow) used in this study centered over Fairbanks with a horizontal resolution of 10 km. The cities of Barrow (now Utqiagvik; 71.18^∘ N,156.44^∘ E) and Fairbanks (64.83^∘ N,147.77^∘ E) over which the F12, F13, F21 and F29 flights took place are also shown with orange dots.

(Milbrandt and Yau, 2005a, b)(Iacono et al., 2008)(Iacono et al., 2008)(Berg et al., 2015)(Janjić, 1994)(Janjić, 1994)(Chen and Dudhia, 2001)(Zaveri et al., 2008)(Zaveri et al., 2008)(Wild et al., 2000)

Table 3 Parameterizations and options used for the WRF-Chem simulations.

For the chemical module, the CBM-Z (Carbon Bond Mechanism) photochemical mechanism (Zaveri and Peters, 1999) coupled with MOSAIC is used. CBM-Z has 67 species and 164 reactions in a lumped structure approach that classifies organic compounds according to their internal bond types. Rates for photolytic reactions are derived using the Fast-J photolysis rate scheme (Wild et al., 2000). Eight size bins are used in MOSAIC. Chemical initial and boundary conditions are taken from the global chemical-transport model MOZART-4 (Model for OZone And Related chemical Tracers, version 4) (Emmons et al., 2010). The fire emissions inventory used is the Fire INventory from NCAR (FINN-v1) (Wiedinmyer et al., 2011). FINN-v1 provides emissions on a per fire basis based on event count information from the Moderate Resolution Imaging Spectrometer (MODIS). The anthropogenic emissions come from the inventory developed within the POLARCAT Model Intercomparison Model Project (POLMIP), which includes SO₂ from both eruptive and noneruptive continuous degassing volcanism (Fisher et al., 2011; Jouan et al., 2014). During winter and spring 2008, sustained eruptive activity was recorded at Kamchatka and the Aleutian Islands (Fisher et al., 2011; Jouan et al., 2014; Atkinson et al., 2013; Burton et al., 2012). Noneruptive activity was common throughout our simulation period (Fisher et al., 2011; Jouan et al., 2014; Atkinson et al., 2013). Soil-derived (dust) and sea salt aerosol emissions are computed online into WRF-Chem respectively based upon the wind erosion formulation of Shaw et al. (2008) and the GOCART (Global Ozone Chemistry Aerosol Radiation and Transport model) sea salt emission module (Chin et al., 2000). For biogenic emissions, the Model of Emissions of Gases and Aerosols from Nature (MEGAN) (Guenther, 2007) computes them online using characteristics of the surface (class of vegetation, soil humidity and temperature, for instance).

4 Results and discussion

This section presents comparisons of WRF-Chem simulations (REF, MYKE2 and MYKE4) against observations, followed by a discussion of the results. Although the comparison between simulated results and observations is presented in the following along the entire vertical profile inside the clouds, the discussion focuses on the altitudes above the 500 hPa level, where heterogeneous nucleation is the most important process. According to Jouan et al. (2012), most of the differences between TIC1 and TIC2 events were confined at cloud top where ice nucleation mostly occurs and air is supersaturated with respect to ice. To compare simulations with observations along the ISDAC flight tracks, simulated results are averaged in a grid box of 10 by 10 km² centered on the location of the flight. ISDAC in situ measurements have been averaged every 20 s, corresponding to a vertical resolution of ∼45 hPa (∼450 m), during ascents and descents through clouds. Simulated WRF outputs are linearly interpolated to the pressure levels of these observations and temporally averaged over a 3 h period, encompassing the area of ISDAC flights. Some statistics are computed using the same method. First, we present some meteorological and chemical properties, followed by an analysis of cloud microphysical properties.

4.1 Temperature and relative humidity over ice

Table 4 presents biases (Bias), Pearson correlation coefficients (Cor) and root mean square errors (RMSE) for the temperature T and relative humidity over ice RH_i for the three simulations (REF, MYKE2 and MYKE4) and above the 500 hPa level. According to Jouan et al. (2012), the uncertainties on the measurements are estimated at ±0.5 ^∘C for T and ±11 % for RH_i. Note that vertical profiles of T and RH_i for flights F13, F21 and F29 are very close to results obtained by Keita et al. (2019). As expected, due to the nudging, the new heterogeneous ice nucleation parameterization does not significantly impact T and RH_i. The lowest temperatures at the top of the clouds, where the process of heterogeneous ice nucleation is important, are relatively well reproduced by MYKE2 and MYKE4 simulations with similar statistics (Cor ≃0.99, RMSE ≃2, Bias $\simeq -\mathrm{2}$ ^∘C), except along F21 flight (Cor ≃0.82, RMSE ≃3.3, Bias $\simeq -\mathrm{3}$ ^∘C), for which the observed increase in temperature caused by the heat exchanged at cold temperatures is not adequately represented by the model. For that flight, the three simulations underestimate RH_i by ±50 % at the top of the cloud. These biases are consistent with the large-scale GFS FLN fields and result in an underestimation of the altitude of the top of the cloud by the model for F21.

Table 4 Root mean square error (RMSE), bias (Bias) and Pearson correlation coefficients (Cor) of the temperature (T) and relative humidity over ice (RH_i) for the three simulations (REF, MYKE2 and MYKE4).

4.2 Aerosol properties

Figure 3 shows the comparison between observed and simulated (REF, MYKE2 and MYKE4) vertical profiles of total aerosol number concentrations (N_a). The Passive Cavity Aerosol Spectrometer Probe (PCASP) externally mounted under a wing of the Convair-580 aircraft sampled ambient clear air just before entering the cloud regions for all flights except F21. The optical particle counter (PCASP) provided particle size distributions and number concentrations in the geometric diameter size range 0.12–3 µm. To allow a fair comparison between WRF-Chem-simulated and PCASP-measured N_a, the model concentrations are summed over bins 3 to 6, corresponding to sizes between 0.156 and 2.5 µm. According to Shantz et al. (2014), the uncertainty in number concentration measured by the PCASP is approximately 10 %. First, the model does not reproduce the observed vertical variability. This may be due to the small sampling domain and time taken during ISDAC, which make comparisons between model simulations and the observed variability difficult, especially at the low horizontal resolution of 10 km used here. For F13, the air mass is relatively clean, with a weak vertical variability of aerosol number concentrations remaining mostly below 210 cm⁻³ on the whole column with mean concentrations around 73 cm⁻³, very close to the simulation mean of 86 cm⁻³. For F29, the PCASP shows that there is a much higher concentration of aerosol particles in the lower troposphere (more than twice that observed during F13, e.g., larger than 400 cm⁻³), particularly at altitudes above 550 hPa near cloud top where peak concentrations exceeding 1000 cm⁻³ have been measured. Comparing the two flights, between 550 hPa and 400 hPa, the simulated aerosol number concentration is overestimated by a factor of 3 above observations for flight F13 and underestimated by 1 order magnitude for flight F29 (Fig. 3). These discrepancies are consistent with Mölders et al. (2011), who analyzed aerosol concentrations during polar night around Fairbanks and showed an overestimation of aerosol concentrations over the nonpolluted site and an underestimation at an polluted site by using WRF-Chem. They concluded that discrepancies result from uncertainty in emissions, especially at Fairbanks. While most models agree that Arctic aerosols can be attributed to a mixture of anthropogenic sources, mesoscale models have difficulty properly simulating aerosol concentrations over the Arctic (Shindell et al., 2008; Eckhardt et al., 2015; Schwarz et al., 2013; Raut et al., 2017). Moreover, even if the simulated results show the same order of magnitude for N_a above 550 hPa (Fig. 3), whereas observations show a large difference between the two flights, we expect that the differences between simulated results for cloud microphysical properties for these two flights could be mainly explained by a combination of differences in the physicochemical properties of aerosols and the altitude of the simulated cloud top.

Figure 3 Comparison of the observed (red) and simulated (green) WRF vertical profiles of total aerosol number concentrations. Observations were measured by the PCASP in situ sensor onboard the Convair-580 just before entering the clouds for flights F13 (solid lines) and F29 (solid lines with diamond markers) Note that PCASP measurements were not available during flight F21.

Figure 4 presents simulated (REF, MYKE2 and MYKE4) vertical profiles of sulfate (SO₄), ammonium (NH₄) and nitrate (NO₃) molar aerosol concentrations along flights F13, F21 and F29, respectively. Unfortunately, no observation of the aerosol chemical composition was available during the campaign to evaluate those results. Vertical distributions indicate a rather constant structure of aerosol molar concentrations for F13, with a mean value around 6.2 nmol cm⁻³ for both SO₄ and NH₄ and a value of 0.5 nmol cm⁻³ for NO₃ (Fig. 4). For F21 and F29 simulated results show peak aerosol concentrations in the mid-troposphere up to a factor of 2 compared to F13 and a larger vertical gradient, with large and moderate depletion in the boundary layer for F21 and F29 (Fig. 4b and c). F21 and F29 have NH₄ mean values of 8 and 10.2 nmol cm⁻³, respectively, and SO₄ mean values both around 7 nmol cm⁻³. These values and the vertical structures correspond relatively well with mean observed concentrations for NH₄ and SO₄ of 7 nmol cm⁻³ seen during the ARCTAS (Arctic Research of the Composition of the Troposphere from Aircraft and Satellites) and ARCPAC (Aerosol, Radiation, and Cloud Processes affecting Arctic Climate) campaigns in April 2008 (Fisher et al., 2011). Fisher et al. (2011) showed that volcanic sources (Aleutian Islands and Kamchatka) accounted for 12 %–24 % of the sulfate at all altitudes, with a peak contribution in the mid-troposphere. The volcanic source is discharged directly in the free troposphere and is thus less affected by deposition than surface sources. This is also supported by satellite observations from the Ozone Monitoring Instrument (OMI) over the North Slope of Alaska, which shows much larger SO₂ concentrations at the end of ISDAC. Clouds sampled during both F21 and F29 appear to form mostly in air masses containing dust and smoke, possibly with a highly acidic coating.

Figure 4 Vertical profiles of sulfate (a), ammonium (b) and nitrate (c) molar aerosol concentration along flights F13 (green), F21 (red) and F29 (light blue).

Figure 5 presents the vertical profile of the neutralization fraction f_n (full line, see Eq. 18) and the contact angle θ (dashed, see Eqs. 20 and 21) for MYKE2 (Fig. 5a) and for MYKE4 (Fig. 5b) along the top of the three flights F13, F21 and F29. Results obtained with MYKE2 and MYKE4 using the same value of the neutralization fraction are very similar. Results from the two simulations are therefore discussed together. The difference lies in the curve shape of the contact angle θ: MYKE4 simulates a more rapid decrease for $\mathrm{0}<f_{\mathrm{n}}<\mathrm{0.5}$ than MYKE2 (Fig. 1). This prescription substantially increases θ values in MYKE4 more than in MYKE2 along the vertical profile by up to 3^∘, especially at the cloud top where nucleation is the dominant process. This change has a positive impact on the nucleation rate: a smaller contact angle in the MYKE2 simulation indeed tends to decrease the critical Gibbs free energy to form ice embryos (Eq. 15), hence leading to a higher nucleation rate of ice crystals. The θ profile in F13 presents a constant shape with values around 17.5 and 20.5^∘ for MYKE2 and MYKE4, respectively. Focusing on MYKE4 for F21, the large contact angle around 21^∘ corresponds to acid INPs, i.e., a smaller f_n than F13, and a decrease in the nucleation rate. Although F29 also shows significant acidity around 400 hPa (Fig. 4b), with higher concentrations of SO₄ than F13, it tends to neutrality around 500 hPa in relation to the increase in ammonium at this altitude in comparison to higher altitudes and the negligible amount of nitrate in the upper part of the cloud (Fig. 4b and c).

Figure 5 Vertical profiles of the neutralization fraction (f_n, full line) and the contact angle (θ, dashed line) for MYKE2 (a) and MYKE4 (b) along flights F13 (green), F21 (red) and F29 (light blue).

Our results reveal that the model broadly reproduces N_a from the ground to the 500 hPa level, but it has difficulty representing N_a in the upper part, even if observations and model results remain of the same order of magnitude. MYKE2 and MYKE4 simulations show higher θ values at cloud top for F21 and F29 in comparison to F13, thus differencing the acidic from the nonacidic cases as expected. In the following section, we will examine the effect of interactive chemistry on the cloud microphysical variables.

4.3 Cloud microphysical structure

Details of the retrieval of cloud microphysical properties and associated uncertainties from the several cloud probes onboard the Convair-580 aircraft are given in Jouan et al. (2012). Figure 6 presents the comparison of the observed and simulated (REF, MYKE2 and MYKE4) vertical profiles of IWC (uncertainties: ±75 %) along the three flights. Observed IWC vertical profiles for F13 and F29 continuously decrease between 800 and 400 hPa, with values in the range of 10⁻¹ to 10⁻² kg kg⁻¹. For flight F21, observed IWC shows a large variability in its vertical structure. IWC values simulated by both MYKE2 and MYKE4 are very similar, with a slight improvement for MYKE2 simulating more IWC. This agrees with the θ difference between MYKE2 and MYKE4 (Fig. 5). A smaller contact angle in the MYKE2 simulation tends to decrease the critical Gibbs free energy to form ice embryos (Eq. 15), hence leading to a higher nucleation rate of ice crystals and higher IWC. Both MYKE2 and MYKE4 broadly capture observed values with a low bias: $+\mathrm{1.2}\times {\mathrm{10}}^{-\mathrm{2}}$ and $+\mathrm{8.1}\times {\mathrm{10}}^{-\mathrm{3}}$ g kg⁻¹ for F13; $-\mathrm{3.2}\times {\mathrm{10}}^{-\mathrm{3}}$ and $-\mathrm{3.5}\times {\mathrm{10}}^{-\mathrm{3}}$ g kg⁻¹ for F21; $-\mathrm{2.1}\times {\mathrm{10}}^{-\mathrm{3}}$ and $-\mathrm{8.1}\times {\mathrm{10}}^{-\mathrm{3}}$ g kg⁻¹ for F29. In contrast, REF strongly underestimates IWC values with a negative bias of 0.01 g kg⁻¹ for F13 and 0.03 g kg⁻¹ for F29. Note that REF does not have any noticeable IWC cloud at these levels in flight F21.

Figure 6 Comparison of the observed (red) and simulated (REF in green, MYKE2 in purple and MYKE4 in cyan) vertical profiles of IWC along flights F13 (solid lines), F21 (dashed lines) and F29 (solid line with diamond markers).

Figure 7 presents a comparison between observed and simulated (REF, MYKE2 and MYKE4) vertical profiles of ice number concentration (N_i) (uncertainties: ±50 %) in the upper part of the cloud where heterogeneous ice nucleation processes are dominant above 500 hPa during flights F13, F21 and F29. The airborne ISDAC vertical profile for the TIC1 observed during F13 varies between 70 and 200 L⁻¹ and is rather constant with altitude. The REF simulation strongly underestimates N_i by 2 orders of magnitude, corresponding rather to a TIC2. MYKE2 and MYKE4 reproduce the observed N_i within the ranges of uncertainties well, while MYKE4 is slightly closer to observations with a bias of 25 L⁻¹. The TIC2 cloud type observed along F21 and F29 flight tracks is characterized by a small concentration of ice crystals ranging between 1 and 30 L⁻¹. For F21, while REF is not able to simulate a persistent cloud, both MYKE2 and MYKE4 show a cloud with N_i close to observations typical of TIC2 under 450 hPa in the range of uncertainties (±50 %). As expected, due to the biases of temperature and relative humidity over ice, the model underestimates the cloud-top altitude for F21. For F29, both MYKE2 and MYKE4 show an increase in N_i compared to REF, which has the best statistics, while MYKE2 and MYKE4 simulations are overestimated by 1 order of magnitude. However, it is reasonably close to satellite observations as analyzed by Keita et al. (2019). Their analysis revealed a large discrepancy in N_i between ISDAC flights and satellite estimations for F29 in the upper part of the cloud. We can notice here that the order of magnitude of N_i for F29 estimated from satellites can call into question the classification of F29 as a TIC2, especially as Jouan et al. (2012), using a flight track above Barrow (now Utqiagvik) instead of Fairbanks, classified this cloud as a TIC1. This discrepancy between airborne measurements, simulated results and satellite observations can be due to the small sampling domain during ISDAC versus the low resolution of satellite products and of the model grid.

Figure 7 Comparison of the observed (red) and simulated (REF in green, MYKE2 in purple and MYKE4 in cyan) vertical profiles of N_i along flights F13 (solid lines), F21 (dashed lines) and F29 (solid line with diamond markers).

Figure 8 presents the comparison of the observed and simulated (REF, MYKE2 and MYKE4) vertical profiles of the mean ice crystal radius (R_i) with uncertainties of ±97 %) along the F13, F21 and F29 flights. Observations show that, although having the same IWC magnitude (Fig. 6), the TIC1 and TIC2 differ in their N_i (Fig. 7) and R_i values. Flight F13 (TIC1), with a large N_i concentration, has R_i values around 25 µm, while both F21 and F29 refer to TIC2 with a low N_i and R_i at least a factor of 2 larger. The INP acid coating in TIC2 inhibits the ice nuclei properties of the INPs, slowing the rate of ice nucleation in comparison to uncoated N_i. Subsequently, this decrease in the nucleation rate increases the amount of available supersaturated water vapor and allows the rapid growth of activated ice crystals. It could explain the persistence of low N_i and the large R_i. For flight F13, MYKE2 and MYKE4 simulate the TIC1 formation above 450 hPa relatively well in the observation range, while below 450 hPa they both overestimate R_i by a factor of 2. For this TIC1 cloud, MYKE2 and MYKE4 give the smallest error in comparison to REF. For flight F21, MYKE2 and MYKE4 improve the comparison of simulated R_i against observations, showing large ice crystals even if the cloud-top altitude is underestimated. For flight F29, observed values of R_i are even larger. MYKE2 and MYKE4 show a little improvement in comparison to REF but only above around 450 hPa, with larger simulated ice crystals than REF. For flights F21 and F29, MYKE2 and MYKE4 underestimate the observed R_i by a factor of 2.

Figure 8 Comparison of the observed (red) and simulated (REF in green, MYKE2 in purple and MYKE4 in cyan) R_i along flights F13 (solid lines), F21 (dashed lines) and F29 (solid line with diamond markers).

4.4 Discussion

Our analysis shows the poor performance of the original REF parameterization in representing ice heterogeneous nucleation with low IWC and reveals that the MYKE parameterization can significantly improve the representation of the IWC at all vertical levels in polluted or unpolluted air masses. Along the three flights, RH_i is therefore lower in the MYKE2 and MYKE4 simulations than in the REF run at cloud top. This may be due to the new parameterization promoting ice nucleation through a reduction of the available supersaturated water vapor. The new parameterization, with the variation in time and space of A_d and N_t, better represents N_i and R_i values at the top of TICs for flights F13 and F21 where the nucleation occurs. The pronounced slope of observed R_i above the 500 hPa level in TIC2 cases (Fig. 8) indicates rapid growth of the ice crystals, which consume supersaturated water vapor faster than it is made available in the model. Finally, for flight F29, the new parameterization slightly improves R_i at the top of the clouds, while under around the 450 hPa level, simulated results show better agreement for the REF simulation. The reason for that is not clear. However, Fig. 5 shows a decrease in θ with altitude between 450 and 500 hPa in connection with an increase in the ammonium molar concentration (Fig. 5b), which leads to a more efficient heterogeneous nucleation of ice at this altitude with smaller ice crystals and larger concentrations.

Finally, from the comparison of the three simulations, we can assess the ability of the new scheme to discriminate TIC1 and TIC2 clouds. For F13, while REF results in a TIC2 cloud, MYKE2 and MYKE4 simulations produce a TIC1 in agreement with observations. As shown before, the orders of magnitude of N_a at the top of the cloud for F13 and F29 are similar, but the neutralization fraction f_n shows more acidic aerosols for F29. For both cases, close values of IWC allow us to compare MYKE results for N_i and R_i. Looking at the top of the cloud (above the 440 hPa level), N_i is lower for F29 than for F13 and R_i is larger for F29 than for F13, responding to acid aerosol through the variation of the contact angle. Within the limit of our calculation, the new parameterization significantly improves the representation of nucleation in TIC1 for F13 versus TIC2 for F29 at the cloud tops, despite the model bias of simulated aerosols by WRF-Chem over the Arctic (Mölders et al., 2011). The comparison between simulations of F21 and F13 cases with MYKE is not so clear. Even if, at the top of the cloud, N_i is lower for F21 than for F13 as expected, R_i is smaller for F21 than for F13, which is not consistent with TICs. However, the comparison of the f_n fraction at the cloud tops shows similar values for F21 and F13 near acid neutrality. This result highlights the importance of a consistent simulation of aerosol physicochemical properties to create a valuable simulation of microphysical ice cloud properties with our new parameterization of heterogeneous ice nucleation.

In general, regarding overall simulated results, MYKE4 shows better agreement with observations than MYKE2 for TIC1 and TIC2 clouds. It is well known that the effect of acid coating on INPs is to reduce their ability to form ice crystals, and this effect increases with the amount of acid (Sullivan et al., 2010; Yang et al., 2011). Moreover, our results suggest that even a low acidity on INPs leads to an important decrease in the heterogeneous ice nucleation rate because, for MYKE4, θ increases more rapidly when acid coating increases, i.e., decrease in the f_n fraction (Fig. 1).

5 Conclusions

A new parameterization of ice heterogeneous nucleation for water-subsaturated conditions, based upon the CNT approach and coupled with real-time chemistry information, is proposed within the WRF-Chem model. The coupling with chemistry helps to link the contact angle θ to the aerosol neutralization fraction, which is a good proxy for the acidity of aerosols. This new parameterization is implemented in the Milbrandt and Yau (2005a, b) two-moment cloud microphysical scheme available in WRF-Chem. It is particularly designed to simulate Arctic ice clouds. In the Arctic, ice clouds are separated into two classes: (1) TIC1 clouds characterized by large concentrations of very small crystals and TIC2 clouds characterized by low concentrations of larger ice crystals. TIC2 clouds induce significant ice crystal precipitation or so-called diamond dust, a notoriously deficient variable to simulate in polar atmospheric models despite its significant contribution to annual snowfall that is generally reported as “trace” in station observations. The model including the original Milbrandt and Yau (2005a, b) scheme and the modified one are applied to three test cases observed during ISDAC: one TIC1 and two TIC2 clouds. For each case, results are analyzed in terms of meteorology, chemistry and cloud microphysical properties by comparison between the new (MYKE2 and MYKE4) and original (REF) parameterization of ice nucleation within the cloud microphysical scheme and with available observations.

Our results show the poor performance of the REF parameterization in representing Arctic ice cloud types at low IWC and underline the fact that the MYKE2 and MYKE4 parameterizations significantly improve the representation of IWC, especially in the top region of the clouds where nucleation dominates, in both polluted and unpolluted air masses. MYKE2 and MYKE4 simulations are in better agreement with observations for the three flights. In contrast, REF always strongly underestimates IWC values with a negative bias and does not see any noticeable IWC cloud at these levels during flight F21.

Aerosol number concentrations are simulated with the same order of magnitude as observations under the 550 hPa level, whereas above the 550 hPa level, the simulated value is overestimated by a factor of 3 for flight F13 and underestimated by 1 order magnitude for flight F29. Despite known difficulties in simulating aerosol concentrations in WRF-Chem over the Arctic region (Mölders et al., 2011), our parameterization achieves the proper representation of cloud types TIC1 for flight F13 versus TIC2 for flights F21 and F29 in the nucleation region at cloud top. Values and vertical structures of ammonium and sulfate molar aerosol concentrations for flights F21 and F29 correspond fairly well to mean observed concentrations, i.e., 7 and 5.5 nmol cm⁻³, during the ARCTAS and ARCPAC campaigns, respectively, with known contributions from volcanic sources peaking in the mid-troposphere. MYKE2 and MYKE4 simulations are similar, showing higher θ values at cloud top for flights F21 and F29 in comparison to flight F13, thus differencing the acidic from the nonacidic cases as expected, and a low sensitivity to the arbitrarily parameterized curve shape.

For the TIC1 case, REF strongly underestimates the ice crystal number concentration by at least 2 orders of magnitude and overestimates the mean radius, resulting in the false representation of an ice cloud corresponding rather to a TIC2. In contrast, the new parameterization captures the cloud type well, with representative microphysical structure (IWC, ice crystal mean radius and ice crystal number concentration) at the top of the cloud where nucleation occurs. TIC2 clouds observed along the F21 and F29 flight tracks are characterized by a small concentration of ice crystals ranging between 1 and 30 L⁻¹. MYKE2 and MYKE4 simulate those ice crystal number concentrations within the range of observation uncertainties. For flight F21, REF is not able to simulate a persistent cloud, while both the MYKE2 and MYKE4 simulations show a cloud with an ice crystal concentration close to observations. Corresponding values are typical of TIC2 cloud under the 450 hPa level even if the model underestimates the cloud-top altitude as a result of biases in the simulated temperature and relative humidity over ice. MYKE2 and MYKE4 also improve the ice crystal mean radius, showing larger ice crystals than REF. For flight F29, both MYKE2 and MYKE4 show an increase in the ice crystal concentration compared to REF, which has the best statistics, but the MYKE2 and MYKE4 results are still overestimated by 1 order of magnitude. MYKE2 and MYKE4 slightly improve the representation of the ice crystal mean radius in comparison to REF above the 450 hPa level, with larger simulated ice crystals than REF. For both TIC2 flights, MYKE2 and MYKE4 nevertheless underestimate the observed mean radius by a factor of 2. Comparing the two versions of the parameterization for the three cases, in general, MYKE4 presents a slight improvement compared to MYKE2, in agreement with θ dependency. Because this difference is small, the dependency of the contact angle on the aerosol neutralization fraction under a concave form should be considered a sufficient condition to improve the representation of heterogeneous ice nucleation in Arctic ice clouds.

In our simulations, secondary organic aerosol (SOA) formation is not considered. However, the concentration of their precursor species, mainly biogenic and aromatic volatile organic compounds, should be low in the ISDAC region and period as suggested by the WRF-Chem simulation. However, results obtained later during the NETCARE campaign (2015) show a potential contribution of SOA to the total mass of Arctic aerosols, but their precursors are not yet identified in the Arctic, which is a new challenge in simulating their formation (Abbatt et al., 2019). Moreover, as our parameterization is dedicated to the simulation of Arctic ice cloud types, we are confident that the combination of CBM-Z and MOSAIC is appropriate even if CBM-Z is a relatively simple gas-phase mechanism and if SOA formation is not considered. Indeed, our results suggest that it is enough to consider the chemical impact on heterogeneous ice nucleation though the degree of aerosol acidity acting as INPs. Despite the huge challenge, our parameterization seems promising. Further studies will help with validations against satellite data and future campaigns. In particular, future flight campaigns should include simultaneously measurements of cloud microphysics properties, aerosol number size distribution, aerosol chemical composition and ice nuclei number concentrations. The next step will be to extend simulations to quantify the role of the ice nucleation of acid pollution in radiation, the atmospheric water balance and, ultimately, the Arctic climate.