Simulating Arctic Ice Clouds during Spring Using an Advanced Ice Cloud Microphysics in the WRF Model

Abstract

Two Types of Ice Clouds (TICs) have been characterized in the Arctic during the polar night and early spring. TIC-1 are composed by non-precipitating small ice crystals of less than 30 µm in diameter. The second type, TIC-2, are characterized by a low concentration of large precipitating ice crystals (>30 µm). Here, we evaluate the Weather Research and Forecasting (WRF) model performance both in space and time after implementing a parameterization based on a stochastic approach dedicated to the simulation of ice clouds in the Arctic. Well documented reference cases provided us in situ data from the spring of 2008 Indirect and Semi-Direct Aerosol Campaign (ISDAC) campaign over Alaska. Simulations of the microphysical properties of the TIC-2 clouds on 15 and 25 April 2008 (polluted or acidic cases) and TIC-1 clouds on non-polluted cases are compared to DARDAR (raDAR/liDAR) satellite products. Our results show that the stochastic approach based on the classical nucleation theory, with the appropriate contact angle, is better than the original scheme in WRF model to represent TIC-1 and TIC-2 properties (ice crystal concentration and size) in response to the IN acidification.

1. Introduction

Over recent decades, temperatures in the Arctic have increased at more than twice the global rate, largely as a result of growing concentrations of greenhouse gases and complex feedback processes [1]. One of the main factors influencing the Arctic rapid warming is how clouds interact with aerosols and radiation [2]. Clouds play a fundamental role in regulating the radiative balance by influencing both solar radiation and infrared radiation both at the surface and at the top of the atmosphere. One of the main features of the Arctic radiative environment is the absence of solar radiation for a significant portion of the year. Low temperatures and high relative humidity during cold seasons make the Arctic a very suitable environment for the formation of ice clouds [3], both in the upper and lower troposphere [4,5]. Shupe [5] notes that ice-only clouds are more common than mixed-phase clouds (MPCs) with similar longevity based on mean-annual statistics over three observation sites (Barrow, Eureka, SHEBA). However, the existence of MPCs is important in these regions during winter, when surface warming resulting from cloud infrared radiation is not compensated by cooling from the increased cloud shortwave albedo via Twomey effect [6,7].

Representation of cloud microphysical processes in climate models is challenging because fundamental microphysical processes are often poorly represented as compared to more complex limited-area models [8]. This is particularly true when considering a higher level of complexity through aerosol-cloud interactions [9,10]. For example, despite the evolution of our knowledge on the processes of nucleation, growth and dissipation of ice crystals [11], the performance of respective parameterizations in models might differ due to the modelling scale and scientific objective. Parameterizations designed for lower latitudes are nevertheless often used on the global scale without being properly tested at high latitudes. The “parameterization way of thinking” turns out to be a relevant method to improve models using observations [12]. The formation of ice crystals can be triggered through homogeneous freezing at very cold temperatures, but also through different heterogeneous nucleation processes upon aerosol particles acting as ice nuclei (IN) [13]. Over the past 10 years, several parameterizations of homogeneous and heterogeneous ice nucleation have been developed to reflect the various physical and chemical properties of aerosols in addition to temperature and/or ice supersaturation and distribution functions [10,14]. Keita et al. [15] tested different parameterizations for the ice nucleation based on stochastic and deterministic approaches [16] in the limited-area version of the Global Multiscale Environmental Model (GEM-LAM) [17]. After comparing simulations of dynamics and microphysics variables against in situ observations from the Indirect and Semi-Direct Aerosol Campaign (ISDAC) [18], they found that the stochastic scheme was the most suitable parameterization for simulating ice clouds in the Arctic. But the validation of the model performance, both in space and time, had not been fully performed yet. Heterogeneous nucleation requires less supersaturation and is favored in low level cloud systems. Acidification of particles has however been raised as potentially playing a significant inhibition role [19], and evidence of such effect has been investigated from airborne and space observations [20].

During ISDAC, research flights were coordinated with CloudSat and CALIPSO (Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation) satellite overpasses. This colocation of satellite and in situ field data provides an efficient way to study Arctic clouds processes, to develop and to evaluate new parameterizations. Satellite instruments have the advantage of providing data over much larger spatial scales and the cloud vertical structure. Ten years satellite observations from the A-train, and particularly with the spaceborne CloudSat radar and CALIOP (Cloud-Aerosol Lidar with Orthogonal Polarization) lidar onboard the CALIPSO platform have deeply improved Arctic cloud-climate studies. Since polar clouds are geometrically and optically thinner than mid-latitude clouds, those active sensors return signals of the cloud vertical structure that are less attenuated in the vicinity of the Arctic surface [12]. As a result, the vertical distribution of aerosol, cloud and hydrometeors structure is observed by the tandem CloudSat and CALIPSO in polar regions with much deeper penetration through the atmosphere than at lower latitudes [21,22]. Figure 1 shows the study location where ISDAC campaign took place and several satellite orbits that will be relevant to the current study. In particular, the Lower panel of Figure 1 clearly illustrates the formation of thin ice clouds (semitransparent to the lidar) where ice nucleation occurs directly on volcanic plumes. These plumes where traced back [20,23] to several active volcanoes in Kamchatka and Aleutian Islands during the study period. Although, no composition measurements were available from ISDAC, it is well known [24], that such volcanic emission are typically composed largely of sulfuric and nitric acids aerosols.

The main objective of this study is to improve the ice nucleation parameterization in the Weather Research and Forecasting model (WRF) model using one of the best heterogeneous ice nucleation parameterizations developed by Keita et al. [15]. From previous studies [25,26], two types of ice clouds (TIC) have been distinguished. The first type (TIC-1) is composed by a large concentration of small ice crystals whereas the second type (TIC-2) is composed of a smaller concentration of larger and precipitating ice crystals. As a reference, we simulate cloud distributions and structures corresponding to four well documented [20] cloud cases obtained during ISDAC. The observed clouds are characterized by comparable levels of Ice Water Content (IWC), but different air mass origins. In particular, the occurrence of various acid aerosols concentrations in those air masses impacted the ice nuclei properties and the cloud microstructure [20,23]. Based on these local studies, we attempt to generalize the effect of IN acidification on polar cold clouds and then verify this extrapolation against a large ensemble of observations. Thus, our simulations are evaluated through a comparison with both ISDAC in situ data and DARDAR (raDAR/liDAR) extended satellite products. Our results represent an effort to link regional modeling, remote sensing and laboratory studies mentioned above with in situ observations of ice clouds. In Section 2, the model experiments and ice nucleation parameterizations are described. Section 3 presents the observational dataset. Results are presented and discussed in Section 4, and conclusions are drawn in Section 5.

