Aerosol-cloud-turbulence interactions in well-coupled Arctic boundary layers over open water

Late springtime Arctic mixed-phase convective clouds over open water in the Fram Strait as observed during the recent ACLOUD field campaign are simulated at turbulence-resolving resolutions. The main research objective is to gain more insight into the coupling of these cloud layers to the surface, and into the role played by interactions between aerosol, hydrometeors and turbulence in this process. A composite case is constructed based on data collected by two research aircraft on 18 June 2017. The boundary conditions and large-scale forcings are based on weather model analyses, yielding a simulation 5 that freely equilibrates towards the observed thermodynamic state. The results are evaluated against a variety of independent aircraft measurements. The observed cloud macroand microphysical structure is well reproduced, consisting of a stratiform cloud layer in mixed-phase fed by surface-driven convective transport in predominantly liquid phase. Comparison to noseboom turbulence measurements suggests that the simulated cloud-surface coupling is realistic. A joint-pdf analysis of relevant state variables is conducted, suggesting that locations where the mixed-phase cloud layer is strongly coupled to the surface by 10 convective updrafts act as hot-spots for invigorated interactions between turbulence, clouds and aerosol. A mixing-line analysis reveals that the turbulent mixing is similar to warm convective cloud regimes, but is accompanied by hydrometeor transitions that are unique for mixed-phase cloud systems. Distinct fingerprints in the joint-pdf diagrams also explain i) the typical ringlike shape of ice mass in the outflow cloud deck, ii) its slightly elevated buoyancy, and iii) an associated local minimum in CCN. 15

westerly flow, featuring considerable variations in cloud cover and temperature in the mid-troposphere. Airborne observations were made during a wide range of cloud conditions, including both stably stratified and convective regimes. Both single-layer and multi-layer clouds structures were observed, but also clear-sky conditions.

RF20
This study exclusively focuses on Research Flight RF20 by the P5 and P6 aircraft on 18 June 2017. This flight took place in 95 the Fram Strait west of Svalbard (see Figure 1). As described by Knudsen et al. (2018), during this period the mid-troposphere experienced drying. High cirrus clouds were present during the previous day but disappeared on 18 June. The air mass in the Fram Strait was slow-moving on this day, being situated over relatively warm open water. Figure 1 illustrates that the cloud situation in the Fram Strait as encountered by the aircraft was relatively complex. However, a closer look suggests that a rough north-south regional division in cloud character can be made. In the north, over the sea ice 100 and its margin, the clouds were absent or very thin, allowing 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 anti-clockwise flight paths from their base in Longyearbyen (LYR) and visited these three 105 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 race-track", both aircraft flew back and forth between C5 and C6. While the P5 doubled back once and maintained altitude above the boundary layer inversion (at about 1.5 km height), the P6 doubled back twice, staying below inversion height. It maintained constant altitude for five brief segments to allow covariance calculations at multiple heights, as shown in Appendix A. In-situ in-cloud 110 measurements were also made by the P6 during this period (see Fig. A1). The enhanced and targeted sampling during the southern race-track 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.

