Aerosol impacts on the entrainment efficiency of Arctic mixed-phase convection in a simulated air mass over open water

Springtime Arctic mixed-phase convection over open water in the Fram Strait as observed during the recent ACLOUD (Arctic CLoud Observations Using airborne measurements during polar Day) field campaign is simulated at turbulence-resolving resolutions. The first objective is to assess the skill of large-eddy simulation (LES) in reproducing the observed mixed-phase convection. The second goal is to then use the model to investigate how aerosol modulates the way in which turbulent mixing and clouds transform the low-level air mass. The focus lies on the low-level thermal structure and lapse rate, the heating efficiency of turbulent entrainment, and the low-level energy budget. A composite case is constructed based on data collected by two research aircraft on 18 June 2017. Simulations are evaluated against independent datasets, showing that the observed thermodynamic, cloudy, and turbulent states are well reproduced. Sensitivity tests on cloud condensation nuclei (CCN) concentration are then performed, covering a broad range between pristine polar and polluted continental values. We find a significant response in the resolved mixed-phase convection, which is in line with previous LES studies. An increased CCN substantially enhances the depth of convection and liquid cloud amount, accompanied by reduced surface precipitation. Initializing with the in situ CCN data yields the best agreement with the cloud and turbulence observations, a result that prioritizes its measurement during field campaigns for supporting high-resolution modeling efforts. A deeper analysis reveals that CCN significantly increases the efficiency of radiatively driven entrainment in warming the boundary layer. The marked strengthening of the thermal inversion plays a key role in this effect. The low-level heat budget shifts from surface driven to radiatively driven. This response is accompanied by a substantial reduction in the surface energy budget, featuring a weakened flow of solar radiation into the ocean. Results are interpreted in the context of air–sea interactions, air mass transformations, and climate feedbacks at high latitudes.

1 Introduction

The ongoing accelerated warming of the Arctic climate involves various processes and feedback mechanisms, many of which are still poorly understood. Recent research has highlighted the role of warm-air intrusions (Bennartz et al., 2013; Pithan et al., 2018) as well as the lapse rate feedback (Pithan and Mauritsen, 2014; Lauer et al., 2020). Clouds play a sophisticated role in these mechanisms. For example, cloud presence in warm-air intrusions significantly affects the downward longwave radiative flux at the surface (Liu et al., 2018). Radiative cooling at the liquid cloud top also causes turbulence, which in turn drives entrainment that counteracts larger-scale subsidence, together maintaining the low-level inversion and lapse rate (e.g., Neggers et al., 2019). Arctic clouds form through a variety of processes acting on a broad range of scales (e.g., Mauritsen et al., 2011; Morrison et al., 2012). Gaining further insight has motivated intense research on Arctic air masses and clouds, including a number of field campaigns at high latitudes (Perovich et al., 1999; Tjernström et al., 2012; Knudsen et al., 2018; Wendisch et al., 2019; Shupe et al., 2021).

This study focuses exclusively on mixed-phase convective clouds in relatively stagnant air masses over open water. Previous studies on marine cold-air outbreaks (CAOs) have shown that the strong surface–atmosphere temperature difference over open water can drive intense cloudy convection, which is efficient in vertically mixing the lower atmosphere (Chlond, 1992; Atkinson and Wu Zhang, 1996; Müller et al., 1999; Gryschka and Raasch, 2005; Fletcher et al., 2016). In comparison, weaker convection in more stagnant air masses has received far less attention. However, such air masses occur frequently and might occur even more in a warmer future Arctic featuring a slower polar jet stream (Screen et al., 2013; Barnes and Screen, 2015). Slow-moving air masses also have much more time to adjust to local conditions, which potentially makes the vertical mixing more efficient. Finally, the ongoing shift in Arctic climate is arguably strongest felt in areas where the sea ice disappears (Liu et al., 2012; Overland et al., 2014). The marine areas adjacent to the sea ice also act as gateways for injections of aerosol (Browse et al., 2014; Ito and Kawamiya, 2010), moisture, and heat (Vázquez et al., 2016; Rinke et al., 2017; Pithan et al., 2018) into the high Arctic.

These reasons motivate taking a closer look at such stagnant marine air masses, in particular concerning clouds and turbulent mixing. Recent years have seen an increased use of large-eddy simulation (LES) to study these processes, often based on Arctic field campaign data. LES can supplement the observational data record and act as a virtual research laboratory (e.g., Ovchinnikov et al., 2014). At the same time, independent measurements of mixed-phase cloud properties can be used to evaluate the simulations (Neggers et al., 2019; Kretzschmar et al., 2020; Ruiz-Donoso et al., 2020). This approach has led to demonstrable progress in understanding Arctic clouds. A few recent papers have investigated aerosol impacts on mixed-phase clouds (de Roode et al., 2019; Stevens et al., 2018). However, no LES study has yet examined these impacts in more stagnant air masses over open water near the ice edge. In addition, no LES study has yet used in situ cloud condensation nuclei (CCN) measurements at the cloud level to constrain simulations of this cloud regime.