2. WRF

2.1. Model Description and Set-up

The Advanced Research WRF (Weather Research and Forecasting) model Version 3.5.1 Version 3.5.1 is used in this study to simulate the evolution of clouds during ISDAC campaign, from 1 to 30 April 2008, over the domain shown in Figure 1 [27]. The domain is based on a Lambert projection centered on Barrow, Alaska over 160 × 100 grid cells with a horizontal resolution of 10 km and 55 vertical levels between the surface and 50 hPa. Three simulations have been carried out in this study (Section 2.3). The first 4 days of the simulations (1–4 April included) are used for model spin-up. WRF options and parameterizations used in this study are summarized in

Table 1. Parameterizations and options used for the WRF simulations.

Meteorological Option	Selected Option
Microphysics	Milbrandt and Yau [29]
Shortwave radiation	Rapid Radiative Transfer Model for Global climate models (RRTMG) [30]
Longwave radiation	RRTMG [30]
Cumulus parameterization	Kain-Fritsch-CuP (KF-CuP) [31]
Planetary boundary layer	Mellor-Yamada-Janjic [32]
Surface layer	Monin-Obukhov Janjic Eta scheme [32]
Land surface	Unified Noah land-surface model [33]

. Meteorological initial and boundary conditions of the atmospheric fields are provided by NCEP (National Centers for Environmental Prediction) Global Forecast System (GFS) Final Analysis (FNL) data at 1° resolution, with 26 pressure levels and updated every 6 h. The surface conditions are initialized with NCEP-archived 0.5° sea-surface temperatures fields (updated every 6 h) and sea ice data (updated once a day). The simulation is nudged to FNL winds, temperature, and humidity updated every 6 h [28] above the planetary boundary layer (PBL).

2.2. Modification of the Microphysical Scheme

In this study, we use the two-moment version of Milbrandt and Yau [29] cloud microphysical scheme (Table 1). It includes the following prognostic variables: the mixing ratio and the number concentration of cloud liquid water, cloud ice water, rain, snow, hail and graupel. In the model, all ice particles are assumed to be just irregular particles under the form of bullet rosettes. In this section, the emphasis is put on the modification of the parameterization of ice crystal nucleation, central in our study.

In the atmosphere ice crystals form through heterogeneous and homogeneous ice nucleation. Homogeneous freezing is the spontaneous freezing of a water (or haze) droplets. According to Pruppacher and Klett [13], the homogeneous freezing rate of cloud droplets is significant at temperatures below ~−38 °C. The parameterization for homogeneous freezing of water droplets follows DeMott et al. [34] in the range −30 °C to −50 °C. Heterogeneous nucleation involves solid substrates called IN initiating the formation of ice crystals. Several mechanisms are known by which aerosol particles promote the formation of the ice phase in clouds: deposition nucleation, immersion freezing, condensation-freezing and contact freezing. In the original version of the Milbrandt and Yau [29] microphysical scheme, deposition and condensation-freezing are functions of water vapor supersaturation with respect to ice following Meyers et al. [35]. 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. Contact freezing is parameterized following Young [36] where the number concentration of contact IN is a function of temperature following Meyers et al. [35]. In the contact formation mode, ice nucleation occurs on a solid particle reaching contact with a water droplet. Immersion freezing of activated rain and cloud water droplets follows the parameterization of Bigg [37]. In the immersion mode, ice nucleation occurs on a solid particle immersed in either an aqueous solution or a haze particle in sub-saturated air with respect to liquid water or in an activated cloud water droplet.

Several modifications have been made to the original version of the microphysical scheme to account for the heterogeneous ice nucleation on both uncoated and sulfuric acid-coated IN. Deposition ice nucleation on uncoated IN and immersion freezing of sulfuric acid-coated IN are parameterized following Girard et al. [19]. These parameterizations are based on the classical nucleation theory (CNT) [13]. The CNT is a stochastic approach in which the nucleation rate depends on the contact angle θ between the ice embryo and the IN. In the parameterization used in Girard et al. [19], the contact angle is derived from the heterogeneous nucleation rate, which is in turn derived from ice onset measurements of Eastwood et al. [38,39]. Their laboratory studies used uncoated (θ = 12°) and sulfuric acid-coated (θ = 26°) kaolinite particles to derive a parameterization based on a contact angle between ice and a substrate. Using a parameterization of ice nucleation on kaolinite particles is supported by the favorable meteorological situation encountered during the ISDAC period. Field observations have highlighted the presence of a large occurrence of haze, smoke and dust particles during ISDAC period [23,40,41,42]. Besides, in the Arctic region, a large fraction of the aerosol particles (including insoluble aerosols such as mineral dust) can be coated with acidic sulfate [43]. Eastwood et al. [39] results support the idea that anthropogenic emissions of SO2 and NH3 may influence the ice nucleating properties of mineral dust particles by increasing the relative humidity required for ice nucleation. This shift in ice nucleation conditions may influence the frequency and properties of ice clouds.

In a given time step (Δt) the number concentration of nucleated ice crystals (N_f) is given by:

$$N_f\left(\Delta t\right)=N_0\left[1-exp\left(-J{A}_{kaolinite}\Delta t\right)\right]$$

(1)

where $A_{kaolinite}$ is the surface area of the kaolinite particles, $N_0=$ ${10}^4{\mathrm{m}}^{-3}$ is the total number concentration of available IN and J is the nucleation rate of embryos per unit surface of particles and is defined as:

$$J({\mathrm{cm}}^{-2}{\mathrm{s}}^{-1})=Bexp\left(\frac{-{\Delta \mathrm{G}}^{*}}{kT}\right)$$

(2)

where B = 1.521 ${10}^{37}{\mathrm{cm}}^{-2}{\mathrm{s}}^{-1}$ is the pre-exponential factor [13], k the Boltzmann constant in ${\mathrm{JK}}^{-1}$, $T$ is the temperature in kelvin (K), ${\Delta \mathrm{G}}^{*}$ the critical Gibbs free energy for the formation of an ice embryo in joules (J), defined as:

$$\mathrm{where} \Delta G^{*}=\frac{16\pi {\sigma }_{iV}^{3}f\left(\cos \mathsf{\theta }\right)}{3{\rho }_i^2{R}_v^2{T}^{2}l{n}^2{S}_i}$$

(3)

${\sigma }_{iv}$ = 106.5 ${10}^3$ J·${\mathrm{m}}^{-2}$ the surface tension between ice and water vapor, ${\rho }_i$ = 0.9 g·${\mathrm{cm}}^{-3}$ the bulk ice density, R_v = 461.5 J·k${\mathrm{g}}^{-1}{\mathrm{K}}^{-1}$ the gas constant for water vapor, T the temperature and S_i the saturation ratio with respect to ice. The $f\left(cos\theta \right)$ is a monotonic decreasing function of the cosine of the contact angle θ as defined by Pruppacher and Klett [13] for an infinite plane surface.

2.3. Simulation Description

Three versions of the microphysical scheme are used for our simulations and detailed in

Table 2. Ice nucleation parameterization schemes used for the simulations. The original microphysical scheme (WRF) is based on the singular hypothesis, while WRF_np and WRF_p are based on the stochastic theory, which includes a contact angle (θ = 12° and θ = 26°, respectively).

Scheme Version	Author	Species	Ice Nuclei State	Ice Nuclei Contact Angle
WRF_np	[38]	Kaolinite	Uncoated	θ = 12°
WRF_p	[39]	Kaolinite	Coated	θ = 26°
WRF	[35]	---------	---------	---------

. WRF refers to the original version of the Milbrandt and Yau [29] microphysical scheme without any modification. Hence, it does not refer to any particular type of air masses. WRF_p and WRF_np, respectively “polluted” (acidic) and “not polluted” (clean) air masses, are the modified versions of the scheme based on the CNT and the parameterizations described in Section 2.2. and are associated to the deposition and immersion freezing nucleation modes. WRF_np and WRF_p use a constant contact angle of θ = 12° and θ = 26° representative of uncoated and sulfuric acid-coated kaolinite particles, respectively.

3. Observation Data

We present in this section the different observation datasets used in this study: the ISDAC campaign dataset (Section 3.1) and the DARDAR (raDAR/liDAR) satellite products (Section 3.2).

3.1. ISDAC Campaign

Arctic clouds simulated by WRF are evaluated against observations from the U.S. Department of Energy (DOE) Atmospheric Radiation Measurement (ARM) ISDAC campaign [18], which was conducted at North Slope of Alaska during April 2008. ISDAC primary aim was to examine the influence of the Arctic aerosol on cloud microphysical properties and on the surface energy balance. During the field campaign, the National Research Council of Canada (NRC) Convair-580 was equipped with 41 instruments including a complete set of dynamic, thermodynamic, radiation, aerosol and microphysical sensors. A set of 27 scientific flights was performed between April 1 and April 30, 2008 under various weather conditions. Flight profiles consisted of spiral, ramped ascents and descents through clouds, and constant altitude legs within and outside clouds between Barrow (72°N, 157°W) and Fairbanks (65°N, 148°W). The strategy was to coordinate flights with multiple aircrafts (NASA DC8, P-3B, B200, NOAA WP-3D) and satellite overpasses (A-Train).

The Convair-580 flew 27 sorties and collected data using an unprecedented 41 state-of-the-art cloud and aerosol instruments for a total of more than 100 h of flight time on 12 different days. Among the ISDAC in situ measurements, the following variables are used in this study: cloud ice content (IWC), number concentration of ice crystals (Ni), temperature (T) and relative humidity with respect to ice (RH_i=S_i x 100). The mean ice crystal radius (Ri) is calculated from IWC and Ni assuming spherical shape of ice and monodisperse size distributions. The instruments used to derive the variables and the way the bulk data were derived are given in Jouan et al. [23]. They have performed a comprehensive study of the ISDAC measurements derived by a variety of instruments and discussed the uncertainties in measurements/observations. They found that the uncertainties on the measurements are the following: ±11% for RH_i, ±50% for Ni, ±75% for IWC and ±0.5 °C for T. Using the IWC and Ni uncertainties, it can be shown that the uncertainty on Ri is ±97%. Details on uncertainties and biases related to these measurements are described in detail in McFarquhar et al. [18], Jouan et al. [23], and Keita et al. [15].

Table 3. Flight number, location, dates, time, Ni, Ri and type of air mass of the ice clouds investigated in this study.

Flight Number Location	Date (In 2008)	Time Profile	Ni/Ri	Air Mass
F12 Barrow	5 April	01:15:35 –01:34:05	(~3000 L−1)/(10–30 um)	Pristine
F13 Barrow	5 April	20:35:26 – 21:00:05	(100–200 L−1)/(25–50 um)	Pristine
F21 Barrow	15 April	00:55:40 – 01:17:24	(<30 L−1)/(50–140 um)	Polluted
F29 Fairbanks	29 April	04:08:22 – 04:27:51	(5–30 L−1)/(80–180 um)	Polluted

presents a brief description of the four specific ISDAC flights that have been used to validate our ice nucleation parameterization. Those flights have been selected according to the approach described in Keita et al. [15]. We have defined an observation box between (65°N, 142°W) and (72°N, 162°W) to identify the area where all flights have been performed. This flight area is schematically reported in Figure 1.

Clouds observed on 5 April 2008 during the F12 (01:15:35 – 01:34:06 UTC) and the F13 (20:35:26 – 21:00:06 UTC) flights are characterized by a relatively large crystal concentration Ni and small ice crystal radius Ri throughout the cloud. These clouds are of TIC-1 clouds according to the criteria of Jouan et al. [23], formed in a pristine environment. During our study, we analyzed measurements from the Passive Cavity Aerosol Spectrometer Probe (PCASP) instrument developed by Particle Measuring Systems (PMS). The PCASP is an optical particle counter used for measuring the concentration and size distribution of particles, assuming a spherical shape and refractive index between 1.52 and 1.59 in the shortwave domain. These particles are in the size range 0.12–3 μm. Figure 2 presents the aerosol concentrations observed in clear air just before entering the cloud regions. It shows that these air masses were relatively clean with a weak vertical variability of aerosol concentrations, staying mostly below 200 cm⁻³ from 800 to 400 hPa [41]. Average concentrations are around 100 cm⁻³.