Observational datasets
The observational data from the ACLOUD campaign used in this study are summarized in Table 1. Figure 2a shows detailed 115 measurements of the clouds in the target area as obtained with the Microwave Radar and radiometer for Arctic Clouds radar onboard 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 for 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.,120 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 the P5 three times, at around 16:20, 16:35 and   The eastern part of the target domain that is centered around dropsonde DS04 is selected as the area to be simulated (indicated 125 by the green box in Fig. 1 and the green line in Fig. 2). This choice is motivated by the following considerations. First, our main objective is to study interactions between convection, mixed-phase clouds and aerosols, which did occur in this area.
Second, the combination of warm water and cold air implies large surface latent and sensible heat fluxes, making the nearsurface 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 130 (LES). 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 hours (Fig. 1)  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.
Cloud measurements are of key importance for this study, given its focus on understanding interactions between micro-145 physics, turbulence and aerosol in mixed-phase convective clouds. In particular, the observational data on hydrometeor occurrence and mass as collected by the P5 and P6 are useful for evaluating LES experiments. While the MiRAC and AMALi data onboard P5 provide information about cloud top heights, vertical structure, and liquid water path, the Nevzorov probe onboard P6 provides in-situ samples of cloud liquid water content. The LaMP CIP probe provides information on ice water content using Brown and Francis (1995) mass diameter relationship on non-spherical particles only. Finally, the high-frequency (100 Hz) 150 turbulence measurements collected by the P6 noseboom (Hartmann et al., 2018) allow calculating (co)variances of temperature and vertical velocity in the boundary layer, even for relatively short flight segments.
During RF20 the P6 carried two Ultra-High Sensitivity Aerosol Spectrometers (UHSAS) to sample submicron-size particles (Zanatta and Herber, 2019). These instruments measured the number size distribution of particles with diameters between 60 nm and 1000 nm by detecting scattered laser light, using the method described in detail by Cai et al. (2008). Figure 4   155 shows that a relatively wide range of aerosol concentrations was encountered. The data suggests a marked difference in aerosol loading below and above the cloud layer, with the lower values found below the clouds. These measurements as documented  by Mertes et al. (2019) form the basis for the initial CCN profile as adopted in the simulations described in Section 3.4. The skill of the model in reproducing the observed vertical structure of aerosol will be assessed.

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) 165 and high latitudes (de Roode et al., 2019;Neggers et al., 2019;Egerer et al., 2020;Neggers, 2020a, b). The code of DALES is open source and maintained online at https://github.com/dalesteam/dales. 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, w}, total water specific humidity q t , liquid water potential temperature θ l , as well 170 as the number and mass of various hydrometeor species. The time-integration makes use of a 3rd-order Runge-Kutta approach. Scalar advection is represented using an upwind scheme, while the centered difference method is applied for the three velocity components. The subgrid-scale transport of heat, moisture and momentum is dependent on a prognostic TKE model. For 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.

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 (2006) (hereafter SB06) was recently implemented, adding a further three prognostic hydrometeors (cloud ice, snow and graupel). First DALES results with this scheme, including an evaluation against 180 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) are particularly relevant for this study.
A key difference with SB06 is the prognostic treatment of the number concentration of Cloud Condensation Nucleii (CCN).
This first applies to activation of CCN in saturated grid cells. The CCN concentration is conserved during nucleation of 185 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 a 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 190 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 (INP) and ignoring any removal (see Appendix C). 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 195 Hallet-Mossop process (Hallett and Mossop, 1974), occurring during the riming of ice hydrometeors in the temperature ranges between 265 K 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 200 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 205 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 C and Table C1.

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 UTC and the dropsonde launch (see Fig. 1). This reflects the 210 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 twelve-hour 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.

215
(2019), a combination of analyses (available every 12 hours) and short-range forecasts (available every 3-hours) is used, effectively yielding a four-dimensional dataset of the atmospheric state variables {θ l , q t , u, v} at 3-hourly temporal resolution and 0.1x0.1 degree 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 0.5 • × 0.5 • -wide column 220 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, timeaveraging is then applied over 12 hours to obtain the composite forcing dataset used to drive the LES. These profiles are 225 described in more detail in Appendix B), 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 initial state, amalgamated with the composite large-scale model data. Where available, the sonde data is averaged onto the vertical grid of the composite ECMWF forcings. At grid layers where no DS04 data is available, the composite model state itself is used instead. Figure 3 shows the initial profiles thus obtained, illustrating that 230 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, 235 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 short wave radiative flux at the lower boundary. The incoming short wave radiative flux at the model ceiling depends on i) seasonality and time of