The main science goal of this study is to use LES to gain more insight into how and to what extent aerosol variations in a slow-moving Arctic air mass over open water can modulate its transformation by low-level turbulent/convective mixing and clouds. Of particular interest are aerosol impacts on the efficiency of radiatively driven entrainment in warming the boundary layer, as well as the associated shifts in the heat budgets of the boundary layer and the surface. While the entrainment efficiency has previously been investigated for warm turbulent clouds in the subtropics (Stevens et al., 2005), this is not yet the case for mixed-phase clouds at high latitudes. What is also still unclear is how CCN concentrations might affect this efficiency (Garrett et al., 2002; Douglas and L'Ecuyer, 2021). Given the importance of low-level warming and aerosol variations in Arctic amplification, in particular in the context of the lapse rate feedback, gaining more insight into this process and its sensitivities is crucial.

To achieve these objectives a composite LES case is constructed based on the extensive data collected by the Polar 5 and Polar 6 aircraft of the German Alfred Wegener Institute in the Fram Strait during the ACLOUD campaign (Arctic CLoud Observations Using airborne measurements during polar Day) (Wendisch et al., 2019). Research flight RF20 on 18 June 2017 sampled mixed-phase clouds as embedded in a stagnant air mass off the sea ice edge. The sampled convection was significant but still relatively weak compared to typical cold-air outbreak conditions. The boundary conditions and large-scale forcings for the simulations are based on weather model data, while the initial state is based on dropsonde data as well as in situ aerosol data at the cloud level. First the control experiment will be evaluated against independent aircraft data on clouds and turbulence, seeking agreement on basic bulk properties that are well observable. Based on this control run, sensitivity tests are then performed on CCN and levels in the air mass. Compared to previous LES studies of this kind, a much broader CCN range is covered to capture the large observed variation in Arctic air masses between pristine polar and polluted continental values.

Section 2 describes details of the ACLOUD field campaign, including the weather situation, the research flights, and the observational datasets collected. The model configuration adopted in this study is described in detail in Sect. 3, including the LES code, the treatment of microphysics, the case configuration including forcings and boundary conditions, and the experimental setup. The presentation of the results is subdivided into two parts. Part I describes the basic behavior of the control experiment, including an evaluation of cloud and turbulence statistics against ACLOUD observational datasets (Sect. 4). Part II then focuses on the aerosol sensitivity experiments (Sect. 5). The obtained results are further interpreted in Sect. 6, and the main conclusions and outlook are summarized in Sect. 7.

2 Observations

2.1 The ACLOUD field campaign

The ACLOUD field campaign took place from 23 May to 26 June 2017 in the vicinity of Svalbard. ACLOUD and its sister campaign PASCAL (Physical feedbacks of Arctic planetary boundary level Sea ice, Cloud and AerosoL; Macke and Flores, 2018; Wendisch et al., 2019) were part of the ongoing (AC)³ research program (Arctic Amplification: Climate Relevant Atmospheric and Surface Processes and Feedback Mechanisms; Wendisch, 2017). Both campaigns focused on clouds in the lower troposphere in the northern Fram Strait during Arctic spring. ACLOUD featured collocated airborne observations (Ehrlich et al., 2019b) performed by the aircraft Polar 5 (P5) and Polar 6 (P6) of the German Alfred Wegener Institute (Wesche et al., 2016). An overview of the synoptic conditions during ACLOUD is provided by Knudsen et al. (2018). Airborne observations were made during a wide range of cloud conditions, including both stably stratified and convective regimes.

2.2 RF20

This study exclusively focuses on research flight RF20 by the P5 and P6 aircraft on 18 June 2017 (see Fig. 1). As described by Knudsen et al. (2018), on this day the mid-troposphere experienced drying. High cirrus clouds were initially present but disappeared during the day. The air mass in the Fram Strait was slow-moving and situated over relatively warm open water.

Figure 1 (a) The location of the RF20 mission over Fram Strait west of Svalbard (b) with a close-up of the target area: MODIS TERRA true-color image at 250 m effective resolution during RF20 on 18 June 2017. The flight tracks of Polar 5 and Polar 6 are shown in red and yellow, respectively. For P5 the waypoints (C), RV Polarstern (PS) and Longyearbyen airport (LYR) are also indicated, with the dropsonde (DS) locations shown as blue dots. The 24 h air mass trajectory intersecting with DS04 is shown as a solid green line. The target domain including the southern racetrack section is indicated by the dotted orange box, of which the simulated part is indicated by the dotted green box. MODIS data obtained through NASA Worldview (https://worldview.earthdata.nasa.gov/, last access: 22 July 2022).

Figure 1 shows that the cloud situation in the Fram Strait as encountered by the aircraft was relatively complex. A rough north–south regional division in cloud character can be made. In the north, over the sea ice and its margin, the clouds were absent or very thin, allowing for good visibility of the sea ice from the satellite and the aircraft. In the western and middle parts of the Fram Strait the clouds were thicker but still only weakly convective, visible in Fig. 1 as vague but not completely opaque cloud patches. In the southern part the clouds were truly convective, being broken, thicker, and more opaque.

The P5 and P6 aircraft followed counterclockwise flight paths from their base at Longyearbyen airport (LYR) and visited these three regimes consecutively. This study focuses exclusively on the convective clouds in the southern areas, as sampled during the last eastbound flight leg between waypoints C5 and C6. In this section, also referred to as the “southern racetrack”, both aircraft flew back and forth between C5 and C6. While P5 doubled back once and maintained altitude above the boundary layer inversion (at about 1.5 km height), P6 doubled back twice, staying below inversion height and maintaining constant altitude for five brief flight segments. In situ in-cloud measurements were made by P6 during this period. The enhanced and targeted sampling during the southern racetrack sections, as well as the occurrence of significant mixed-phase convection, motivates adopting this area as the target domain of this study, as indicated by the orange box in Fig. 1.

2.3 Observational datasets

The observational data from the ACLOUD campaign used in this study are summarized in Table 1. Figure 2a shows detailed measurements of the clouds in the target area as obtained with the Microwave Radar/radiometer for Arctic Clouds radar on board P5 (MiRAC; Mech et al., 2019). Because of the doubling back between C5 and C6, the same cloud structure appears three times, once in each panel. Cloud top height varies significantly along the flight track in this area, which is typical of broken convective cloud fields. The maximum cloud top height is at approximately 1.4 km, which is consistent with the thermal inversion height visible in the DS04 profile (see Fig. 3a) and the Airborne Mobile Aerosol Lidar (AMALi; Stachlewska et al., 2010) measurements (also included in Fig. 2). The MiRAC flight sections between 7–8.5^∘ E feature significant but narrow convective precipitation that also reached the surface. This area was visited by P5 three times, at around 16:20, 16:35, and 17:00 UTC. The DS04 dropsonde was also launched into this area. The freezing level is situated at about 350 m height, which is well below the maximum cloud top.

Figure 2 Time–height cross section of P5 MiRAC radar reflectivity Z_e during RF20 (contour shading) and the indicated liquid cloud top from AMALi (black dots). The displayed longitude range corresponds to the orange target domain as shown in Fig. 1, while the lime green horizontal line indicates the simulated domain. The black arrows indicate the flight direction of the aircraft for each leg, with the start and end times indicated at the sides. The location of the dropsonde DS04 is indicated by the dotted dark-green vertical line.

Ehrlich et al. (2019a)Kliesch and Mech (2019)Neuber et al. (2019)Mertes et al. (2019)Zanatta and Herber (2019)Chechin (2019)Dupuy et al. (2019)J. Hartmann et al. (2019)Savtchenko et al. (2004)

Table 1 Overview of the observational datasets used in this study. Data in rows 1–7 are accessible through the PANGAEA database, while the MODIS data are available through the NASA Worldview^* interface.

^* https://worldview.earthdata.nasa.gov (last access: 20 July 2022).

The eastern part of the target domain that is centered around dropsonde DS04 is selected as the area to be simulated (indicated by the green box in Fig. 1 and the green line in Fig. 2). This choice is motivated by the following considerations. First, marine mixed-phase convection did occur in this area. Second, the combination of warm water and cold air implies large surface latent and sensible heat fluxes, making the near-surface convection vigorous and potentially well coupled to the cloud layer. Such convective conditions often occur in this region, for a large part controlled by the wind direction. Third, convective clouds are well resolved in large-eddy simulations. A further advantage is that the cloud-bearing low-level air mass in this area was slow-moving and almost stagnant. This is evident from (i) the dropsonde profiles of u and v (Fig. 3) and (ii) the trajectory staying in the proximity to the dropsonde location for about 24 h (Fig. 1). This broadens the time span that the simulation results can be justifiably compared to relevant measurements. The P5 and P6 aircraft visited the area between waypoints C5 and C6 multiple times, which enhances the sample size. Finally, the decision to only simulate the eastern part of the target domain is mainly based on the need to avoid averaging large-scale forcings over a too wide area so that the simulated convection remains optimally representative of the local (convective) conditions surrounding DS04.

Figure 3 Dropsonde profiles from RF20 of (a) virtual potential temperature θ_v; (b) water vapor specific humidity q_v; (c) relative humidity RH; and the (d) zonal and (e) meridional wind speeds u and v, respectively. The DS04 sounding is shown in grey. The idealized initial (init) profiles are shown in orange.

The DS04 dropsonde data provide detailed insight into the thermodynamic structure of the marine boundary layer in the simulation domain. Figure 3a and b show that an inversion layer is present between 1.3–1.5 km height and can be recognized in most state variables. A well-mixed sub-cloud layer of about 350 m depth is capped by a relatively deep cloud layer, characterized by high relative humidity values and a conditionally unstable thermodynamic vertical structure. The inversion features a θ_v jump of about +3 K. In contrast, the jump in water vapor specific humidity q_v is relatively small, indicating the presence of significant water vapor above the inversion. There is a notable gap in wind measurements between 1.6–2.5 km, where samples have been removed after quality checks.

The observational data on hydrometeor occurrence and mass as collected by P5 and P6 are used to evaluate the LES experiments. While the MiRAC and AMALi data on board P5 provide information about cloud top heights, vertical structure, and liquid water path, the Nevzorov probe on board P6 provides in situ samples of cloud liquid water content. The Laboratoire de Météorologie Physique (LaMP) cloud imaging probe (CIP) provides data on both liquid and ice water content, using a Brown and Francis (1995) mass diameter relationship on non-spherical particles only. Finally, the high-frequency (100 Hz) turbulence measurements collected by the P6 nose boom (Hartmann et al., 2018) allow for calculating (co)variances of temperature and vertical velocity in the boundary layer, even for relatively short flight segments.

A unique aspect of RF20 is the operation of two ultra-high-sensitivity aerosol spectrometers (UHSASs) on board P6, both sampling submicron-size particles (Zanatta and Herber, 2019). These instruments measured the number size distribution of particles with diameters between 60 and 1000 nm by detecting scattered laser light, using the method described in detail by Cai et al. (2008). Figure 4 shows that a relatively wide range of aerosol concentrations was encountered. The data, as documented by Mertes et al. (2019), suggest a marked difference in aerosol loading below and above the cloud layer, with the lower values found below the clouds. The availability of in situ aerosol data greatly helps to constrain the LES, given the well-known sensitivity of mixed-phase clouds to this parameter. Accordingly, the UHSAS data form the basis for the initial CCN profile as adopted in the simulations described in Sect. 3.4. We hypothesize that the removal of aerosol by precipitation is the main reason for the lower values below the clouds; if this is the case, it should also show up in the simulations.

Figure 4 The probability density function (PDF) of aerosol concentrations (for aerosol particles of diameters between 80 and 1000 nm) as observed by the UHSAS-2 mass spectrometer on P6 during the southern racetrack section of RF20. The southern racetrack is shown in Fig. 1b, which was in this region between 15:40 and 17:10 UTC. The time resolution of the observation is 3 s. The red dashed line indicates the initial CCN concentration used in the simulations (Zanatta and Herber, 2019; Mertes et al., 2019).

3 Model configuration

3.1 DALES

The simulations in this study are carried out with the Dutch Atmospheric Large-Eddy Simulation model (DALES; Heus et al., 2010). DALES has been successfully applied to simulate observed turbulent/convective boundary layers and clouds in many climate regimes, including the tropics (Vilà-Guerau de Arellano et al., 2020), the subtropics (Van der Dussen et al., 2013; de Roode et al., 2016; Reilly et al., 2020), mid-latitudes (Neggers et al., 2012; Corbetta et al., 2015; Van Laar et al., 2019), and high latitudes (de Roode et al., 2019; Neggers et al., 2019; Egerer et al., 2021; Neggers, 2020a, b). The code of DALES is open source and maintained online at https://github.com/dalesteam/dales (last access: 8 November 2021). The governing equations, numerical aspects, and the various subgrid physics packages are described in detail by Heus et al. (2010), and accordingly only a brief summary is provided here. At the foundation of the model are the Ogura–Phillips anelastic equations for a set of prognostic variables including the three velocity components u, v, and w; total water specific humidity q_t; and liquid water potential temperature θ_l, as well as the number concentration and mass concentration of various hydrometeor species. The time integration makes use of a third-order Runge–Kutta scheme. Scalar advection is represented by the centered difference method, which is also applied for the three velocity components. The subgrid-scale transport of heat, moisture, and momentum is dependent on a prognostic turbulence kinetic energy (TKE) model. For longwave and shortwave radiation a multi-waveband transfer model is used in combination with a Monte Carlo approach (Pincus and Stevens, 2009). A new mixed-phase cloud microphysics scheme was recently implemented, as described in more detail in the next subsection.

3.2 Microphysics

The control version of DALES (Heus et al., 2010) includes a double-moment microphysics scheme for warm clouds featuring two hydrometeors, cloud water, and rain (Seifert and Beheng, 2001). To simulate Arctic mixed-phase clouds the mixed-phase extension described by Seifert and Beheng (2006a) (hereafter SB06) was recently implemented, adding a further three prognostic hydrometeors (cloud ice, snow, and graupel). This is a full two-moment implementation; the mass concentration and number concentration of each of five hydrometeors are thus prognostic variables. The first DALES results with this scheme, including an evaluation against RV Polarstern cloud measurements during the PASCAL campaign in 2017, are described by Neggers et al. (2019). In principle the implementation in DALES closely follows SB06. In this section some details of the implementation are described that are either (i) different from the SB06 description or (ii) particularly relevant for this study.

A key difference with SB06 is the prognostic treatment of the number concentration of cloud condensation nuclei (CCN). This first applies to activation of CCN in saturated grid cells. The CCN concentration is conserved during nucleation of cloud droplets, their condensational growth, and evaporation. The sedimentation of cloud droplets contributes together with the convection to the vertical transport of CCN. The self-collection of cloud droplets and precipitating processes act as sinks for CCN. For simplicity the collection of cloud droplets by ice particles and the freezing of cloud droplets are also treated as CCN sink terms. The glaciation of clouds does not cause CCN depletion because in the Wegener–Bergeron–Findeisen regime the vapor deposition on growing ice crystals evaporates liquid water but leaves the surrounding CCN unaffected (Schwarzenböck et al., 2001).

The primary ice production in SB06 accounts for ice nucleation, as well as freezing of cloud droplets and raindrops. Ice nucleation follows the parameterization proposed by Reisner et al. (1998), prescribing a constant number concentration of the available ice-nucleating particles (INPs) and ignoring any removal (see Appendix B). The freezing of liquid hydrometeors is described by the stochastic model proposed by Bigg (1953). The secondary ice production accounts for ice multiplication by the Hallett–Mossop process (Hallett and Mossop, 1974), occurring during the riming of ice hydrometeors in the temperature ranges between 265 and 270 K (Griggs and Choularton, 1986; Beheng, 1982). Other mechanisms of secondary ice production are not considered. Processes modifying the number and mass of ice hydrometeors include deposition, riming, aggregation of snow, self-collection of snow, partial conversion of snow and ice crystals to graupel, collection of snow by graupel, sublimation, melting, evaporation, and enhanced melting (i.e., melting due to collisions with liquid hydrometeors in temperatures above freezing point). The contributions to number and mass tendencies by these microphysical processes are calculated in the order established by Seifert and Beheng (2006b).

The majority of parameters in the DALES microphysics scheme follow the control setup of SB06, with the exception of the values of coefficients for shape and velocity of cloud ice. These were adjusted to the same values as adopted in the recent intercomparison study on a marine cold-air outbreak by de Roode et al. (2019) to better reflect conditions in Arctic low-level clouds. This decision is also motivated by the fact that both cases describe conditions in more or less the same region. The full setting of microphysical parameters adopted in this study is provided in Appendix B and Table B1. As will be shown in Sect. 4, the SB06 scheme carries sufficient complexity to satisfactorily reproduce the observed characteristics of the mixed-phase cloud layer.

3.3 Initialization, boundary conditions and composite forcing

The back trajectory calculated from the time and location of the DS04 dropsonde indicates that the air mass did not move within a degree of this location in the time period between 00:00 UTC and the dropsonde launch (see Fig. 1). This reflects the approximately stagnant wind conditions as also detected by the DS04 dropsonde (see Fig. 3). In addition, the large-scale conditions did not change much on this day. These conditions motivate adopting a time-composite case setup that reflects large-scale conditions in the region as averaged over the 12 h period leading up the DS04 launch.

Large-scale data from the Integrated Forecasting System (IFS) of the European Centre for Medium-range Weather Forecasts (ECMWF) are used to represent the impacts of larger-scale phenomena during the simulation. Following Van Laar et al. (2019), a combination of analyses (available every 12 h) and short-range forecasts (available every 3 h) is used, effectively yielding a four-dimensional dataset of the atmospheric state variables θ_l, q_t, u, and v at 3-hourly temporal resolution and 0.1×0.1^∘ spatial resolution. In this study these are calculated at 3-hourly points along the back trajectory, a method previously adopted by Neggers et al. (2019). The forcing profiles are time and height dependent. Horizontal advective forcings are represented as prescribed advective tendencies, calculated through horizontal averaging within a $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$-wide column around the location. The tendency due to large-scale subsidence relies on a prescribed profile of pressure velocity that acts on the evolving vertical structure in the LES. Forcings in the momentum equation include the Coriolis term and the pressure gradient term, in combination expressed as the departure of the model wind from the prescribed geostrophic wind. The latter is calculated from the pressure field. Given these time- and height-dependent forcing profiles at the trajectory points, time averaging is then applied over 12 h to obtain the composite forcing dataset used to drive the LES. These profiles are described in more detail in Appendix A) and are characterized by a persistent low-level subsidence, a weak advective cooling and moistening tendency above 1 km height, and negligible geostrophic forcing of the wind.

The DS04 dropsonde profiles are used as an initial state, amalgamated with the composite large-scale model data. Where available, the sonde data are averaged onto the vertical grid of the composite ECMWF forcings. At grid layers where no DS04 data are available, the composite model state itself is used instead. Figure 3 shows the initial profiles thus obtained, illustrating that the method successfully yields profiles that are continuous and do not include huge jumps that could result from mismatches between sonde and ECMWF data.

The surface boundary conditions include a prescribed skin temperature and humidity. The latter is calculated by assuming oceanic ice-free conditions so that the associated saturation specific humidity can be used. These skin values are then used to interactively calculate the surface fluxes of heat, moisture, and momentum, using prescribed roughness lengths for heat, moisture, and momentum. The calculation of the bulk drag and exchange coefficients relies on the stability functions that are native to the DALES code (Heus et al., 2010). A prescribed surface albedo is used to calculate the upward shortwave radiative flux at the lower boundary. The incoming shortwave radiative flux at the model ceiling depends on (i) seasonality and time of day and (ii) the composite large-scale state above the model ceiling. In doing so we follow the method adopted by Van Laar et al. (2019). This composite large-scale state above the model ceiling is also used to determine the downward longwave radiative flux.

All large-scale forcings are time constant and are applied in a horizontally homogeneous way, being the same in every LES grid column. In addition, horizontally periodic boundary conditions are applied. The resulting simulation can thus be interpreted as a statistical downscaling of the dropsonde profile, with the LES acting as a generator of the small-scale variability existing in the domain around it.

3.4 Experiment overview

This study makes use of one control simulation, designed to match the observed thermodynamic and cloudy state as closely as possible. This experiment is to serve as a benchmark simulation for planned subsequent studies (not covered in this paper). Table 2 summarizes the main characteristics of this experiment. The size of the simulated domain is considered wide and high enough to accommodate the typical width of convective structures observed in the Arctic (Müller et al., 1999). The duration of the experiment is 72 h, in order to provide the turbulent boundary layer with enough time to equilibrate and to cover three full diurnal cycles. Radiation is interactive with all five hydrometeors. Above the turbulent domain the composite ECMWF profile is used in the calculation of the downward fluxes at the model ceiling. The solar inclination is time and latitude dependent, introducing a diurnal cycle in the radiation.

Table 2 Summary of defining characteristics of the LES experiments for RF20.

A sponge layer is applied in the top 1 km of the computational domain to prevent reflection of waves off the rigid top boundary (Shepherd et al., 1996; Heus et al., 2010). Additionally, continuous nudging towards the initial profile is applied above the boundary layer inversion to prevent excessive model drift in this height range. To this purpose a Newtonian relaxation term (Neggers et al., 2012) is included in the prognostic equations for u, v, θ_l, and q_t, adopting a timescale of 3 h. Inversion height is calculated interactively, defined as the level at which the vertical gradient in θ_l is strongest. A 300 m deep transition zone is included above the inversion across which the nudging intensity increases linearly with height. In this configuration, the resolved turbulence and convection below the inversion remains unaffected by the free-tropospheric nudging and can freely equilibrate in response to the initial and boundary conditions and the prescribed time-constant forcings. Such nudging has successfully been applied in previous LES studies in which equilibration played an important role (Sandu et al., 2009), motivating adopting the same technique here.

The observed statistical distribution of aerosol concentration as shown in Fig. 4 informs the initial CCN profile adopted in the control experiment, chosen here to be constant with height at 100 cm⁻³. This choice is motivated by the following arguments. Firstly, it is safe to assume that only a fraction of the observed aerosol can act as CCN, here assumed to be approximately 50 %–90 %. Secondly, CCN are treated prognostically in the model and can evolve freely during the simulations, and no external sources of CCN are considered for simplicity. With the convection and clouds gradually removing aerosol below the inversion (Bigg and Leck, 2001), a choice of 100 cm⁻³ should after some time result in a vertical structure resembling the observed one as shown in Fig. 4, with concentrations below the clouds being about half of the values above. The sensitivity to CCN is tested by means of two additional experiments. One experiment adopts 10 cm⁻³, representing pristine air as often found in the high Arctic. The other adopts 1000 cm⁻³, representing polluted air from continental origins. Both extremes are frequently observed in Arctic air masses (Bigg et al., 1996; Stohl et al., 2006; Garrett et al., 2010).

The initial cloud droplet number concentration in supersaturated areas is set accordingly: the initial mean droplet size must be higher than the advised threshold ${\stackrel{\mathrm{‾}}{x}}_{\min }=\mathrm{4.2}\times {\mathrm{10}}^{-\mathrm{12}}$ g (Seifert and Beheng, 2006a), and at most $\mathrm{1}/\mathrm{2}$ of the initial CCN number concentration is activated. The motivation for this initial concentration is straightforward: the simulation can start with neither too large nor too small droplets and spins up without encountering a lack of CCN. The number concentration of INPs is assumed time constant and prescribed at 1000 m⁻³, therein following Reisner et al. (1998) and M. Hartmann et al. (2019).

4 Results I: evaluation

The main goal of the evaluation of the control experiment is to test to what extent the bulk characteristics of the observed boundary layer and clouds are reproduced. This approach is in line with previous LES studies of Arctic mixed-phase clouds (e.g., Ovchinnikov et al., 2014; Stevens et al., 2018; Neggers et al., 2019). Bulk statistics here include the mean vertical structure as well as various heights, such as the cloud boundaries and inversion. The availability of in situ cloud and turbulence measurements during RF20 allows us to go one step further and also evaluate profiles of liquid and frozen hydrometeor mass as well as turbulent (co)variances.

4.1 Time evolution

Figure 5a and b document the time development of the domain-averaged cloud structure and phase during the control simulation of the composite RF20 case. After a spinup period of about 12 h the boundary layer more or less equilibrates, staying close to this state for the remainder of the simulation. The cloud layer is in the mixed phase, with liquid cloud coexisting with ice cloud. Both specific mass concentrations reach a maximum immediately below the inversion. A weak diurnal cycle in cloud mass, cloud thickness, and inversion height can be distinguished, superimposed on the long-term equilibrium (better visible in time series data shown later in Fig. 13). Such diurnal signals are well known from warm marine stratocumulus and are driven by daytime absorption of shortwave radiation alternating with nighttime cloud top cooling (Rozendaal et al., 1995; Wood et al., 2002; Duynkerke et al., 2004). The presence of this signal here indicates that in early June the solar radiation is already strong enough at these latitudes to drive a boundary layer response.

Figure 5 Time–height contour plots of domain-averaged variables during the control experiment. Specific mass of (a) cloud liquid water q_l and (b) cloud ice water q_i. The averaging time is 30 min. The dashed and dotted lines reflect the lowest base and highest top of liquid clouds in the domain, respectively.

Figure 6a and b show the total (resolved + subgrid) fluxes of virtual potential temperature and total specific humidity. Significant boundary-layer-deep transport is present, reflecting a high degree of coupling between the cloud layer and the surface. This aspect makes this case distinctly different from turbulent mixed-phase clouds over homogeneous sea ice, which are often fully decoupled (e.g., Solomon et al., 2014). The evolution of the humidity flux shows that it takes about 12 h for the surface-driven turbulence to properly spin up after initialization. Once this has occurred the transport is continuous, with the local minimum in the θ_v flux near the liquid cloud base indicating the presence of a shallow transition layer (see Fig. 6a). Such stable layers with CIN (convective inhibition) are a well-known feature of cumulus-capped boundary layers (Albrecht et al., 1979) and decoupled stratocumulus (Nicholls, 1984). The transport intensity is intermittent at times, which could reflect the impact of subsampling of convective events due to a still limited domain size.

Figure 6 Same as Fig. 5 but now showing (a) the turbulent θ_v flux and (b) turbulent q_t flux.

4.2 Vertical structure

The simulated vertical profiles of three state variables at three time points are compared to the dropsonde sounding data in Fig. 7. In general the simulated profiles agree well with the observations concerning the vertical structure of the boundary layer, with a relatively well-mixed layer up to about 500 m capped by a cloud layer that features a conditionally unstable thermodynamic lapse rate. This cloud layer with high relative humidity extends up to about 1.4 km. The inversion layer of a few hundred meters deep is also reproduced well, featuring a θ_v jump of about +3 K and a negligible q_v jump. During the 48 h simulation this vertical structure does not change much, with the boundary layer only deepening by about ∼150 m.

Figure 7 Domain-averaged vertical profiles during the control experiment. (a) Virtual potential temperature θ_v. (b) Water vapor specific humidity q_v. (c) Relative humidity RH. The DS04 dropsonde sounding is shown in grey. Three subsequent time points are shown.

Figure 8 shows domain-averaged vertical profiles of all five hydrometeor species of the DALES microphysics scheme at t=24 h, grouped as suspended and falling species. This time point was chosen based on the good agreement of the inversion height with the observations (see Fig. 7). Cloud liquid water q_l has a distinct mode near the inversion but still has significant values below, reflecting the presence of rising convective updrafts. Cloud ice q_i also peaks near the inversion but disappears a few hundred meters above the liquid cloud base. These amplitudes in q_l and q_i compare well to previous studies of convective mixed-phase clouds over open water in the Arctic (Klein et al., 2009; Morrison et al., 2009). Snow q_s and graupel q_g masses peak in the middle of the cloud layer and turn into rain q_r below the freezing-level height at z∼400 m. Significant precipitation mass is lost on the way down due to sublimation and evaporation (at the altitudes below the freezing level). Rain shafts reaching the surface are consistent with the observations; they are visible in the radar measurements shown in Fig. 2.

Figure 8 The specific masses of all five hydrometeor species in the control simulation at t=24 h. (a) Suspended species, including liquid water q_l and ice water q_i. (b) Precipitating species, including rain q_r, snow q_s, and graupel q_g. Liquid species are shown in blue; frozen species are shown in pink. The vertical range of AMALi cloud top heights sampled within the target domain as shown in Fig. 2 is shaded grey.

4.3 Cloud boundaries and hydrometeors

Figure 8 also indicates the range of cloud top heights as sampled by the AMALi instrument on board P5 (data shown in Fig. 2). This range only reflects AMALi data sampled within the modeling domain, which was visited three times in short succession. Because the cloud deck was not homogeneous and contained numerous gaps, the data include sufficient hits at lower heights to yield a representative estimate of cloud layer depth. The model liquid water fits well between the AMALi cloud layer boundaries in the control experiment at this time point (t=24 h), consistent with the good agreement on inversion height as detected by the dropsonde.

Figure 9 Simulated and observed hydrometeor masses. Panels (a) and (b) show liquid water content (LWC) sampled by the Nevzorov and CIP probes, respectively, while panel (c) shows the ice water content (IWC) measured by the CIP probe. Only samples during the southern racetrack segment of RF20 are included. The hourly mean LES results at t=24 h are shown in color and represent conditional averages over the area fraction covered by the associated hydrometeors. The thick colored line represents the sum of all liquid-phase hydrometeors (blue, q_l+q_r) and ice-phase hydrometeors (pink, $q_{\mathrm{i}}+q_{\mathrm{s}}+q_{\mathrm{g}}$) in the model. The thin dashed line represents only the suspended species (q_l and q_i, respectively). The grey box–whisker plots represent the P6 measurements in the racetrack section, with the 5th–95th percentile range (black line) and the interquartile range (grey box) indicated. The Nevzorov data are analyzed on five separate level flight legs, while the CIP data are vertically gridded at 250 m spacing.

Figure 9 evaluates the hydrometeor mass in the LES against the in-cloud measurements by the Nevzorov and CIP probes on board P6. The hourly mean LES results at t=24 h are shown in color, while the observational data are shown as box–whisker plots. A complication with hydrometeor evaluations of this kind is that their definition in the microphysics scheme might not necessarily match that of the observation system. For example, for ice hydrometeors it is hard to separate between suspended (cloud) and precipitating particles. To avoid such problems the comparison is made as straightforward as possible, by only making comparisons based on phase. This means that in each panel the simulated and observed values include all hydrometeor species in the phase of interest. For the LES this means $(q_{\mathrm{l}}+q_{\mathrm{r}})/\mathit{\rho }$ for liquid and $(q_{\mathrm{i}}+q_{\mathrm{s}}+q_{\mathrm{g}})/\mathit{\rho }$ for ice (shown as thick lines in Fig. 9). For the CIP probe the full size distribution of ice particles can thus be included in the data. Note that the LES microphysics scheme does include secondary ice production, allowing for the formation of large and heavy ice hydrometeors in the simulation (Sullivan et al., 2017; Georgakaki et al., 2022). For reference, the mass of only the suspended model hydrometeors is also indicated (dashed line).

Figure 9a and b compares the liquid-phase hydrometeor mass of the LES to the Nevzorov and CIP measurements, respectively. Note that the vertical gridding of the data are slightly different for each instruments; for the Nevzorov probe the box–whisker plots represent data on five separate straight flight legs, while for the CIP probe a vertical gridding of 250 m is used. The model data are conditionally averaged over the area where the hydrometeors occur; this way, the contribution by hydrometeor-free air to the average is excluded, and the comparison with the in-cloud observational data are fair. While slight differences exist in liquid water content (LWC) between the two probes, in general they agree on the magnitude and vertical structure. The observed mean LWC peaks in the middle of the cloud layer at ∼800 m height and again near the inversion at ∼1400 m. The LES always sits within the observed 5th–95th percentile range in the cloud layer and is reasonably close to the interquartile range. It also reproduces the vertical structure, which is encouraging. Some LWC is detected below the AMALi cloud base (∼400 m), which is probably rain. Including rain q_r in the LES data shows that the model is consistent with this feature.

Figure 9c then evaluates the mass of all frozen hydrometeors against the CIP data. The data show that large spread exists but that on average the IWC mass is considerably smaller in magnitude compared to the observed LWC. Also note that some ice is observed below the cloud base. The model is again situated within the 5th–95th percentiles. This time, a considerable difference exists between the suspended and falling hydrometeors in the LES; suspended cloud ice q_i is only present near the inversion (dashed line), while snow q_s and graupel q_g contribute most of the frozen hydrometeor mass at lower heights. Including those falling species is needed to explain the observations at lower heights.

We conclude from this analysis that the model is representative of the observed clouds, both for the liquid and ice phase. It should be noted in this respect that the observed sample size of hydrometeors is still limited; only relatively few clouds were sampled by P6 on the racetrack, also during a limited period. This introduces some uncertainty, preventing us from drawing any conclusions beyond the bulk (mean) state. Any higher-order evaluations, for example concerning spatial variability or particle size distributions, require much more substantial datasets and therefore go beyond the scope of this study.

Figure 10 Simulated and observed turbulent (co)variances, including (a) the vertical velocity variance $\stackrel{\mathrm{‾}}{{w^{\prime }}^{\mathrm{2}}}$, (b) the temperature variance $\stackrel{\mathrm{‾}}{{T^{\prime }}^{\mathrm{2}}}$, and (c) the associated temperature flux $\stackrel{\mathrm{‾}}{w^{\prime }{T}^{\prime }}$. Similar to Fig. 9 the box–whisker plots again indicate P6 measurements, this time showing nose boom data during the five level flight sections on the southern racetrack. The mean LES profile at t=24 h is shown in red.

4.4 Turbulence

The measurements at 100 Hz of temperature and vertical velocity made by the sensors in the P6 nose boom during RF20 (M. Hartmann et al., 2019) allow for calculating variances of both variables, as well as the turbulent heat flux. On the RF20 southern racetrack five segments at constant heights within the lowest 2 km were included for this purpose, adopting the well-known method as explored in previous classic studies (Nicholls and Lemone, 1980). Although this area is located slightly to the west of the DS04 target area, the near-surface turbulence is expected to be similar in both. A running average of 10 min is calculated, based on which time series of the (co)variances can be estimated.

Figure 10 shows the distribution of (co)variances on each of the five level flight legs. Observational data are again shown in the box–whisker style, similar to Fig. 9. LES data again reflect hourly averages around t=24 h. Below the cloud base (∼450 m) the LES profiles exhibit some well-known features of turbulent mixed layers, including a convex $\stackrel{\mathrm{‾}}{{w^{\prime }}^{\mathrm{2}}}$ structure, a concave $\stackrel{\mathrm{‾}}{{T^{\prime }}^{\mathrm{2}}}$ structure, and a linear decrease in the heat flux $\stackrel{\mathrm{‾}}{w^{\prime }{T}^{\prime }}$ towards slightly negative values near the layer top. The latter is consistent with the local minimum in the buoyancy flux seen in Fig. 5d. These vertical structures are not clearly visible in the observations; however, a large spread exists at each flight leg, and the model is almost always situated within the observed 5th–95th percentile range. The exception is $\stackrel{\mathrm{‾}}{{w^{\prime }}^{\mathrm{2}}}$ at the lowest and highest flight legs, which we speculate could be caused by the strong variability in w in the close vicinity of convective cells. In the middle of the cloud layer (∼900 m) the model data show a second maximum, reflecting the impact of latent heat release on turbulence. Near the inversion the temperature variance peaks as expected, due to the close vicinity of a strong temperature gradient. This feature in the model is also supported by the measurements.

The general outcome of this evaluation is that the LES turbulence profiles are mostly situated within the observed 5th–95th percentile range, with a few small exceptions. The data suggest that the agreement is best for the temperature variance and the vertical heat flux, with the LES sitting reasonably close to the interquartile range. This supports the conclusion that the amplitude and vertical structure of both the intensity and transport of the convection and turbulence in this case are reasonably well captured by the model.

Figure 11 Time–height contour plots of liquid cloud fraction (cfrac) during the three CCN sensitivity experiments. (a) Pristine high-Arctic conditions (N=10 cm⁻³). (b) Control conditions observed during ACLOUD RF20 (N=100 cm⁻³). (c) Polluted continental conditions (N=1000 cm⁻³). Liquid cloud fraction is calculated at each vertical layer as the fraction of grid points where the cloud liquid water content is above the threshold (set to 0.01 g m⁻² for cloud liquid water).

5 Results II: aerosol impacts

With a satisfactory agreement between the control run and the measurements established, the next step is to assess the sensitivity of the results to the aerosol levels in the simulated air mass. A few recent LES intercomparison studies on Arctic mixed-phase clouds have investigated this dependence. Stevens et al. (2018) investigated three LES codes for a tenuous mixed-phase stratocumulus case observed over sea ice, reporting that the lack of CCN can seriously limit the liquid water path (LWP). Higher CCN concentration levels in a marine cold-air outbreak case were studied by de Roode et al. (2019), finding that this leads to increased LWP. The RF20 case studied here is similar to the CAO case but also differs considerably, in particular in the much weaker convection and the stagnancy of the air mass. Our goal is to investigate and understand the impact of CCN also for these conditions. Specific focus lies on the efficiency of radiatively driven entrainment, the role played by the ice phase, and the impacts on the energy budgets of the boundary layer and the surface. A much wider CCN concentration range is covered compared to the previous LES studies. This allows for interpretation of the consequences of CCN impacts for air mass transformations and air–ocean interactions in this region of the Arctic.

5.1 CCN

The two additional experiments with N={10, 1000} cm⁻³ are designed to reflect observed extreme CCN conditions in the area of interest. The lower value represents pristine conditions typical of the high Arctic, while the higher value represents polluted continental conditions as sometimes encountered in warm and moist air intrusions from the mid-latitudes (Bigg et al., 1996). The N=100 cm⁻³ value of the control experiment is based on in situ P6 measurements in the air mass during ACLOUD RF20. Accordingly, together these values span a realistic CCN range.

Figure 12 Profiles of (a) liquid water potential temperature θ_l, (b) total water specific humidity q_t, and (c) the CCN concentration N for the three CCN experiments. Data are time-averaged over the last 48 h. Colors represent the initial CCN concentration (cm⁻³). The modes of the measured concentrations below (×) and above (+) the cloud layer (shown in Fig. 4) are also indicated, for reference.

Figure 11 shows time–height contour plots of the area fraction of liquid cloud water in three CCN sensitivity runs. The depth of the boundary layer shows a strong response, deepening from 1.2 km to above 2 km for the continental setup. The capping liquid cloud layer situated immediately below the inversion also thickens considerably. Note that only the control N=100 cm⁻³ experiment is close to the observed cloud top height. During the first 48 h of the continental experiment the deepening is so aggressive that the stratiform cloud layer shows signs of decoupling from the surface-driven convection below, probably because of large top entrainment rates. This decoupling is indicated by a steepening of the virtual potential profile and a gap in the cloud fraction forming above 1100 m, nearly separating the cloud layer into two distinct layers in Fig. 5c (32–48 h interval).

Figure 12 documents the impact of CCN on the thermodynamic structure and depth of the convective boundary layer. For larger N values, the convective mixing reaches deeper into the air mass in which it is embedded, lifting the thermal inversion. An important side effect is that it also makes the thermal inversion stronger, as expressed by the increased jump in θ_l across the inversion layer. The boundary layer in the continental CCN experiment has also warmed and moistened significantly compared to the other two runs. The vertical structure of N reveals its time evolution, with concentrations below the inversion gradually decreasing due to removal by precipitation. During the simulated period this decrease is only limited, with CCN levels in the convective layer still retaining substantial values. This vertical structure in the model is also consistent with the P6 observations, with the N=100 cm⁻³ experiment sitting closest to the observed amplitude above and below the cloud layer.

Figure 13 compares the evolutions of vertically integrated or projected cloud properties during the three CCN sensitivity runs. The liquid cloud cover a_l is in the model defined as the fraction of the domain where the vertically integrated cloud liquid water exceeds the minimum threshold of 0.01 g m⁻². The liquid cloud cover goes from about four octas for pristine conditions to persistent full cover for continental conditions. The liquid water path (LWP) also increases with CCN, temporarily reaching very high levels during the spinup of the high-CCN run, before settling at a lower value. The ice water path (IWP) increases more or less monotonically. None of the experiments have fully equilibrated at the end of the simulation, with the N=1000 cm⁻³ experiment seeing the largest drift. A weak diurnal cycle is visible in all three variables. Interestingly, before the appearance of cloud ice the three cases evolve almost identically; it is only afterwards that significant differences develop. Although the ice initiation is not dependent on the cloud droplet size, the large droplets nevertheless contribute to the development of the ice phase. The reasons are the following: firstly, larger cloud droplets are more likely to undergo heterogeneous freezing; secondly, larger cloud droplets are more likely to collide with ice particles and freeze. The ice particles then join more mass, leading to increased fall velocity and thus further collisions. Additionally, when the collision occurs, the freezing large droplet produces more splinters (Hallett and Mossop, 1974). On the other hand in the simulation with high CCN concentrations, the cloud water is distributed over higher-number cloud droplets, leading to weaker riming tendencies and thus slower development of the ice phase (see Fig. 15b). This suggests that ice formation plays a key role in how CCN impacts the boundary layer clouds.

Figure 13 Time series of domain-averaged microphysics properties during the three CCN sensitivity experiments. (a) Liquid cloud cover a_l. (b) Cloud liquid water path LWP. (c) Cloud ice water path IWP.

Figure 14 then summarizes the impact of CCN on a selection of bulk properties of the convective boundary layer, presented as a bar plot of time averages that cover the last 48 h of the simulation. This is done to effectively remove any diurnal cycle in these signals. The changes shown in the top four panels indicate that the response of the boundary layer clouds to a CCN increase as detected in this case is in principle similar in sign compared to the cold-air outbreak case studied by de Roode et al. (2019). Boundary layer depth z_t is here calculated as the height at which ${\partial }_z\stackrel{\mathrm{‾}}{{\mathit{\theta }}_{\mathrm{l}}}$ is maximum. Depth z_t increases along with the liquid and frozen cloud amount LWP and IWP, while the surface precipitation P decreases, consistent with the second cloud–aerosol indirect effect.

Figure 14 Bar plots of diurnally averaged bulk properties of the convective boundary layer, for all three CCN experiments as shown in Fig. 13. (a) Boundary layer depth z_t. (b) LWP. (c) IWP. (d) Surface precipitation rate P. (e) Net longwave radiative flux divergence across the liquid cloud layer ΔF_rad. (f) Top entrainment rate ϵ_t. (g) Large-scale vertical velocity at boundary layer top W_t. (h) Jump in liquid water potential temperature Δθ_l across the inversion. The time averaging covers the last 48 h of the simulations.

Figure 15 Bar plots of the effects of cloud processes. (a) The comparison of the shortwave downward flux estimated from the Twomey formula (shown as bars with the pattern) and the shortwave downward flux from the from the simulation output (shown as full colored bars); (b) vertically integrated tendencies in (total) ice water mass budget within the cloud layer: primary ice nucleation and freezing of cloud droplets (Pr_i, significantly smaller than other tendencies), the combined effect of deposition and sublimation (Dep + Sub), riming of cloud droplets (Rim), and sedimentation (Sed). The time averaging covers the last 48 h of the simulations.

Another indirect effect of the aerosols in the Arctic is the change in the shortwave downward radiative flux due to changes in the optical properties of the clouds (Markowicz et al., 2021; Im et al., 2021). We estimate the Twomey effect using the simple formula (Twomey, 1977) for the optical depth τ:

$$\begin{array}{}\text{(1)}& \mathit{\tau }=\mathrm{2}\mathit{\pi }N{\stackrel{\mathrm{‾}}{r}}^{\mathrm{2}}h,\end{array}$$

where N is the number concentration of the cloud droplet per unit of volume, $\stackrel{\mathrm{‾}}{r}$ is the mean radius of cloud droplets, and h is the depth of the cloud. Based on the τ calculated from the mean values over the cloud layer, we estimate the transmitted shortwave downward radiative flux and compare it with the shortwave downward flux from the simulation output (see Fig. 15a). There is a reasonably good agreement for the pristine conditions (N=10 cm⁻³) but a significant underestimation for the control case. Finally for the continental case, the estimated Twomey effect differs from the simulation by more than an order of magnitude. However, this outcome is not surprising in the Arctic (Yang and Liu, 2022) given the non-uniform structure of the mixed-phase clouds with regions of lower cloud fraction in the middle and lower portion of the cloud layer (Fig. 11c).

The lower four panels of Fig. 14 provide more insight into the impact of CCN on the turbulent mixing at the inversion. As one expects, the net longwave radiative flux divergence across the liquid cloud layer ΔF_rad increases with LWP. It seems to saturate towards very large CCN, reflecting shifts in both cloud top cooling and cloud base warming (Stevens et al., 2005). As is well known, the cloud top cooling generates turbulence, which in turn drives entrainment of overlying dry and warm air into the boundary layer. The associated top entrainment rate ϵ_t at z_t is here calculated as the residual from the boundary layer mass budget:

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial z_{\mathrm{t}}}{\partial z}}={\mathit{ϵ}}_{\mathrm{t}}+W_{\mathrm{t}},\end{array}$$