In the second half of the campaign, ice clouds observed on 15 April 2008 during F21 flight (00:55:40–01:42:07 UTC) and on 25 April 2008 during F29 flight (04:08:22 – 04:27:51 UTC) are characterized by a small ice crystal concentration Ni, corresponding to the definition of TIC-2 clouds according to Jouan et al. [23]. It should be noted that there were no PCASP measurements during F21 flight. Aerosol concentrations measured during F29 flight are also reported in Figure 2. The PCASP instrument showed that there was a much higher concentration of aerosol particles in the lower troposphere (more than twice those observed during F12 and F13, e.g., larger than 400 cm⁻³) and more particularly at pressures below 550 hPa, near cloud top where peak concentrations exceeding 1000 cm⁻³ have been measured. This increase in aerosol particles is 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 the ISDAC campaign. These large deviations have been traced to volcanic activities in Kamchatka and Aleutian regions [20,23] that are known to be acidic. Clouds sampled during both F21 and F29 appear to form mostly in air masses containing dust and smoke, possibly with a highly acidic coating. Jouan et al. [20] found that the large SO₂ concentrations over the North Slope of Alaska during this period were due to volcanic emissions in the Aleutians. Furthermore, IN measurements on 15 April in Barrow show a very low IN concentration, compatible with the hypothesis of acidic inhibition during this period [20,23,41].

3.2. DARDAR Satellite Products

The DARDAR dataset is derived from a combination of the CloudSat, CALIPSO and MODIS (MODerate resolution Imaging Spectroradiometer) observations [44,45]). Cloud parameters such as vertical profiles of thermodynamic cloud phase, IWC, extinction and effective radius (Re) are obtained from a variational synergistic algorithm. In regions of the cloud detected by just one of these two instruments, their algorithm retrieves parameters from radar or lidar measurements alone [45]. The ice particle mass is obtained assuming the mass-diameter relationship proposed by Brown and Francis [46], and the effective radius Re is derived using the results of Foot [47]. The RMSE (root mean square error) in IWC and Re derived from DARDAR are retrieved according to the methodology of Delanoë and Hogan [44]. Jouan et al. [23] found that within TICs, the normalized RMSE of IWC and Re from DARDAR product ranged from 10 to 30% for IWC and less than 20% for Re. Retrieved properties are reported at the CloudSat horizontal resolution (about 1.4 km) and CALIPSO upper tropospheric vertical resolution (60 m). The lidar is highly sensitive to aerosol, haze and thin clouds like sub-visible cirrus. The backscattered signal can be observed down to cloud base for thin clouds with an optical depth less than about 3, where the signal is not completely attenuated. This apply to most ice clouds in the high latitudes, especially during the cold season. When the optical depth exceeds 3, the radar is still capable of penetrating the cloud down to its base. Delanoë et al. [23] have demonstrated the performance of DARDAR to retrieve cloud properties in cirriform layers and precipitating areas showing simultaneously thin and very thick ice clouds to highlight the radar–lidar complementarity in their algorithm. Comparative studies of microphysical properties in cirrus clouds between DARDAR products and in situ observations have shown the reliability of DARDAR retrievals [48] for a vast range of cloud properties found in the Arctic.

DARDAR V1.0 data are used in this study. Four DARDAR tracks as close as possible (in terms of Euclidean distance) to the four synchronous ISDAC flights have been selected to validate our model simulations over a larger area at the synoptic scale (Figure 1): T12 (4 April 2008, 06:39:28–07:25:38UTC), T13 (5 April 2008, 06:36:58–07:027:43 UTC), T21 (14 April 2008, 06:36:30–07:27:51 UTC) and T29 (24 April 2008, 06:38:53–07:26:42 UTC) associated to flights F12, F13, F21 and F29, respectively. The DARDAR mask was used to identify ice cloud phase [23]. On each satellite overpass, the number of crystals Ni is retrieved from both DARDAR-derived IWC and Re according to Jouan et al. [23]. The comparative analysis is focused in the part of the cloud where both radar and lidar observations are detecting ice to better constrain retrievals (above 4 km).

4. Results and Discussion