Control experiment 245
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 hours, to provide the turbulent boundary layer with enough time to equilibrate and to cover three full 250 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.
A sponge layer is applied in the top 1 km of the computational domain to prevent reflection of waves off the rigid top boundary. In addition, continuous nudging towards the initial profile is applied above the boundary layer inversion to prevent 255 excessive model drift in this height range. A Newtonian relaxation term is included in the prognostic equations for {u, v, θ l , q t } to this purpose, adopting a timescale of 6 hours. A 200 m deep transition zone is included above the inversion across which the nudging intensity increases linearly with height. Note that 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   (Stevens et al., 2001;Stevens et al., 2020). It is also distinctly different from turbulent 280 mixed-phase clouds over homogeneous sea ice, which are often fully decoupled (e.g. Solomon et al., 2014). What sets these convective boundary layer cloud systems apart from their equivalents in warmer climes is the presence of cloud ice (not shown in this rendering). from warm marine stratocumulus, and is driven by daytime absorption of short wave radiation alternating with nighttime cloud 290 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.

Time evolution
Figure 6c-d 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. The evolution of the humidity flux shows that it takes about 12 hours for the surface-driven turbulence to properly spin up after 295 initialization. Once this has occurred the transport is continuous, with the local minimum in the θ v flux near liquid cloud base indicating the presence of a shallow transition layer (see Fig. 6c). 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.   Morrison et al., 2009). Some rain q r is present below liquid cloud base, but snow q s and graupel q g masses are much larger, peaking at a slightly lower height compared to the suspended hydrometeors. This feature, in combination with the almost discrete change at z ∼ 400 m from snow/graupel to 315 rain, suggests that significant precipitation forms near the inversion and then forms shafts extending to the surface. Precipitation melts when it crosses freezing level, losing mass on its way down due to sublimation and evaporation. Such shafts are actually visible in the radar measurements shown in Fig. 2. The results indicate that the maximum height at which hydrometeors are observed is well reproduced by the model. This is consistent with the good agreement on inversion height as detected by the dropsonde. The shape of the cloud top distribution is less well reproduced, although some liquid cloud tops are present in the model above z ≈ 500 m. If this is a model shortcoming, or rather results from the small observational sample size due to the sampling of only one convective event, remains a question that can not be conclusively answered at this point. For example, the convective structure that was observed might have been 330 rare, so that the observations are skewed towards the associated cloud properties. are avoided. Accordingly, the full size distribution of ice particles as detected by CIP is used. Note that the microphysics scheme does include secondary ice production, allowing the formation of large and heavy ice hydrometeors in the simulation. Figure 10a shows that in general the measured values of liquid hydrometeor masses are situated inside the pdf of LES values.

345
Also the simulated vertical structure seems realistic, with the largest values found near the inversion. Figure 10b indicates that the same applies to frozen hydrometeor mass, with the CIP measurements situated well within the simulated bandwidth, except for a single spike near the inversion. We conclude from these results that, apart from a few outlying samples, in general the simulated structure and amplitude of mixed-phase hydrometeor mass is realistic.

350
The measurements at 100 Hz of temperature and vertical velocity made by the sensors in the P6 noseboom during RF20  allow calculating variances of both variables, as well as the turbulent heat flux. Usually a dedicated flight strategy is adopted to this purposes, featuring segments at which the aircraft flies at a fixed altitude (Nicholls and Lemone, 1980). During RF20 the southern racetrack featured five segments at constant height within the lowest 2km, as shown in   Both the simulation and the noseboom observations suggest the existence of a textbook turbulent mixed-layer below cloud base, featuring a concave w 2 structure, a convex T 2 structure and a linear decrease in the associated heat flux towards slightly negative values just below cloud base (at about z = 500 m). The latter is consistent with the local minimum in the buoyancy flux seen in Fig. 6d. At the flight segment in the middle of the cloud layer w 2 and T w have a second maximum, which the LES also reproduces and reflects the impact on turbulence of latent heat release due to condensation. In general the domain 370 averaged LES profiles are situated within the observed 5-95 percentile range, which is encouraging. The measured ranges do We conclude from these results that the amplitude and vertical structure of both intensity and transport by convection and turbulence between surface and inversion is reasonably well captured by the LES.

Results II: Joint-pdf analyses
With a satisfactory agreement between model and measurements established, we now proceed and use the numerical results to investigate the interactions between aerosol, clouds and turbulence in more detail. To this purpose we will focus on single, 380 well-developed convective cells as appear in the last 24 hours of the simulation.