with W_t being the large-scale vertical velocity at z_t. Figure 14f confirms that the entrainment rate increases in a similar tread as ΔF_rad. Note that the prescribed subsidence increases with height in this case (see Fig. A1). As a result, the CCN-richer boundary layer also equilibrates at a higher z_t. The slight inequality of ϵ_t and −W_t again reflects that the equilibration has not yet completed at the end of the simulated period.

Figure 16 Scatter plots of terms in the radiative entrainment efficiency α, as defined by Eq. (3). The net longwave flux difference across the cloud layer ΔF_rad is plotted against (a) the top entrainment rate ϵ_t and (b) the entrainment flux ϵ_tΔθ_l in energy units. Diagnostics at 2 h⁻¹ frequency are shown as small semi-transparent dots, while the 48 h mean is shown as an opaque big dot. In panel (a) the dotted line represents the linear fit through the means, while in panel (b) it represents the line at which α=1.

Figure 17 Same as Fig. 14 but now showing various surface (sfc) energy fluxes. (a) All four components of the net radiative energy flux Q^*, separated between longwave (LW) and shortwave (SW) bands as well as upward (u) and downward (d) directions. (b) All terms in the surface energy budget (SEB), including net downward radiation $-Q^{*}$, upward sensible heat H, upward latent heat E, and downward ground heat −G. The values are horizontally averaged over the whole surface; the time averaging then covers the last 48 h of the simulations.