This section is dedicated to present comparisons of model simulations (WRF, WRF_p and WRF_np) against observations, followed by a discussion of the results. The vertical profiles of the simulated and in situ observed variables are compared over a region of 10 by 10 km centered around the location of ISDAC along the flight traks. ISDAC in situ measurements have been averaged every 20 s, corresponding to a vertical resolution of ~45 hPa, during ascents and descents through clouds. Modeled WRF outputs are linearly interpolated to the pressure levels of these observations and temporally averaged over a three hour period encompassing the area of ISDAC flights. DARDAR products are averaged over the aircraft flight range identified in Figure 1 (box between (65°N, 142°W) and (72°N,162°W) and serve as a basis for comparison against model results over the same region. On each vertical cloud profile, we examine cloud microphysical properties in three broad height regions, defined as follows: low-level clouds when pressure is higher than 650 hPa (altitude lower than ~3.4 km), high-level clouds when pressure is lower than 500 hPa (altitude higher than ~5.8 km), and mid-level clouds in-between, following Mason et al. [49].

4.1. Meteorological Profiles

Figure 3 and Figure 4 present the comparison of the observed and simulated (WRF, WRF_np and WRF_p) vertical profiles of T and RH_i for the four flights F12 (A), F13 (B), F21 (C) and F29 (D). The measurements variability is also reported by the standard deviation of all points at each altitude level. As expected, we do not observe a significant difference on T and RH_i between the different simulations. The temperature profile is indeed not significantly affected by the modification of IWC due to a change in the ice nucleation parameterization. The observed temperature profile continuously decreases between 800 hPa and 400 hPa by ~30 °C. The model reasonably underestimates the temperature for the four flights when the pressure is higher than 650 hPa with a bias lower than 2 °C.

The lowest temperatures at the top of the clouds are also well reproduced, except along F21 flight. For that flight, the three runs underestimate T by ~5 °C and RH_i by ~40% at the top of the cloud. The vertical profile of RH_i is rather constant with height in the different simulations, whereas it is increasing in the observations. Although the mean value is very well represented along the four flights, the model is unable to represent the substantial enhancement of RH_i with height by ~20–40% between 800 and 400 hPa. The differences between the different runs in terms of RH_i are not really significant and are the results of second-order feedbacks; a modification of the ice nucleation rate will affect the total number of ice crystals Ni but also the rate of crystal growth by riming or by condensation through the Wegener–Bergeron–Findeisen effect. Subsequently, the growing ice crystals consume water vapor faster than it is made available by the cooling of the air. Consequently, the supersaturation with respect to ice decreases in the model due to this process. The higher the IWC, the lower the RH_i. Along the four flights, RH_i is therefore lower in the WRF_np simulation with higher IWC, than in the WRF_p run with lower IWC.

4.2. Spatial Structure of IWC in Large Scale Clouds

The vertical profiles of IWC obtained from the three simulations extracted along tracks T12 (A), T13 (B), T21 (C) and T29 (D) (Figure 1) are compared to those derived from the DARDAR observations. The corresponding curtain plots are given in Figure 5 and Figure 6. The red rectangles represent the flight regions (Figure 1). Along tracks, DARDAR data show an IWC vertical profiles continuously decreasing from 10⁻¹ g/kg at 2 km to less than 10⁻² g/kg above 9 km. Comparisons shown in Figure 5 and Figure 6 indicate that the vertical structures of IWC are broadly captured by both WRF_np and WRF_p simulations at all vertical levels, with a RMSE lower than 10⁻² g/kg. WRF run systematically underestimates the IWC for mid- and low-level clouds with a low mean bias of –3.10⁻² g/kg and may miss some clouds presenting low IWC values. On the contrary, WRF_np run always simulates clouds with larger IWC values than those obtained from WRF_p (37% difference) and to a larger extent than from WRF-derived values (ratio close to 3).

The curtain plot of IWC values extracted along T12 (Figure 5A) shows that the horizontal structure with two denser columns near 67°N and 70°N is well captured. The agreement is not as good for T13 (Figure 5B). In that case, the vertical cloud structure retrieved from satellite observations is more representative of two layers, except near the center of the domain (about 69°N) where a more homogeneous profile can be noticed on the vertical column, likely caused by precipitation.

Along T21 and T29 tracks (Figure 6A,B), WRF simulation (IWC~10⁻⁷-10⁻⁶ kg/kg) strongly underestimates DARDAR-derived IWC (IWC~10⁻⁵–10⁻⁴ kg/kg), as clouds are even not developed in the model, except below 4km for T21. WRF simulation is representative of a clear-sky situation between 66°N and 74°N for T29 as it totally fails to model clouds at all vertical levels. In all case studies, WRF run hardly represents low IWC values as horizontal structures evidenced by satellite observations are poorly reproduced.

Along all tracks, WRF_np and WRF_p simulations better capture the IWC in high- and mid-level clouds (e.g., between 4 and 10 km, corresponding to temperatures below −25 °C) but underestimate it in low-level clouds, below 4 km with respect to DARDAR retrievals. Our analysis reveals that the use of the CNT primary ice nucleation (WRF_p and WRF_np runs) clearly improves the representation of the IWC in the reference clouds. In addition, the difference on the simulated IWC using a parameterization with different contact angles (12° and 26°) (<0.5 10⁻² g/kg) is smaller than the discrepancy between the base WRF run and observations (low bias of ~−3.10⁻² g/kg).

4.3. IWC Profiles Along Flight Tracks

Comparison of the observed (both in situ measurements and DARDAR products) and simulated (WRF, WRF_np and WRF_p) local and average vertical profiles of IWC along F12 (A), F13 (B), F21 (C) and F29 (D) flights is presented in Figure 7. In addition to simulated local profiles extracted in the vicinity of in situ observations, Figure 7 also shows modeled and satellite-derived profiles averaged (labeled_mean) over the area of the flights (red rectangles in Figure 5 and Figure 6). For the sake of clarity, Figure 7 only represents the modeled average profile showing the best statistical score according to results discussed in Section 4.2 at larger scale. Along F12 and F13 flights, the WRF_np_mean is therefore represented, whereas the WRF_p_mean profile is shown along F21 and F29 flights. The following figures (Figure 8B, Figure 9B, Figure 10B, Figure 11B, Figure 12B, Figure 13B, Figure 14B and Figure 15B) will also use this convention. Significant differences between in situ and average DARDAR (DARDAR_mean) profiles are found. They strongly depend on times and locations where clouds have been preferably sampled by in situ measurements.

As expected from curtain plots shown in Figure 5 and Figure 6, WRF run does not predict any noticeable IWC for the F12 and F21 flights. IWC is strongly underestimated in the WRF simulation (−10⁻² g/kg) along F13 (−4.10⁻² g/kg) and F29 (−8.10⁻² g/kg) flights at all vertical levels with respectively two and one orders of magnitude. This quantifies the poor performance of the base run of the WRF model in representing ice clouds with low IWC (about 0.01 g/kg or below) and shows the general underestimation of IWC of the base parameterization. In contrast to the base run (WRF), the WRF_np and WRF_p simulations significantly improve the simulated IWC for all flights. Along F12 flight, IWC derived from the WRF_p simulation remains much too low. Only the IWC predicted by the WRF_np run turns out to be close to the observations at the cloud top (−8.10⁻³ g/kg), indicating an improved performance of the WRF_np run in that “not polluted” air mass. Apart from that case, both WRF_np and WRF_p simulations give results within a factor of 2–5 at low IWC, and a factor of 2 at large IWC. The IWC is systematically predicted to be larger in the WRF_np than in the WRF_p simulation. A smaller contact angle in the WRF_np run indeed tends to decrease the critical Gibbs free energy to form ice embryos (Equation (3)), hence leads to a higher nucleation rate of ice crystals. This is in agreement with the larger IWC values found in the WRF_np simulation and associated with a slightly lower RH_i.

The averaged DARDAR profiles in the region of the flights (DARDAR_mean profiles) are of the same order of magnitude as the in situ measurements. The agreement is especially very good over the full range of altitudes for F21, above 4 km for F12 and below 4 km for F13. It is however strongly overestimated when the IWC obtained from airborne measurements is very low (< 10⁻³ g/kg below 4 km for F12 flight) and underestimated by 1-2 orders of magnitude when the IWC measured by in situ probes is larger (10⁻²–10⁻¹ g/kg). Different times and locations between the aircraft and the satellites together with natural cloud variability might explain such discrepancies. In addition, Dong et al. [50] note that generally DARDAR retrievals are thought to slightly overestimate in situ IWC measurements. Nevertheless, DARDAR_mean profiles generally report the slight decrease of IWC from 10⁻² to 10⁻¹ g/kg with altitude in accordance with in situ observations. Given the challenge to faithfully evaluate IWC, the average vertical profiles (WRF_np_mean for F12 and WRF_p_mean for F21) show a good performance and both simulate DARDAR_mean with a RMSE close to 10⁻² g/kg, i.e., a relative error lower than 10%. For F13 case, both the model and the satellite derived IWC follow the same vertical variability with a distinct minimum of approximately 10⁻³ g/kg between 4.5 and 6 km. The corresponding Pearson correlation coefficient of WRF_np_mean versus DARDAR_mean is fairly good (0.65), yet the model underestimates the satellite profile by one order of magnitude. On the F29 case, even if the model does a good job in representing the in situ observations, its performance is much worse when considering the comparison to the satellite product showing the difficulty to have a clear and coherent picture from independent observations. A correlation coefficient of 0.30 is found between WRF_p_mean and DARDAR_mean and the simulation underestimates IWC by one order of magnitude at the top of the cloud, with very low IWC in the lower troposphere up to 650 hPa.

In summary, results at the local (in situ airborne measurements) and regional (remote sensing observations) scales both confirm the very large discrepancies between the original microphysical scheme of the WRF model and the observations in the representation of the IWC in the Arctic. Building a parameterization based on the contact angle improves substantially the vertical structure of ice clouds in terms of layering and IWC. Section 4.4 will further investigate the microphysical parameters of ice clouds in TIC-1 and TIC-2 clouds.

4.4. Microphysical Parameters of Ice Clouds

In this section, we examine the impacts of the various microphysical scheme in terms of their resulting number concentration (Ni) and ice crystal radius (Ri). The profiles derived from the different simulations are compared with airborne ISDAC campaign and DARDAR satellite data. Section 4.4.1 presents results when the sampled vertical profiles have been influenced by clean air masses (F12 and F13 case studies). Section 4.4.2 is devoted to the influence of polluted air masses (F21 and F29 case studies).

4.4.1. Clean Air Masses

Figure 8A and Figure 9A present the cross-sections of DARDAR-derived and simulated ice crystal number concentration (Ni) extracted along T12 and T13, respectively. Figure 8B and Figure 9B represent the vertical profiles of Ni simulated and measured during ISDAC F12 and F13 flights, respectively. Figure 8B and Figure 9B also include mean vertical profiles of Ni from DARDAR (DARDAR_mean) and from WRF_np (WRF_np_mean) averaged over the box between (65°N, 142°W) and (72°N, 162°W) area reported in Figure 8A and Figure 9A. Note that the DARDAR_mean predicts any variability on the vertical structure of the clouds with value around 200/L and 50/L on T12 and T13, respectively.

Although the simulations show some skill to represent the vertical cloud structure (Figure 8A and Figure 9A), large discrepancies exist on the quantitative values of Ni in large cloud systems. On T12 track (Figure 8A), WRF Ni is strongly underestimated at all vertical levels by three orders of magnitude. The skill of the WRF_np and WRF_p simulations improves significantly over the original WRF. WRF_p reproduces the vertical structure of Ni relatively well, except in low-level clouds where it is often underestimated. WRF_np shows an improvement of the representation of Ni in low-level clouds, but overestimates Ni values at higher altitudes by about a factor of 2. Along T13 track (Figure 9A), the performance of WRF_np and WRF_p simulations is rather good, but overestimating Ni of high-level clouds in regions near 53°N, 54°N, 67°N and 72°N. WRF run is characterized by a significant underestimation with a RMSE of 158 L⁻¹ at low level, and an overestimation in high level thin ice clouds between 63°N and 73°N. In the lowest layers, the best performance is given by the WRF_np results with a RMSE of 114 L⁻¹.

The airborne ISDAC vertical profile for F12 (Figure 8B) shows an average value of 150 L⁻¹ in broad agreement with the average DARDAR retrieval on T12 at altitudes below 6.4 km (pressures larger than 450 hPa). It is increasing to about 10³ L⁻¹ within a 550 m thin layer (~ 50 hPa) at cloud top, but remains close to DARDAR on average in this layer. The WRF_np profile averaged along the DARDAR track overestimates Ni by a factor of 5, but the average WRF_p profile underestimates it by about the same factor. When compared to ISDAC measurements averaged over the upper cloud layer (450–400 hPa), all schemes underestimate Ni except the WRF_np scheme within a factor of 2, still smaller than the factor 10 to 50 observed for WRF and WRF_p. In the layer below (between 450 and 550 hPa), all schemes underestimate Ni values within two orders of magnitude except for WRF_np. Average WRF_np results are overestimated by a factor of 2, but still remain better than the WRF_p scheme results. The vertical microstructure of the cloud with Ni increasing from 10⁻³ L⁻¹ at 2 km to 10³ L⁻¹ at 6.4 km is therefore improved in the WRF_np run.

The observed vertical profile of Ni for the F13 ISDAC flight (Figure 9B) is rather constant with altitude, and ranges from 70 to 200 L⁻¹, even in the upper part of the cloud layer (Figure 9A). The average profile of DARDAR (DARDAR_mean) presents a fair correlation with the ISDAC mean profile (R² ~0.45), even if a systematic weak underestimation is seen (mean bias of 37 L⁻¹). The original WRF simulation strongly underestimates Ni by at least two orders of magnitude in the upper part of the cloud, above an altitude of 4 km (600 hPa). In contrast, our parameterization with WRF_np run relatively overestimates Ni by a factor of 3 on average in this part of the cloud. WRF_p simulation reproduces Ni relatively well at the cloud top, but remains within a factor 5 smaller on average. Under 4 km, WRF run does not detect any noticeable Ni and strongly underestimates observations by four orders of magnitude. WRF_np and WRF_p simulations underestimate Ni in this altitude range by two orders of magnitude, the best performance being given by WRF_np. In comparison to DARDAR_mean profile, the skill of the WRF_np_mean is rather good except in a thin layer between 500 and 550 hPa. This is mostly because the model properly represents two distinct cloud layers (Figure 9A), which cannot be discerned when observations are averaged in the vicinity of the flights.

Figure 10A and Figure 11A present the cross-sections of DARDAR-derived (Re) and simulated (mean radius of particles Ri) ice crystal size extracted along T12 and T13, respectively. Figure 10B and Figure 11B represent the vertical profiles of Ri obtained during ISDAC F12 and F13 flights, respectively, including average DARDAR (DARDAR_mean) and WRF_np_mean profiles. Along T12 track, DARDAR predicts an effective radius continuously decreasing from 80 µm at 2 km to less than 50 µm above 9 km. WRF respectively strongly underestimates Ri by a factor 5 above 6.5 km and overestimates it under 6.5 km by a factor of 2. Using a contact angle of 12° leads to a reduced size of ice crystals. WRF_p and WRF_np indeed underestimate Ni with a RMSE of 20 and 25 µm, respectively. The average crystal size measured in situ is about 20-40 µm (Figure 10B). Along F12 flight, WRF and WRF_p runs give Ri close to 0 through the whole cloud, in agreement with negligible IWC simulated in those runs (Figure 7A). In contrast, all values given by WRF_np, WRF_np_mean, DARDAR_mean and in situ observations are in better agreement. A very small mean bias of 7 µm is noticed between WRF_np and the ISDAC profiles. The underestimation is a bit larger in the vertical column between WRF_np_mean and DARDAR_mean, but the performance of this parameterization can be considered as much improved given the high sensitivity of the variables.

Despite an underestimation of the size of ice crystals along T13 track (Figure 11A), the discrepancies between the three simulations are not significant. The smallest error is given by WRF_np with a RMSE of 44 µm, followed by WRF_p (RMSE of 54 µm). The largest error is obtained with the standard version of the WRF model, reaching a RMSE of 57 µm. Values extracted along F13 flight (Figure 11B) are somewhat larger than along F12 flight (about 30-50 µm). Overall, a much better performance is found with WRF_np on F13 case in cloud layers above 4 km, (P <620 hPa) associated to clean air masses. Below this altitude level, there is an overestimation of Ri by factor 2 in all simulations. It is linked to the underestimation of the number of ice crystals in this region. A fairly good performance is obtained using the WRF_np_mean profile with a RMSE close to 15 µm. As a summary, the Ni is better represented in clouds under clean air influences when the WRF_np run is used.

4.4.2. Polluted Air Masses

Figure 12A and Figure 13A present the cross-sections of DARDAR-derived and simulated Ni extracted along T21 and T29, respectively. Figure 12B and Figure 13B represent the vertical profiles of Ni obtained during ISDAC F21 and F29 flights, and averaged over T21 and T29 tracks, respectively.

All simulations demonstrate a skill to represent the vertical structure of the presence of ice crystals in the atmospheric column (Figure 12A and Figure 13A). However, large discrepancies exist on the quantitative values of Ni in these two large cloud systems. On T21 and T29 tracks (Figure 12A, Figure 13A), WRF_np largely overestimates Ni at all vertical levels, with a RMSE of 303 and 1057 L⁻¹, respectively. In contrast, WRF simulation presents a strong underestimation of Ni for these two case studies: the RMSE is 100 and 165 L ^–1 for T21 and T29, respectively. Only WRF_p reproduces relatively well Ni over T21 (RMSE of 9 L⁻¹) and T29 (RMSE of 15 L⁻¹) tracks. The improvement from the WRF_p parameterization with respect to observations is indeed essentially dependent on the contact angle $\theta$. When $\theta$ increases from 12° to 26°, then $\cos \left(\theta \right)$ decreases and $f\left(\cos \left(\theta \right)\right)$ increases. It follows that ΔG* increases, meaning that the energy barrier for nucleation becomes stronger and therefore difficult to reach. As a consequence, the nucleation rate J decreases and Ni decreases and becomes closer to observations. It should be noted that WRF_p run does not detect a TIC-1 cloud at latitudes lower than 56°N, while this cloud is present in both WRF and WRF_np simulations. This confirms the inability of the WRF_p run and the rather good performance of the WRF_np prescription to represent ice clouds in clean air masses (Section 4.4.1). According to Jouan et al. [23], it is consistent with the fact that this southern region of our domain is not affected by the polluted air mass in both F21 and F29 cases further North of Alaska. Satellite observations from the Ozone Monitoring Instrument (OMI) support that result by showing an increased concentration of SO₂ on 15 April (F21) and on 25 April (F29) over the North Slope of Alaska [20]. In this region under the influence of volcanic aerosols, the best agreement between model and observations is provided by the WRF_p simulation (Figure 12A and Figure 13A).

Both TIC-2 clouds observed along F21 and F29 flight tracks are characterized by a small concentration of ice crystals varying between 1 and 30 L⁻¹ from cloud top to cloud base (Figure 12B and Figure 13B). The average profile derived from the spatial observations (DARDAR_mean) is about a factor of 2 larger than Ni observed during the ISDAC on both F21 and F29 flights. WRF and WRF_np runs respectively strongly underestimate and overestimate Ni by one to two orders of magnitude against in situ observations (Figure 12B and Figure 13B). Only WRF_p, suitable for acid aerosol, simulates reasonably well Ni (within a factor of 2 on average) and produces the best statistics along F21 flight (Figure 12B). On Figure 13B, comparisons give similar results, but this time the WRF_p simulation overestimates to a larger extent (factor 10) Ni detected from the F29 measurements. WRF-derived Ni are the closest to the observed values, indicating that our parameterization, although generally improving the results, still need further investigation to better understand the sources of remaining discrepancies.

Figure 14 and Figure 15 present observed cross sections and vertical profiles of Re and Ri on April 15 and 25. Figure 14B and Figure 15B include the average DARDAR (DARDAR_mean) and WRF_p_mean profiles. Along F21 flight, ISDAC profile shows a rapid growth of the ice crystal size at cloud top (~60 µm) gradually increasing to about 100 µm at cloud base. Those values are a factor of 2 larger than radii obtained along F12 and F13 flights. This F21 case refers to TIC-2 clouds characterized by a smaller number of larger ice crystals. An excellent agreement is found between DARDAR and ISDAC profiles. Along T21 track and F21 flight, all the simulations underestimate DARDAR and ISDAC observations in the range of uncertainty except WRF base run which strongly underestimates clouds at all vertical levels. This is caused by the corresponding large underestimation of IWC in the base WRF simulation (Figure 7C). Simulations based on a contact angle also underestimate observed crystal radii with a low bias of –42 and –28 µm for WRF_np and WRF_p runs, respectively. The same result is obtained at large scale: WRF_p_mean underpredicts DARDAR_mean Re values with a modest RMSE of 22 µm for T21 case.

Values of Ri are even larger along F29 flight. They indeed range from 90 µm at cloud top to 170 µm at cloud base. The WRF_np simulation strongly underestimates observed in situ measurements with a mean bias of –69 µm. A better agreement is found with the WRF run (bias of -29 µm) and WRF_np run (bias of -50 µm), (Figure 15B). DARDAR vertical profiles give smaller values in the atmospheric column than ISDAC observations. The DARDAR_mean profile present Re values from 50 to 90 µm. A good agreement is obtained with the WRF_p_mean simulation (RMSE ~35 µm). The comparison of the vertical profiles of Ri and Re retrieved along F21/F29 flights and T21/T29 tracks suggests that the best performance is obtained in the WRF_p simulation, known to be more representative of air masses containing more acidic aerosols and therefore TIC-2 clouds.

5. Summary and Conclusions

The objective of this study was to improve the heterogeneous ice nucleation parameterization in the Weather Research and Forecast (WRF) model. Several cases combining in situ measurements and remote sensing observations have been identified in clean and polluted situations to perform comparative modeling tests. Two types of ice clouds (TIC) have been distinguished, according to classifications based on extensive measurements from ground-based sites and satellite remote sensing in the Arctic during the polar night and early spring. TIC-1 clouds are composed by non-precipitating very small (unseen by radar) ice crystals whereas TIC-2 clouds are characterized by a low concentration of large precipitating ice crystals. It is hypothesized that TIC-2 formation is linked to the acidification of aerosols, which inhibit the ice nucleating properties of IN. As a result, the IN concentration is reduced in these regions, resulting to a smaller concentration of larger ice crystals. In this work, we successfully simulated cloud distributions and structures corresponding to four well documented case studies as references, obtained during ISDAC. The observed clouds were characterized by approximately the same order of magnitude of IWC (10⁻² g/kg), but differed by the origin of air masses. Numerical simulations with a modified microphysical parametrization (WRF_np in pristine air masses and WRF_p in polluted environments) have been evaluated against in situ aircraft data and satellite remote sensed observations in the vicinity of the flights. The main findings in this study are the following:

• The original WRF model fails to simulate properly most ice clouds, and without discrimination as compared to WRF_np and WRF_p parameterizations. IWC is strongly underestimated in the original WRF simulation along all flights at all vertical levels with several orders of magnitude. This reveals the poor performance of the base run of the WRF model in representing ice clouds with low IWC (about 0.01 g/kg or below).
• Our analysis also reveals that the use of the CNT primary ice nucleation (WRF_p and WRF_np runs) clearly improves the representation of the IWC against both local and large scale references in Arctic thin ice clouds. In addition, the difference on the simulated IWC using a parameterization with different contact angles (12° and 26°) is smaller than the discrepancy between the base WRF run and observations. The results of our modifications in the microphysical scheme in terms of number concentration (Ni) and ice crystal radius (Ri) also show that the use of WRF_np and WRF_p simulations is improving substantially the results both in TIC-1 and TIC-2 clouds, representative of clean and polluted air masses, respectively.
• The WRF_np simulation generally results in smaller ice crystals (Ri~10–80 µm) with larger concentrations (Ni~50-200 L⁻¹). The cloud microstructure observed during flights in clean air masses corresponds to TIC-1 clouds characteristics, as aimed for. In TIC-1 case, it was suggested that the dominant nucleation process was the deposition ice nucleation due to an atmosphere not saturated with respect to liquid water. Compared to the original WRF base case, WRF_np improves IWC by factor 2, Ni up to 1 order of magnitude and Ri by 40% in terms of bias. In contrast, the WRF_np scheme does not represent the observed microstructure of ice clouds observed in polluted conditions (TIC-2 clouds). These results were expected since the WRF_np scheme has been specifically developed to reproduce ice nucleation in a pristine environment.
• On the other hand, the WRF_p scheme can simulate small concentrations (Ni<40 L⁻¹) of larger crystals (Ri~50-180 µm). In the TIC-2 case, either immersion or condensation freezing of haze droplets (coated IN) was hypothesized, due to larger relative humidity with respect to ice, close to saturation with respect to liquid water. Low concentrations of IN combined to the high supersaturated air with respect to ice, lead to an explosive growth of ice crystals by water vapor diffusion. Compared to the original WRF base, WRF_p improves IWC by factor 4, and reduce Ni RMSE by 60% but increase Ri bias by 20%. The WRF_p run was parameterized, based on laboratory experiments [38,39], to reproduce ice nucleation in polluted environments (TIC-2 clouds). We find that it corresponds relatively well to the observed microstructure of the ice clouds observed in air masses under the influence of volcanic plumes. These polluted cases were expected to lead to IN coated with sulfuric acid air masses conditions [20]. Inversely the WRF_p simulations do not reproduce clouds formed in a pristine environment (TIC-1 clouds). The deference was hypothesized to be closely related to the nucleation process [6].
• Results at a larger scale from a comparison to satellite data (DARDAR) confirm the conclusions on vertical profiles sampled during ISDAC, suggesting that our results are not only a particular feature obtained during the ISDAC campaign at a specific location, but can be applicable on general cases of thin ice clouds in the Arctic.

Our parameterization is nevertheless limited in space because it is tested for a fixed concentration of aerosols while, in a dynamic evolution, the concentration of aerosols varies in time, in space and with the acid concentration. A following study will account for variable aerosol and acid concentrations in order to couple the microphysical scheme to a prognostic chemical module. The improved parameterization will therefore take into account the temporal and spatial variation of the aerosol concentrations and their degree of acidity using the proper contact angle. Our parameterization is also limited to the ice nucleation on dust particles only. Further work could consider heterogeneous ice nucleation on other types of particles, e.g., primary biological particles (bacteria, fungal spores and pollens). Taking into account the mixed composition of aerosols to develop specific parameterizations of ice nucleation could improve simulations of ice clouds where aerosol plumes of various types are observed.