Spatial structure
Figure 12 digs deeper into the spatial structure of convection and clouds in the RF20 case. A convective cluster is sliced vertically, revealing moist plumes rising from close to the surface, condensing and subsequently reaching the stratiform mixedphase cloud layer just below the inversion. This is commensurate with the cloud macrophysical structure visible in Fig. 5. An interesting observation is that often no cloud ice is present in these core areas at all. In addition, the maxima in cloud 395 ice occur some distance away from the center of the cluster. As a result, the cloud ice has a distinct ring-like structure that surrounds the center of the cluster. This ring forms some times after cluster onset, and over time radiates out until it covers about half the domain. The ring is not present in liquid water, but does correlate with weak positive anomalies in q t and θ l . As shown in Fig. 12c), the ice ring also correlates with a slightly higher cloud top.

400
While only one convective event is discussed here, these clusters frequently occur and always behave similarly. These findings motivate interpreting such clusters and their outflow area as "convective hotspots", being places where the coupling to the surface by moist updrafts is strongest, where intense hydrometeor transitions take place, and where interactions with thermodynamic state and turbulence are strongest. How these interactions work is investigated next, using joint-pdf analyses based on scatterplots between various relevant variables. An example is shown in Figure 14, for the two-dimensional horizontal slice at 405 the 3d time point as shown in Fig. 13. Each dot represents the state of a single gridbox on the horizontal slice. The width of these probability density functions (pdfs) represents the variance at the level of interest. In addition, the shape and orientation of the joint-pdf directly reflect the underlying physical/dynamical mechanisms that create and maintain them. These diagrams have previously been used to investigate mixing between warm convective clouds and their environment (Paluch, 1979), but their use in investigating cold mixed-phase clouds is less established. Including hydrometeor properties in such analyses as an 410 extra (color) dimension can provide more insight into aerosol-cloud-turbulence interactions that are unique for mixed-phase clouds. Figure 14a shows a scatter plot of two thermodynamic state variables that are conserved for all phase changes of water, including the ice phase. These are the total specific humidity q t = q v +q l +q i and the liquid-ice water potential temperature θ li .
In contrast to its more well-known cousin θ l which only accounts for the latent heating of condensation L v , θ li also corrects  for the latent heating of fusion L i : where Π is the Exner function and c p is specific heat capacity of air at constant pressure. In the absence of other diabatic processes, the process of turbulent mixing directly shows up in these diagrams as a tight and elongated joint-pdf, situated between the two thermodynamic states between which the air is mixed. Figure 14a indicates that this feature, often referred to 420 as the "mixing line", is also present in this case. Comparison to the mean vertical profile (dashed black line) indicates that the mixing takes place between moist and cold air from close to the surface (black triangle) and a point on the mean profile (black plus) that is situated just above the local mean state at this height (black circle). The presence of a well-defined mixing line in this case demonstrates that the convective mixing process in mixed-phase cloud systems acts similar to warm cloud regimes (Neggers et al., 2002). A second cluster of datapoints can be distinguished, situated around the local mean state in a fairly 425 horizontal direction. The widths of this pdf in both directions reflects the turbulence in the stratiform cloud layer that is driven by cloud top cooling.
23 https://doi.org/10.5194/acp-2021-888 Preprint. Discussion started: 5 November 2021 c Author(s) 2021. CC BY 4.0 License. Figure 14b provides more insight into turbulence by considering vertical velocity w. One feature immediately stands out, which is the diagonally oriented tail towards large humidity and velocity values. These are the rising moist updrafts arriving at the inversion level, as also visible in Fig. 12. To trace these data points also in other panels such as Fig. 14a, a subset of 430 points is now defined using the double criterion w > 1 m s −1 and q t > 3.2 g kg −1 (chosen purely for visualization). These are indicated in orange in all panels. Figure 14a illustrates that all these updraft points sit in the top of the mixing line, and are associated with the strongest anomaly values in q t and θ li at this height that are still close to the subcloud layer values (black triangle). Their high correlation indicates strong upward fluxes of heat and moisture, representing turbulent transport that acts to maintain humidity below the inversion in the well-coupled boundary layer.