Entrainment warming and cloud top cooling counteract each other in the heat budget of the boundary layer. Which one wins is expressed by the radiative entrainment efficiency α, defined by Stevens et al. (2005) as

$$\begin{array}{}\text{(3)}& \mathit{\alpha }={\displaystyle \frac{{\mathit{ϵ}}_{\mathrm{t}}\mathrm{\Delta }{\mathit{\theta }}_{\mathrm{l}}}{\mathrm{\Delta }{F}_{\mathrm{rad}}}},\end{array}$$

where Δθ_l is the jump in liquid water potential temperature across the thermal inversion, here defined as the layer between the first maximum and minimum in the second derivative of $\stackrel{\mathrm{‾}}{{\mathit{\theta }}_{\mathrm{l}}}$ below and above z_t, respectively. Figure 16 shows scatter plots of various components in Eq. (3). Although considerable scatter exists at small timescales, a well-defined linear relation exists between the time-averaged ϵ_t and ΔF_rad, confirming that cloud top cooling is driving the entrainment. Figure 16b then shows that cloud top cooling is almost always stronger than the entrainment heating, representing entrainment inefficiency (α<1). However, the efficiency improves towards higher CCN concentrations, approaching α=1 for continental conditions. At this point, entrainment warming is almost fully compensating the cloud top cooling. Comparing both panels indicates that this has to be caused by a non-linear increase in Δθ_l, a measure of thermal inversion strength (shown in Fig. 14h). The strengthening can be explained as a consequence of a deep mixed layer growing into a weakly stable overlying layer. Apparently, the boundary layer responds to a CCN concentration increase by both strengthening the inversion and boosting the warming efficiency of radiatively driven entrainment. This is an interesting outcome in the context of Arctic amplification, given that both processes play a role in the lapse rate feedback.

