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 es-sential 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 het-erogenous ice nucleation is most likely occurring, even with the known bias of simulated aerosols by WRF-Chem over the Arctic.

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;Mc-Farquhar 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 timedependent 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., 2010DeMott et al., , 2015Cirisan 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 −1 . 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 −1 . 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 nu-cleation 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 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 ∈ (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: and 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 where c x = ρ x π 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 −3 . The size spectrum of each category is described by a common generalized gamma distribution function (Schoenberg Ferrier, 1994) of the form 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 where N 0x is the intercept parameter given by 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 −1 ) and the mass mixing ratio q T,x (kg kg −1 ) are respectively given by and where m x (D) is obtained from Eq. (3). 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 ∼ −32 • C. In the range −30 to −50 • C, MY05 follows DeMott et al. (1994) with 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, 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 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): 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).

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 where A d is the total surface area of dust particles and N t is the total number concentration of available INPs. In pre- 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.

Hydrometeor Source terms
Cloud nucleation, condensation and evaporation, self-collection Rain autoconversion, evaporation, accretion of cloud, self-collection, melting of frozen hydrometeors Ice nucleation (contact, deposition, condensation freezing, rime splintering, immersion, homogeneous freezing of cloud), riming of cloud, deposition and sublimation Snow conversion from ice (including ice aggregation), collection of ice and cloud, deposition and sublimation, aggregation (self-collection), collisional freezing with rain Graupel collisional freezing of rain and ice-snow-graupel, conversions from ice and snow, collection of cloud and ice, deposition and sublimation Hail collisional freezing of rain and ice-snow-graupel, collection of cloud-rain-ice-snow, deposition and sublimation, probabilistic freezing of rain, conversion from graupel vious 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 where B = 10 26 cm −2 s −1 is the kinetic coefficient (Pruppacher and Klett, 1997), k is the Boltzmann constant (J K −1 ), 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 where σ iv = 106.5 × 10 −3 J m −2 is the surface tension between ice and water vapor, ρ i = 0.5 g cm −3 is the bulk ice density and R v = 461.5 J kg −1 K −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: where φ = 1 − 2q cos θ + q 2 and q = r n r g , with r g being the critical germ size expressed as 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. (2008Eastwood et al. ( , 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 2 ), 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 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 immersionfreezing 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 icenucleating 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. (2008Eastwood et al. ( , 2009) from heterogeneous nucleation rates on kaolinite particles obtained in laboratory measurements. As a best fit, they found limiting values of θ = 26 • in polluted air and θ = 12 • 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. (2018Kumar et al. ( , 2019a 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: with p larger than 1. We have chosen a quadratic (p = 2) form for simplicity: 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: 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 0 < f n < 0.5 with a more (Eq. 21) or less (Eq. 20) rapid decrease between 0.5 Figure 1. Variation of f n with θ for MYKE2 (blue line) and MYKE4 (green line). 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.   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 Monin-Obukhov Janjic Eta scheme (Janjić, 1994) Land surface Unified Noah land-surface model (Chen and Dudhia, 2001) Chemistry  (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).

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 2 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.   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 −3 on the whole column with mean concentrations around 73 cm −3 , very close to the simulation mean of 86 cm −3 . 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 −3 ), partic- ularly at altitudes above 550 hPa near cloud top where peak concentrations exceeding 1000 cm −3 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 4 presents simulated (REF, MYKE2 and MYKE4) vertical profiles of sulfate (SO 4 ), ammonium (NH 4 ) and nitrate (NO 3 ) 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 in-dicate a rather constant structure of aerosol molar concentrations for F13, with a mean value around 6.2 nmol cm −3 for both SO 4 and NH 4 and a value of 0.5 nmol cm −3 for NO 3 (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 4 mean values of 8 and 10.2 nmol cm −3 , respectively, and SO 4 mean values both around 7 nmol cm −3 . These values and the vertical structures correspond relatively well with mean observed concentrations for NH 4 and SO 4 of 7 nmol cm −3 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 2 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 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 0 < f n < 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 4 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).

Aerosol properties
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.

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 −1 to 10 −2 kg kg −1 . 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: +1.2 × 10 −2 and +8.1 × 10 −3 g kg −1 for F13; −3.2 × 10 −3 and −3.5 × 10 −3 g kg −1 for F21; −2.1 × 10 −3 and −8.1 × 10 −3 g kg −1 for F29. In contrast, REF strongly underestimates IWC values with a negative bias of 0.01 g kg −1 for F13 and 0.03 g kg −1 for F29. Note that REF does not have any noticeable IWC cloud at these levels in flight F21. 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 −1 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 −1 . 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 −1 . 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 hu- midity 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 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.

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).

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 −3 , 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 −1 . 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 MO-SAIC is appropriate even if CBM-Z is a relatively simple gas-phase mechanism and if SOA formation is not consid-ered. 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.
Author contributions. SAK and EG developed and implemented the parameterization with the support of JCR and ML. SAK performed the simulations with the technical support of TO. SAK analyzed results and wrote the paper with the support of JCR, ML and JPB. All authors contributed to the paper and to the analysis.
Competing interests. The authors declare that they have no conflict of interest.
Acknowledgements. As we know that Eric Girard had numerous discussions with our colleague Allan Bertram on the form of the dependency between the contact angle and the aerosol neutralization fraction, we would like to thank him for his invaluable help. We thank NETCARE (Network on Climate and Aerosols: Addressing Key Uncertainties in Remote Canadian Environments) and NSERC (Natural Sciences and Engineering Research Council of Canada) for funding support and ARM (Atmospheric Radiation Measurement Program) for the data collected during ISDAC. Computer modeling benefited from access to IDRIS HPC resources (GENCI allocations A003017141 and A005017141) and the IPSL mesoscale computing center (CICLAD: Calcul Intensif pour le CLimat, l'Atmosphère et la Dynamique). We acknowledge the use of the WRF-Chem preprocessor tools mozbc and fire_emiss provided by the Atmospheric Chemistry Observations and Modeling Lab (ACOM) of NCAR.
Review statement. This paper was edited by Samuel Remy and reviewed by three anonymous referees.