435
Calculating the intersection point of the mixing line with the mean profile yields the height of the air with which the convective updrafts are effectively mixing. This is achieved by least-squares fitting a line to the points with q t > 3.0, which reflects the lower end of the mixing-line cluster of points. The resulting fit is shown in Fig. 14a, yielding an intersection point at the height of z m = 1289 m (indicated by the black plus). This height is about 100 m above the level of diagnosis, which is consistent with earlier studies of this kind for warm convection (Blyth et al., 1988).

Hydrometeor transitions
The behavior of mixed-phase hydrometeors is investigated next, using the same technique. Figure 15a shows that a strong correlation exists between q t and q l , reflecting that points with highest q t also have largest q l . This is typically observed in moist updrafts that cool and condense while they are rising, as indicated by the orange points. The shape of the joint-PDF of q t and cloud ice q i (see Fig. 15b) is still highly correlated but is less tight, and also has a distinct chevron shape with a 445 tail extending towards high q t -low q i values. This tail is where the updrafts are located, as indicated by the orange points.
Apparently the highest cloud ice values are not found in the rising cores; put more strongly, the strongest updrafts are almost free of cloud ice. This result is consistent with the gaps in the ice cloud field visible in Fig. 13c. Accordingly, the chevron shape seems caused by the onset of glaciation in updrafts, with the highest ice masses in the tip on the right reflecting air that has recently experienced this transition in phase.

450
To gain more insight, the points that carry most ice mass are traced in these diagrams by applying the criterion q i > 0.08, a value inspired by this diagram. This subset of points is indicated in light blue in all panels. While the high q i points are still situated on the mixing line (see Fig. 14 a), they are no longer rising see Fig. 14 b) and have lower q l (see Fig. 15 c). This confirms the idea that the process of glaciation mainly takes place in air that was until recently part of an updraft, turning high q l values into ice. Geometrically, this creates a distinct ring of high q i anomalies in the radial outflow region around the updraft 455 core, as visible in Fig. 13b.

Dynamics and aerosol
Dynamics and aerosol characteristics are considered next. Figure 16a shows that the high q i points are relatively tightly clus-465 tered around a single point in (θ v , w) space, which is almost non-moving vertically but still being weakly buoyant at θ v = 0.3 K. Taking into account that these points represent ex-updraft air that has just glaciated, the weak positive buoyancy could have been caused by the latent heat of fusion. This shared experience could partially explain why the points with most ice mass all have a similar buoyancy and are tightly clustered in this diagram, because the latent heat released is directly proportional to the water mass that has changed phase.

470
Making an informed estimate can help in this respect. From Fig. 14a we estimate that, on average, about 0.4 g kg −1 of liquid water has been converted into ice, as measured by the distance on the horizontal axis between the midpoints of the orange and light-blue sets of data points. Multiplying this by the latent heat of fusion, L i = 3.3e 5 J kg −1 , and dividing by the specific heat capacity of air at constant pressure c p yields a temperature change of about 0.13 K. Comparing this to the θ v excess of about 0.3 K suggests that the latent heating can explain about 50 % of this temperature difference. Accordingly, we conclude 475 that other processes must also contribute to this buoyancy level. Snow and graupel formation removes glaciated water mass, a process that also increases gridbox buoyancy. Acting together, both the latent heating of fusion and the reduction in condensate loading due to frozen precipitation formation are thus responsible for the weakly positive buoyancy in the ice ring surrounding the cluster. This positive buoyancy, which is unique for mixed-phase convection, is probably also responsible for the slightly raised cloud top in the ice ring visible in Fig. 12c. Lc and L f stand for latent heat release due to condensation and fusion, respectively. B + indicates a state of heightened buoyancy, while N − s and N − g indicate loss of CCN due to snow and graupel formation, respectively. Existing stratiform cloud ice is indicated as "old ice", while newly glaciated cloud ice is labeled as "new ice". Mixing between updrafts and their environment is indicated in orange. The mesoscale circulation associated with the convective hotspot is shown as blue arrows. Black lines indicate either convection-related clouds (foreground) or existing stratiform ice clouds (background).
The final component of our joint-pdf analysis concerns the behavior of CCN in the cloud layer. Figure 16b shows the joint pdf between θ v and CCN. Note that in this experiment the CCN is prognostic, and can change by many processes. What immediately catches the eye is the highly correlated cluster of points in the low CCN-high buoyancy quadrant. Interestingly, all three subsets of glaciated hydrometeors sit in this tail, with the snow occupying the extremest regions. This behavior suggests that the formation of frozen precipitation is the main source of CCN depletion, with snow being the most efficient. We speculate 485 that the strong correlation with buoyancy is probably not causal, but reflects that CCN depletion and buoyancy boosting are both driven by the same process (glaciation and frozen precipitation formation), and happen at the same location (the ice ring). This is consistent with the local minimum in the CCN profile at this height, as shown in Fig. 17a.