The impact on the surface radiative fluxes is considered in Fig. 17a. While the downward longwave flux sees a relatively small increase, the downward shortwave flux strongly reduces by about 200 W m⁻², reflecting that most of the solar radiation is now reflected at the liquid cloud top. Although this is less that what would be expected in the Twomey effect (as shown in Fig. 15a), it is nevertheless a significant difference. As a result, the downward net radiative flux at the surface $-Q^{*}$ (shown in Fig. 17b) reduces by about 130 W m⁻². This energy flux is the main driver of the surface energy budget (SEB), which, following Stull (1988), can be defined as

$$\begin{array}{}\text{(4)}& -Q^{*}=H+E-G,\end{array}$$

where H is the upward sensible heat flux, E is the upward latent heat flux, and −G is the downward heat flux into the ocean. As is often the case with the SEB over oceans for daytime weakly unstable conditions, H and E are typically much smaller than the net radiation and ocean heat flux. Moving from pristine to continental CCN values reduces both H and E by a small amount. But the reduction in both $-Q^{*}$ and −G is much larger, with the one more or less compensating the other. In other words, the flow of energy through the SEB is diminished, resulting in a reduced energy flux into the ocean.

6 Discussion

With the main results thus described in detail, we now proceed with a more general discussion of their implications, in various contexts.

6.1 Constraining LES with observations in the Arctic

The results show that the control setup yields a simulation that equilibrates close to the observed thermodynamic, turbulent, and cloudy state of the marine boundary layer as probed during ACLOUD RF20. Part of this control setup is the independently sampled CCN value of N=100 cm⁻³ by P6 during RF20 (see Fig. 4). Given the large sensitivity of the results to CCN as found in this study, the good agreement of the control experiment with the RF20 data is not trivial. It implies that the availability of such in situ aerosol measurements is crucial for the successful configuration of subsequent LES experiments that are to realistically reflect observed conditions. This finding is thus a recommendation for future field campaigns in the Arctic.

6.2 Air mass transformation

Figure 12 further illustrates the profound impact of CCN on the boundary layer structure and depth in this cloud regime over open water. These findings indicate that CCN can act as a catalyst for the convective mixing into the air mass. For higher CCN concentrations, the convective mixing reaches deeper into the air mass in which it is embedded, lifting the thermal inversion and making it stronger. In particular radiatively driven turbulent entrainment becomes more efficient compared to surface-driven convection. This CCN impact thus plays a role in how air masses transform as they travel over open water near the sea ice edge. In case it travels onward into the high Arctic, its properties at lower levels will be more diluted as a result, affecting its net impact on the Arctic climate system (Pithan et al., 2018). Numerical simulations covering a much larger domain in the Arctic but still resolving entrainment to a reasonable degree, as well accounting for the impact of CCN on its efficiency, could provide further insight.