Discussion
While the case configuration yields a simulation that equilibrates close to the boundary layer and cloud state as observed 490 during RF20, some aspects of the idealized experiment can have impact on the behavior of the simulated turbulence and clouds. These include the size of the turbulent domain, and the spatial and temporal discretization. The way the larger-scale forcing is represented can also affect the results. For example, one could adopt heterogeneous forcing in a nested setup, which allows representing advection of mesoscale features into the domain that are now ignored. We speculate that this can affect the circular structure of the convective clusters. Another important model component is the microphysics scheme, for which one 495 expect particular sensitivity. Closely related to microphysics is the treatment of CCN and INP. In a follow-up study the authors explore some of these potential sensitivities in the LES results for this ACLOUD case, focusing on microphysics-related aspects and interpreting these in the context of the ongoing warming of the Arctic climate.
What emerges from the joint-pdf analyses is a conceptual picture of how turbulence, microphysics and aerosol interact in well-coupled Arctic boundary layers over open water. This idea is schematically illustrated in Fig. 17b. Relatively long-lived 500 and strong convective clusters inject moist air into the mixed-phase stratiform cloud deck. When the associated liquid clouds stop rising, they gradually glaciate while fanning outwards, creating distinct "ice rings" in the stratiform cloud layer. The associated latent heating of fusion, in combination with the formation of frozen precipitation, then cause a weak buoyancy increase, which helps to further lift the outflow cloud. The precipitation that forms, consisting of both graupel and snow, is effective in locally removing aerosol. This probably causes the local minimum in the CCN profile immediately below the 505 thermal inversion (see Fig. 17a), which is then mixed across the whole convective boundary layer. This vertical structure is in agreement with the ACLOUD observations (see Fig. 4), as well as other observational studies (Fitch and Garrett, 2020). The intensity of these processes is largest inside these mixed-phase convective clusters and their outflow region, which motivates interpreting them as "hotspots" or engines for enhanced aerosol-cloud-turbulence interactions in this climate regime.

510
In this study an LES experiment is performed for conditions as observed during Research Flight 20 of the ACLOUD campaign in June 2017 over open water west of Svalbard. The simulation features a mixed-phase cloud layer that is well-coupled to the surface by relatively long-lived convective clusters in liquid phase. Evaluation against in-situ and remote-sensing measurements collected by the P5 and P6 aircraft indicates that the thermodynamic, kinematic, turbulent and cloud structure are well reproduced. Simulation output is then subjected to a joint-pdf analysis, giving insight into aerosol-cloud-turbulence interactions 515 that take place inside the mixed-phase convective clouds. The analysis confirms that the mixing-line paradigm as well-known from warm convection also holds for mixed-phase clouds, suggesting that the mixing source of rising convective updrafts is about 100 m above them. We find that the glaciation of ex-updraft air creates distinct ring-like patterns in the stratiform ice cloud layer. The associated heat of fusion drives local buoyancy, which acts to further lift cloud top in these areas. This process is unique for mixed-phase cloud layers fed by moist convection. Snow forms where cloud ice mass is large, while graupel 520 predominantly forms in updrafts. Both forms of frozen precipitation further contribute to positive buoyancy in the stratiform layer. Precipitation is also efficient in locally removing CCN, creating a distinct minimum in its profile.
An important conclusion from this study is that the joint-pdf analysis in conserved variable space, best known from its previous application in studies of entrainment in warm convective clouds, can be similarly effective when applied to cold or mixed-phase clouds. It seems particularly suited to investigate interactions between aerosol, hydrometeors and turbulent 525 dynamics. This motivates the further testing of this method for other cold cloud regimes, such as clouds at high latitudes and/or the tenuous cloud regime. Another follow-up activity is to find out if joint-pdf analyses of purely observational datasets, for example high-frequency in-cloud measurements of multiple variables by aircraft, exhibits the same fingerprints. This is for now considered a future research topic. We hope the results obtained in this modeling study can provide a useful context for this effort.