6.3 Climate impacts

It is interesting to interpret the impact of CCN on the entrainment efficiency as found in this study in the context of Arctic amplification, in particular the lapse rate feedback (LRF). This would mainly work through the low-level heat budget. As defined by Eq. (3), the entrainment efficiency reflects the ratio of two terms in this budget, namely entrainment warming and longwave radiative cooling. But the full heat budget contains a few more terms, which might also play a role. To find out, we now consider the budget for the layer of air below the base of the boundary layer thermal inversion, here taken as the base height of the layer over which the jump Δθ_l was calculated (shown in Fig. 14h). The vertically integrated bulk budget of sensible heat can then be written as

$$\begin{array}{}\text{(5)}& c_{\mathrm{p}}h\langle \mathit{\rho }\dot{\mathit{\theta }}\rangle =c_{\mathrm{p}}h〈\mathit{\rho }{\dot{\mathit{\theta }}}_{\mathrm{p}}+\mathit{\rho }{\dot{\mathit{\theta }}}_{\mathrm{f}}〉,\end{array}$$

where h is layer depth. A dot indicates a tendency, and the brackets indicate a vertical average over depth h. A distinction is made between the contributions by physical processes and large-scale forcings, indicated by subscripts “p” and “f”, respectively. The physics term can be further expanded as

$$\begin{array}{}\text{(6)}& c_{\mathrm{p}}h〈\mathit{\rho }{\dot{\mathit{\theta }}}_{\mathrm{p}}〉=H+P+\mathrm{\Delta }{F}_{\mathrm{LW}}+\mathrm{\Delta }{F}_{\mathrm{SW}}+F_{\mathit{ϵ}},\end{array}$$

where ΔF_LW and ΔF_SW represent the difference in longwave and shortwave flux across the layer and F_ϵ stands for the entrainment flux at the layer top. Note that by considering only the layer below the base of the thermal inversion, the reduction in h (or loss of mass) due to subsidence within the inversion layer is excluded.

Figure 18 shows all budget terms as averaged over the last 48 h. By convention, a term is positive when it contributes heat to the layer. Longwave cooling ΔF_LW is always the main sink, counteracted by surface fluxes H and P, entrainment warming F_ϵ, and shortwave absorption ΔF_SW. These physics approximately balance the sink term due to layer-internal large-scale forcing (including prescribed horizontal advection and subsidence). Adding CCN significantly shifts this balance, with surface fluxes becoming less important and radiative and entrainment terms becoming more dominant. In other words, the turbulence shifts from surface driven to radiatively driven. The entrainment flux contributes most to the enhanced net heating tendency, which implies that the boosted entrainment efficiency as previously discussed in Sect. 5.1 drives the net response of the boundary layer energy budget.

Figure 18 Analysis of the bulk (BL) heat budget of the boundary layer below the base of the capping inversion. (a) Bar plot of all physics terms in the sensible heat budget, plotted as energy fluxes. (b) The net contribution by all physics terms as shown in panel (a) and the large-scale forcing $h{\dot{\mathit{\theta }}}_{\mathrm{LS}}$.

What these results imply is that the impact of CCN on the low-level energy budget through the entrainment efficiency is substantial and can not be ignored for understanding the LRF. In this respect it is worthwhile to remember that our simulations reflect clouds over open water. Recent studies have emphasized the importance of considering the partitioning between marine and ice-covered areas in understanding the LRF (Jenkins and Dai, 2021). In newly formed open-water areas where sea ice has disappeared, one expects the lower-level stability to decrease. However, our results suggest that this does not automatically imply a weakened inversion, in the case that the melt is accompanied by a simultaneous increase in CCN (Browse et al., 2014).

6.4 Atmosphere–ocean interactions

We find that for high CCN levels in the air mass the SEB is greatly reduced, featuring a much smaller flux of energy into the ocean. This has a cooling effect on the oceanic mixed layer, which in the long term might percolate to greater depths where the main ocean circulation takes place. Note that in our simulations the ocean skin temperature was prescribed such that G is calculated as a residual of the SEB. In other words, the changes in net radiation do not feed back into the sensible and latent heat fluxes, which only depend on the lower-atmospheric state. Additional LES runs including a simple but interactive oceanic bulk mixed layer would allow this two-way air–ocean interaction to take place. This could give insight into the typical timescale of the adjustment process, the sign and magnitude of the air–ocean feedback, and the amount of energy lost to deeper parts of the oceanic mixed layer. Such simulations are considered a future research effort. Another factor to be considered when interpreting the impact of the large radiative flux differences on the Arctic climate system is that the actual frequency of occurrence of this cloud regime in the area of interest is not yet considered in this study.

6.5 Further sensitivities

Some aspects of the experimental design can affect the behavior of the simulated turbulence and clouds. These include the size of the turbulent domain, the spatial and temporal discretization, and the way the larger-scale forcing is configured. For example, one could adopt heterogeneous boundary forcing in a nested setup, which allows for representing advection of mesoscale features into the domain that are now ignored. While this might aid realism, it would make the simulation more complex and perhaps harder to interpret in terms of the response of convection to CCN. Another key model component is the microphysics scheme, for which one also expects sensitivity. In a follow-up study the authors explore some of these potential sensitivities for this ACLOUD case.

7 Conclusions and outlook

Observational data collected by the P5 and P6 polar research aircraft in a relatively stagnant air mass over the Fram Strait during ACLOUD RF20 were used to test the skill of LES in reproducing observed key characteristics of mixed-phase convection over open water near the sea ice edge. A unique aspect of this campaign is the availability of in situ aerosol measurements directly above and below the cloud deck. These allow for realistic initialization of CCN in the LES and create opportunities for gaining insight into how aerosol content in the air mass affects the mixed-phase convection. Our main conclusions can be summarized as follows:

1. The observed thermodynamic, kinematic, turbulent, and cloud structure are well reproduced by the model when initialized with the observed CCN concentration.

2. Changing the CCN concentration in the air mass from pristine to continental values substantially alters the cloud amount, the boundary layer depth, the thermal inversion strength, and the amplitude of the surface energy budget.

3. The efficiency of radiatively driven entrainment in warming the boundary layer, relative to other processes, is found to substantially increase with CCN.

4. As a result, the turbulence shifts from surface driven to radiatively driven, with the surface energy budget being significantly reduced.

5. The strengthening of the thermal inversion plays a key role in the CCN impact on the entrainment efficiency.

The first two conclusions form a recommendation for future field campaigns in the Arctic that have a modeling spinoff in mind. Only when detailed measurements of in situ aerosol are taken throughout the boundary layer can numerical models satisfactorily reproduce the observed clouds in a realistic way. The absence of such aerosol datasets would introduce uncertainty in the simulations that makes it much harder to interpret their actual representativeness of observed conditions. The last three conclusions contribute to our understanding of air mass transformations, air–sea interactions, and feedbacks in Arctic climate. They imply that in all of these processes, the role of CCN can not be ignored.

The marine mixed-phase convection case as defined and examined in this study could serve as a prototype scenario for further studies, given the completeness of the RF20 dataset in covering thermodynamics, aerosol, and clouds. For example, microphysics schemes in LES models could be tested and critically assessed and compared to in situ cloud data. Other possible investigations could focus further on the lapse rate feedback. What would still be instructive is to investigate how often the stagnant convective conditions actually occur in the region. This could be achieved by establishing the frequency of occurrence of this regime in reanalysis or climate model data.

Appendix A: Composite forcings

Figure A1 shows the vertical structure of the time-constant forcings adopted for the LES experiments of ACLOUD RF20, as calculated from ECMWF analysis and short-range forecast data. Figure A1a and b show the prescribed tendencies of temperature ${\dot{t}}_{\mathrm{adv}}$ and humidity ${\dot{q}}_{\mathrm{adv}}$ due to large-scale advection. Figure A1c and d shows the zonal and meridional geostrophic wind speeds u_g and v_g, respectively. Figure A1e shows the pressure velocity Ω. The shadings indicates the 1st–99th and 25th–75th percentile spread among included ECMWF profiles, while the median is shown as a black dotted line.

Figure A1 Profiles of the prescribed large-scale forcings based on ECMWF data.

Appendix B: Microphysics scheme

Table B1 Overview of setting microphysical parameters for hydrometeors.

The overview covers the setting of microphysical parameters for the size and velocity of hydrometeors, as well as for particle mass distribution of hydrometeors under the assumption of a generalized gamma distribution.

B1 Nucleation of ice crystals

Ice nucleation in the bulk microphysics scheme of Seifert and Beheng (2006a) combines two previous approaches. Firstly, the number of activated ice nuclei is a function of supersaturation with respect to ice, as proposed by Meyers et al. (1992):

$$\begin{array}{}\text{(B1)}& N_{\mathrm{IN}}=N_{\mathrm{M}}\exp \left(a_{\mathrm{M}}+b_{\mathrm{M}}{S}_i\right)\mathrm{if}{S}_i>\mathrm{0}\mathrm{and}T<T_{\mathrm{M}},\end{array}$$

where T is the absolute temperature; S_i is the supersaturation with respect to ice surface; the values of parameters follow number parameter N_M=10³ m⁻³, the intercept coefficient $a_{\mathrm{M}}=-\mathrm{0.639}$, the linear coefficient b_M=12.96; and finally the threshold T_M=268.15 K limits below which temperatures does the ice nucleation occur. Secondly, in order to avoid very low number concentrations, the nucleation is limited to be within 1 order of magnitude from the modified Fletcher's formula (Fletcher, 1962; Reisner et al., 1998):

$$\begin{array}{}\text{(B2)}& N_{\mathrm{IN},\mathrm{F}}={\mathrm{10}}^{-\mathrm{2}}\exp \left(\mathrm{0.6}\left(T_{\mathrm{0}}-\max \left(T,T_{\min }\right)\right)\right)\mathrm{if}T<T_{\mathrm{M}},\end{array}$$

where T₀ is the freezing point of water and T_min=246 K limits production at extremely low temperatures.

B2 Freezing of hydrometeors

The timescale in the heterogeneous freezing of supercooled cloud droplets is dependent on their size and temperature, following Pruppacher and Klett (1997) and Khain et al. (2000). The timescale in the homogeneous freezing of cloud droplets is then given by Cotton and Field (2002). The heterogeneous freezing of raindrops is analogous to the heterogeneous freezing of cloud droplets. The only difference lies in classifying the resulting particle as graupel.

B3 Secondary ice production

The secondary ice production includes ice multiplication by the Hallett–Mossop process, occurring during the riming of ice hydrometeors in the temperature ranges between 265 and 270 K (Hallett and Mossop, 1974; Griggs and Choularton, 1986). The number of ice splinters released during the process is dependent on the temperature and the riming rate, following the parameterization of Beheng (1982).

B4 Nucleation of cloud droplets

The nucleation of cloud droplets again follows Seifert and Beheng (2006a) as closely as possible. Firstly, the nucleation rate is calculated explicitly as

$$\begin{array}{}\text{(B3)}& \begin{array}{rl}{\displaystyle \frac{\partial N_{\mathrm{c}}}{\partial t}}{|}_{\mathrm{nuc}}=& \stackrel{\mathrm{̃}}{C_{\mathrm{ccn}}}\mathit{\kappa }{S}_{\mathrm{l}}^{\mathit{\kappa }-\mathrm{1}}w{\displaystyle \frac{\partial S_{\mathrm{l}}}{\partial z}}\mathrm{if}{S}_{\mathrm{l}}\ge \mathrm{0},\mathrm{and}{S}_{\mathrm{l}}\le S_{\max }\\ & \mathrm{and}w{\displaystyle \frac{\partial S_{\mathrm{l}}}{\partial z}}>\mathrm{0},\end{array}\end{array}$$

where C_ccn and κ are CCN parameters, S_l is supersaturation with respect to liquid water surface (expressed in %), w is the vertical velocity of the air, and S_max is a threshold for saturation when all available CCN are activated. Secondly, the nucleation rate is limited so the number of droplets does not exceed N_ccn, the current CCN number concentration. While RF20 reflects maritime conditions, the values of the κ parameter and the saturation threshold are set to constant values of κ=0.462 and S_max=1.1 % (Seifert and Beheng, 2006a). Unlike in the original description, the parameter $\stackrel{\mathrm{̃}}{C_{\mathrm{ccn}}}$ is not the constant C_ccn, but it is instead dependent on the aforementioned variable N_ccn. The consistency with the power law relation for activation spectra is maintained by calculating this parameter as

$$\begin{array}{}\text{(B4)}& \stackrel{\mathrm{̃}}{C_{\mathrm{ccn}}}={\displaystyle \frac{\mathrm{1}}{{\left(S_{\max }\right)}^{\mathit{\kappa }}}}{N}_{\mathrm{ccn}}.\end{array}$$

The values of other important microphysical parameters are shown in Table B1.