530
The obtained results and insights also have a bearing on Arctic Amplification. One wonders how often the stagnant wind conditions as encountered during ACLOUD RF20 actually occur in the region, or more globally speaking, in areas of open water off the sea-ice edge. Such information is needed to fully understand the role in Arctic climate played by the "convective hotspots" for aerosol-cloud-turbulence interaction as identified in this study. The effectiveness of meridional transport of aerosol and pollution into the high Arctic is affected by the occurrence of such convective events on the way. The associated overturning 535 plays a role in the transformation of Arctic air masses, which have been identified as a key process in Arctic climate change (Pithan et al., 2018). Exploring these topics, potentially using dedicated LES experiments based on field campaign data, is for now considered a future research activity that could be inspired by the outcome of this study.
Appendix A: P6 flight details Figure A1 shows various properties recorded by the P6 aircraft within the southern race-track section of RF20, situated inside 540 the target area of this study (as indicated by the orange box in Fig. 1). The section included five flight legs at constant height, used for calculating the observed covariances as shown in Fig. 11. The third panel shows the cloud liquid water content measured by the Nevzorov probe as used to evaluate the LES in Fig. 10. Figure B1 shows the vertical structure of the time-constant forcings adopted for the LES experiments of ACLOUD RF20, as 545 calculated from ECMWF analysis and short-range forecast data. Panels a) and b) show the prescribed tendencies of temperaturė t adv and humidityq adv due to large-scale advection. Panels c) and d) shown the zonal and meridional geostrophic wind speeds u g and v g , respectively. Panel e) shows the pressure velocity Ω. The shadings indicates the 1-99 and 25-75 percentile spread among included ECMWF profiles, while the median is shown as a black dotted line.

C1 Nucleation of Ice Crystals
Ice nucleation in the bulk microphysics scheme of Seifert and Beheng (2006) 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): where T is the absolute temperature, S i is the supersaturation with respect to ice surface, and the values of parameters follow: 555 number parameter N M = 10 3 m −3 , the intercept coefficient a M = − 0.639, and 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 one order of magnitude from the modified Fletcher's formula (Fletcher, 1962;Reisner et al., 1998):  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 generalised gamma distribution.
where T 0 is the freezing point of water and T min = 246 K limits production at extremely low temperatures.

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

C3 Secondary Ice production
The secondary ice production includes ice multiplication by Hallet-Mossop process, occurring during the riming of ice hydrometeors in the temperature ranges between 265 K and 270 K (Hallett and Mossop, 1974;Griggs and Choularton, 1986).

570
The number of ice splinters released during the process is dependent on the temperature and the riming rate, following the parameterization of Beheng (1982).

C4 Nucleation of Cloud Droplets
The nucleation of cloud droplets again follows Seifert and Beheng (2006) as closely as possible. Firstly, the nucleation rate is calculated explicitly as 575 ∂N c ∂t nuc = C ccn κ S κ−1 l w ∂S l ∂z if S l ≥ 0, and S l ≤ S max , and w ∂S l ∂z > 0 (C3) 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, 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 κ parameter and the saturation threshold are set to constant values κ = 0.462 and S max = 580 1.1 % (Seifert and Beheng, 2006). Unlike in the original description, the parameter C 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 The values of other important microphysical parameters are shown in the Table C1.