Model emulation to understand the joint effects of ice-nucleating particles and secondary ice production on deep convective anvil cirrus

Ice crystal formation in the mixed-phase region of deep convective clouds can affect the properties of climatically important convectively generated anvil clouds. Small ice crystals in the mixed-phase cloud region can be formed by heterogeneous ice nucleation by ice-nucleating particles (INPs) and secondary ice production (SIP) by, for example, the Hallett–Mossop process. We quantify the effects of INP number concentration, the temperature dependence of the INP number concentration at mixed-phase temperatures, and the Hallett–Mossop splinter production efficiency on the anvil of an idealised deep convective cloud using a Latin hypercube sampling method, which allows optimal coverage of a multidimensional parameter space, and statistical emulation, which allows us to identify interdependencies between the three uncertain inputs. Our results show that anvil ice crystal number concentration (ICNC) is determined predominately by INP number concentration, with the temperature dependence of icenucleating aerosol activity having a secondary role. Conversely, anvil ice crystal size is determined predominately by the temperature dependence of ice-nucleating aerosol activity, with INP number concentration having a secondary role. This is because in our simulations ICNC is predominately controlled by the number concentration of cloud droplets reaching the homogeneous freezing level which is in turn determined by INP number concentrations at low temperatures. Ice crystal size, however, is more strongly affected by the amount of liquid available for riming and the time available for deposition growth which is determined by INP number concentrations at higher temperatures. This work indicates that the amount of ice particle production by the Hallett–Mossop process is determined jointly by the prescribed Hallett–Mossop splinter production efficiency and the temperature dependence of ice-nucleating aerosol activity. In particular, our sampling of the joint parameter space shows that high rates of SIP do not occur unless the INP parameterisation slope (the temperature dependence of the number concentration of particles which nucleate ice) is shallow, regardless of the prescribed Hallett–Mossop splinter production efficiency. A shallow INP parameterisation slope and consequently high ice particle production by the Hallett– Mossop process in our simulations leads to a sharp transition to a cloud with extensive glaciation at warm temperatures, higher cloud updraughts, enhanced vertical mass flux, and condensate divergence at the outflow level, all of which leads to a larger convectively generated anvil comprised of larger ice crystals. This work highlights the importance of quantifying the full spectrum of INP number concentrations across all mixed-phase altitudes and the ways in which INP and SIP interact to control anvil properties. Published by Copernicus Publications on behalf of the European Geosciences Union. 17316 R. E. Hawker et al.: Effects of mixed-phase ice production on convective anvils


Introduction
Deep convective clouds are an important component of the global hydrological cycle and radiative budget (e.g. Lohmann et al., 2016;Massie et al., 2002). The anvil cirrus cloud they produce can persist in the atmosphere for several hours to a few days and therefore impact outgoing radiation long after the deep convection has decayed (Luo and Rossow, 2004). However, accurately representing the spatial and temporal complexity of large convective systems and therefore convectively generated cirrus presents extensive challenges for atmospheric modelling (Prein et al., 2015).
Deep convective cloud systems extend vertically from the boundary layer to the tropopause and can have a horizontal radius of over 1000 km. They are dynamic and powerful systems with updraught speeds of up to 50 m s −1 (Frank, 1977;Musil et al., 1986;Xu et al., 2001). In addition, a multitude of different thermodynamic and microphysical conditions can exist within the same system. There is also a scarcity of measurements of these climatically important clouds, particularly profile measurements within the convective core (Fan et al., 2016), and thus a scarcity of data with which to validate representations of deep convective clouds in models. The myriad of competing microphysical processes operating within deep convective clouds, along with the difficulty in validating model simulations against observations, cause the simulation of deep convective clouds to be subject to a large number of parametric and structural uncertainties (Johnson et al., 2015;. In particular, mixed-phase microphysics presents a challenge for cloud modelling because it is critical for deep convective cloud properties and poorly understood (Prein et al., 2015).
One of the largest uncertainties in quantifying aerosolcloud interactions and the resultant climate impacts is the amount of, and balance between, liquid and ice in mixedphase clouds. In particular, the representation of microphysical processes affecting cloud phase in tropical convection contributes substantial uncertainty to the simulated climate response to global warming in climate models (Medeiros et al., 2008;Stevens and Bony, 2013). The representation of the amount of ice within deep convective clouds is also important for the representation of the amount and intensity of precipitation, the prediction of which is one of the most socially and economically important roles of numerical weather forecasting (Arakawa, 2004;Prein et al., 2015).
Within the mixed-phase region of deep convective clouds, i.e. the region between 0 and ∼ −38 • C where both liquid and ice can coexist, we hypothesise that three factors controlling ice production strongly influence the partitioning of condensate into cloud liquid and ice, which are listed as follows.
i. The limiting number concentration of ice-nucleating particles (INPs: aerosol particles with the ability to initiate the freezing of cloud droplets at mixed-phase temperatures) at the top of the mixed-phase cloud regime (in this work, this is the INP number concentration at −38 • C). In an aerosol population made up entirely of dust particles, all with ice-nucleating potential, the limiting number concentration of INP would equate to the number concentration of dust.
ii. The temperature dependence of the INP number concentration at mixed-phase temperatures (the rate of increase in ambient INP number concentrations as temperature decreases from ∼ 0 to ∼ −38 • C). This determines the concentration of INP at lower mixed-phase altitudes and therefore the altitude of liquid depletion due to heterogeneous freezing in the lower and middle mixed-phase cloud levels (e.g. Hawker et al., 2021;Takeishi and Storelvmo, 2018).
iii. The efficiency of ice production by secondary ice production (SIP) mechanisms, whereby small ice particles are produced from existing frozen hydrometeors , such as ice crystals frozen heterogeneously, or larger snow and graupel particles.
The limiting number concentration of INPs in the atmosphere is extremely variable and depends on several interacting factors. For example, Saharan dust is an efficient INP at temperatures below −15 • C and the largest component by mass of the global aerosol budget (Tang et al., 2016;Textor et al., 2006). The export of this atmospherically important INP, across the Atlantic Ocean, varies hugely depending on factors such as season (Ridley et al., 2012); desert soil moisture (Laurent et al., 2008); local wind speed (Grini et al., 2005;Laurent et al., 2008); and the occurrence and intensity of convection, wet removal, and dry deposition (Bou Karam et al., 2014;Marsham et al., 2011;Provod et al., 2016), both in source (Heinold et al., 2013) and transport regions (Sauter et al., 2019;Twohy and Twohy, 2015). As a result of variations in dust emission and transport, summertime INP number concentrations in the Saharan Air Layer can vary by up to 4 orders of magnitude at −33 • C (Boose et al., 2016). Variations in INP number concentrations can impact cloud properties and cloud radiative forcing (e.g. Shi and Liu, 2019;Solomon et al., 2018). However, the reported effect of changes to INP number concentrations on cloud properties can be non-linear, counterintuitive, or conflicting depending on the environmental conditions, magnitude of the tested perturbation, or study methodology (Deng et al., 2018;Fan et al., 2010a, b;Gibbons et al., 2018;Hawker et al., 2021;van den Heever et al., 2006;Phillips et al., 2005Phillips et al., , 2007. The temperature dependence of INP number concentration depends, amongst other factors, on the aerosol type providing INP in a given scenario. A large number of aerosol types have the ability to act as INPs, including mineral dust (Atkinson et al., 2013;Niemand et al., 2012;Price et al., 2018;Welti et al., 2018), organic material in sea spray (Mc-Cluskey et al., 2018;Wilson et al., 2015), bacteria (Šantl-Temkiv et al., 2015), and pollen (Diehl et al., 2002). INPs comprised of marine organics emitted with sea spray tend to have a shallower temperature dependence than INP comprised of mineral dusts. This means marine organic INPs tend to have a higher ice-nucleating ability than mineral dust INPs at warm temperatures but lower ice-nucleating ability at colder temperatures (Atkinson et al., 2013;DeMott et al., 2016;McCluskey et al., 2018;Niemand et al., 2012;Vergara-Temprado et al., 2017;Wilson et al., 2015). In numerical weather and climate models, the temperature dependence of INP number concentration can be described by the slope of the INP parameterisation (i.e. d(log 10 N INP )/dT , as described in Hawker et al., 2021). The INP parameterisation slope depends on aerosol type (DeMott et al., 2010;Harrison et al., 2016Harrison et al., , 2019 and any ageing the aerosol has been subjected to (Boose et al., 2016;Brooks et al., 2014) as well as particle properties yet to be fully understood, such as surface morphology . The INP slope of any one aerosol population (composed of different INP types) is extremely uncertain and difficult to accurately predict without specific measurements. Variation in ice nucleation active site densities (n s ) even of materials of similar mineralogy can span several orders of magnitude at any one temperature (Atkinson et al., 2013;Harrison et al., 2016Harrison et al., , 2019. The factors governing active site location and densities are not fully understood but are theorised to be related to features on an INP surface such as surface pits , hydrophilic sites (Freedman, 2015), or lattice mismatches (Kulkarni et al., 2015). Variations in the temperature dependence of INP number concentration can affect the cloud development and the altitude at which cloud glaciation occurs, as was noted by Takeishi and Storelvmo (2018). This difference in glaciation altitude has been shown to cause differences in hail amount, intensity, and size (Liu et al., 2018); anvil ice crystal number concentration (N ice ) (Takeishi and Storelvmo, 2018); and radiative forcing (Hawker et al., 2021) of convective clouds.
Observational campaigns have long documented the existence of ice crystals at concentrations vastly exceeding the concentration of INP in clouds with relatively warm cloud top temperatures, indicating the presence of SIP mechanisms (e.g. Crawford et al., 2012;Field et al., 2017;Huang et al., 2017;Lasher-Trapp et al., 2016). SIP can occur via processes such as rime splintering (i.e. the Hallett-Mossop process), droplet shattering, collision fragmentation, and sublimation fragmentation Korolev and Leisner, 2020). The most well-studied SIP mechanism is the Hallett-Mossop process by which small ice splinters are produced during the riming of liquid drops onto existing frozen hydrometeors (Crawford et al., 2012;Field et al., 2017;Hallett and Mossop, 1974;Phillips et al., 2007). However, even the Hallett-Mossop process is relatively poorly defined and its importance disputed. A recent laboratory study failed to observe rime splintering in conditions designed to stimulate the Hallett-Mossop process (Emersic and Connolly, 2017), and some recent literature suggests that previous observations of ice crystal number concentration attributed to the Hallett-Mossop process may have been indicative of other secondary ice formation mechanisms . Nevertheless, as it is the only SIP mechanism that is currently represented in most numerical weather prediction (NWP) models, we focus on the uncertainty associated with the Hallett-Mossop process in this study.
In addition to the individual uncertainties in INP number concentration, INP temperature dependence, and SIP rates, these three factors can also interact causing non-linear or counterintuitive changes in cloud properties, further motivating the exploration of their combined effects here. For example, intermediate INP number concentrations and intermediate Hallett-Mossop ice production rates have been found to produce higher cloud ice crystal number concentrations (ICNC) than high INP number concentrations or high SIP rates alone due to non-linear interactions between the two freezing mechanisms whereby very high heterogeneous freezing rates affect the availability and efficiency of secondary ice production, and vice versa (Crawford et al., 2012;Sullivan et al., 2017).
We investigate the individual and interacting effects of the INP number concentration, the temperature dependence of INP number concentration across the full spectrum of mixedphase temperatures, and the Hallett-Mossop ice production efficiency on the micro-and macro-physical properties of an idealised deep convective cloud by conducting a large ensemble of idealised simulations using Latin hypercube sampling to select the input parameter combinations and using statistical emulation where appropriate to analyse the ensemble output. We quantify the importance of the three uncertain input parameters and their interactions with one another for the anvil properties of the simulated deep convective cloud. Our methodology proves to be a powerful tool for analysing and understanding the behaviour of complex systems (Johnson et al., 2015;Lee et al., 2011;Marshall et al., 2019; because Latin hypercube sampling enables dense sampling over a defined parameter uncertainty space, allowing more extensive coverage of the defined parameter space than a traditional one-at-a-time test, and statistical emulation allows the production of detailed response surfaces of system behaviours across the three dimensions. This paper is structured as follows: Sect. 2 describes the idealised cloud model and the simulation set-up, as well as the methods used in our analysis. In Sect. 3, we examine the role of the uncertain input parameters in determining the ice crystal number concentration, ice crystal size, and the cloud fraction of the simulated deep convective anvil cirrus. In Sect. 4, we detail the limitations of our study. Section 5 summarises the main findings and implications of this study.

Model set-up and simulation design
This work utilises the Met Office NERC Cloud Model (MONC), which is a large Eddy simulation (LES) model with interactive cloud microphysics and radiation. While, the underpinning science of MONC is based on the Met Office Large Eddy Model (LEM) (Gray et al., 2001), MONC is a complete redesign of the LEM, which incorporates a pluggable component architecture to improve usability and is designed to be highly scalable (Brown et al., 2015(Brown et al., , 2018. Here, MONC is coupled to the Met Office Cloud AeroSol Interacting Microphysics (CASIM) module, which is a multimoment bulk scheme. MONC-CASIM has been used to investigate aerosol-cloud interactions in nocturnal fog (Poku et al., 2019) and low-level clouds during the West African monsoon season . CASIM has also been used with the Met Office Unified Model in regional simulations of coastal mixed-phase convective clouds (Miltenberger et al., 2018a, b), south-east Pacific stratocumulus clouds (Grosvenor et al., 2017), Southern Ocean supercooled shallow cumulus , mid-latitude cyclones (McCoy et al., 2018), cloud condensation nuclei (CCN)-limited Arctic clouds (Stevens et al., 2018), and tropical convective clouds (Hawker et al., 2021).
The simulations presented use a grid box spacing of 250 m (500 × 500 grid boxes) and 138 vertical levels. The model diagnostics are output every 5 min, and the time step is flexible to maintain model stability with a maximum value of 2 s and a minimum value of 0.01 s. MONC-CASIM is configured to be a two-moment scheme in this work. The number and mass concentrations for cloud droplets, rain droplets, ice crystals (or cloud ice), graupel, and snow are prognostic variables. The prognostic aerosol variables utilised in this work are the soluble accumulation-mode aerosol mass and number concentrations and the coarse mode dust mass and number concentrations. The aerosol can be advected around, but in the simulations presented here we choose to switch off scavenging processes, and the aerosol is therefore not incorporated into the cloud droplets when activated. The model boundary conditions are cyclical, and as such scavenging the aerosol would result in a rapid removal of all aerosol from the simulation.
The CASIM model configuration is very similar to that of Hawker et al. (2021). Cloud droplet nucleation is parameterised according to Abdul-Razzak and Ghan (2000). The soluble accumulation mode aerosol is used for cloud droplet activation, and a simplistic CCN activation parameterisation is included for the insoluble aerosol mode that assumes a 5 % soluble fraction on dust. Condensation is represented using saturation adjustment, meaning that where water saturation is exceeded at the end of a time step, the specific humidity is adjusted to be the equilibrium saturation over water and the grid box temperature and liquid mass is adjusted accord-ingly. If only frozen hydrometeors are present in a grid box, saturation is treated explicitly. Collision-coalescence, riming of ice crystals producing graupel, and aggregation of ice crystals producing snow are represented. Rain drop freezing is described using the parameterisation of Bigg (1953). Deposition onto ice and snow is treated explicitly allowing ice particles to grow in ice-supersaturated conditions including in the presence of liquid. Heterogeneous freezing (via immersion freezing) is active between −38 and −3 • C in the MONC-CASIM model used in this work and is described in more detailed in Sect. 2.2.1 and 2.2.2. The INP parameterisations resulting from the perturbations described in Sect. 2.2.1 and 2.2.2 inspect the conditions (temperature, cloud droplet number, ICNC) and aerosol concentrations within a grid box and use that information to predict an ice production rate via heterogeneous freezing. The supercooled droplets are depleted by the freezing parameterisation, but scavenging of INPs is not represented. As stated above, inclusion of scavenging was not possible as due to the cyclical boundary conditions of the simulation, scavenging processes would result in the rapid depletion of aerosol from the domain. The consideration of the number of ice crystals already present as well as the number of INP available when calculating the rate of heterogeneous freezing in a grid box acts as a control on the number of heterogeneous ice crystals forming in the absence of scavenging. Homogeneous freezing of cloud droplets is parameterised according to Jeffery and Austin (1997).
Radiative processes are represented by the Suite of Community RAdiative Transfer codes based on Edwards and Slingo (SOCRATES) (Edwards and Slingo, 1996;Manners et al., 2017), which in this study considers all five hydrometeor types for the calculation of cloud radiative properties. Changes in size and number of cloud droplets are considered. The cloud droplet single scattering properties are calculated from the cloud droplet mass and effective radius in each grid box using the equations detailed in Edwards and Slingo (1996). A fixed effective radius of 30 µm for ice crystals is used in the radiation calculations. For the other hydrometeor types (snow, graupel, rain), SOCRATES considers changes in mass but does not explicitly consider changes in number concentration or size (though changes in number and size will affect mass concentrations which are considered). As the use of SOCRATES in MONC including the radiative effects of ice hydrometeors has not been extensively tested, we do not examine the radiative diagnostic outputs or present them in this paper. Large-scale wind shear is prescribed and constant throughout the simulation. The Coriolis force is inactive in the simulations, and large-scale subsidence is determined by the local column theta.
We simulate a single deep convective cloud using the MONC-CASIM model. The cloud formation is initiated using a single warm bubble with a radius of 20 km, a height of 500 m, and a temperature perturbation of 1.5 • C. The model was initiated using mean profiles (wind velocity and direc- Figure 1. Initial conditions. The potential temperature and specific humidity (a) and wind speed and direction (b) profiles used to initiate the model. The profiles shown were extracted from a Met Office Unified Model simulation of a large deep convective cloud field in the maritime tropical Atlantic (described in Hawker et al., 2021). The profiles were averaged over out-of-cloud areas between 12:00 and 18:00 UTC. tion, potential temperature, specific humidity, and soluble accumulation mode aerosol number and mass concentration) extracted from a Met Office Unified Model simulation of a deep convective cloud field sampled during the "Ice in Clouds Experiment -Dust" flight campaign on the 21 August 2015 (between 12:00 and 15:00 UTC) (Hawker et al., 2021). Details of this simulation including comparisons to observations are available in Hawker et al. (2021). The environmental conditions used to initiate the model are shown in Fig. 1.
The simulation produces a large convective cloud with an extensive anvil (Fig. 2). Figure 3 shows that the cloud evolution for all simulations is similar with a large increase in surface precipitation (Fig. 3a) from 60 to up to 90 min and a decline that begins between 70 and 90 min. Similarly, the maximum cloud top height for most simulations peaks at around 120 min and afterwards declines slightly, indicating a reduction in convective strength (Fig. 3b). It is important for statistical emulation (Sect. 2.4), where one value for each cloud response is extracted from the model, that the clouds in each simulation undergo similar life cycles. We can see from Fig. 3 that this is the case for the simulated deep convective cloud.
When extracting the diagnostic variables and single values to be used for analysis, results from 60 to 180 min into the simulation are used to represent the convective cloud state. Maximum updraught speeds in the convective cloud period range from 30 to 50 m s −1 . Sixty minutes is approximately the time when the cloud top height first reaches the altitudes where freezing can occur (∼ 4 km, Fig. 2b) and therefore where the perturbations to the chosen uncertain input parameters (Sect. 2.2) are expected to start causing divergence between simulations. For focussing on the anvil stage of cloud development, we use the results from between 150 and 240 min in the simulation. There is no change in the model parameters or the forcing between the convective and anvil states, but rather the distinction between the life cycle stages is determined from the cloud evolution. The end of the convective cloud stage is determined by the time when the convective plume has largely decayed (by ∼ 180 min), and the beginning of the anvil cloud stage is determined by the time when a substantial anvil has formed (by ∼ 150 min) (Figs. 2 and 3). Table 1 lists the target output response variables that are investigated and the time period from which they are extracted. Unless a specific altitude is stated, or shown in a figure, the cloud properties shown for hydrometeor number concentrations, ice particle production rates, and cloud condensate values herein and listed in Table 1 refer to the mean integrated column value. For example, anvil ICNC is the mean value of the integrated number of ice crystals in Table 1. Target output variables. List of target output variables discussed in this study and the criteria used to extract their values from the simulation output.

Output variables Criteria
Anvil cloud stage Anvil ice crystal number concentration (ICNC) and size (diagnosed using effective radius defined as the ratio of the third to the second moment of the size distribution).

Convective cloud stage
Ice particle production rates, accretion rates, hydrometeor water paths, updraught speed Cloud condensate mixing ratio > 1 × 10 −6 kg kg −1 (i.e. in-cloud), time period in simulation: 60-180 min. all model columns of the anvil cloud (isolated using the criteria listed in Table 1).

Input parameters and their uncertainty ranges
In this work, we investigate the effect of variations in limiting INP number concentration, INP parameterisation slope, and the efficiency of ice particle production by the Hallett-Mossop process. For the purposes of this study, the magnitudes of these three factors are varied using the following uncertain input parameters.
-The limiting INP number concentration, termed N INP −38 herein, is the total number of aerosol particles capable of nucleating ice at the very top of the heterogeneous freezing regime (i.e. at a temperature of −38 • C).

The value of N INP
−38 is reported for the peak number concentration of the INP layer which is assumed to be transported to all cloud levels due to the strength of the applied warm bubble and resultant updraught.

INP parameterisation slope (λ INP )
λ INP is defined as the change in the log 10 of the INP number concentration per degree Celsius change in temperature as defined in Eq. (1): where N INP is the INP number concentration in m −3 and T is the temperature in degrees Celsius. Equation (1) is applied at temperatures between −38 and −3 • C. The range of λ INP values are calculated by varying the exponent (P , units of m −2 • C −1 ) in Eq.
(2) below, which determines the number of active sites per unit area of an aerosol population at temperature T , from −1.3 and −0.1. For this study, we define the number of active sites, n s , as where i is the intercept of the natural logarithm of n s (in active sites m −2 ) at 0 • C and T is the ambient tempera- 3), which is slightly steeper than that of the Atkinson et al. (2013) parameterisation based on K-feldspar. The maximum (shallowest) value of λ INP is −0.0434 • C −1 (P = −0.1), which is slightly shallower than that of the Meyers et al. (1992) parameterisation. The minimum (steepest) and maximum (shallowest) slopes simulated in this work are shown in Fig. 4a for a dust number concentration of 1 cm −3 with a mean radius of 1 µm. Using n s (T ), we can calculate the INP number concentration in m −3 at temperature T (N INP T ) as follows: where S is the surface area of the available INP in m 2 m −3 and N INP −38 is the limiting INP number concentration in m −3 as defined in Sect. 2.2.1. The minimum function shown in Eq. (3) is not used in this paper at temperatures other than −38 • C due to the parameterisation change described below.
In addition to varying the exponent, the original Niemand et al. (2012) parameterisation is altered in this work to allow the intercept at 0 • C to be flexible to ensure that the INP number concentration declines constantly from 0 to −38 • C.  was necessary to satisfy the assumptions of statistical emulation (Sect. 2.4) and to allow us to determine whether it is the limiting INP number concentration (e.g. total dust number concentration where dust is the only ice-nucleating material present in an aerosol population) or INP efficiency (e.g. whether the aerosol population is made up of marine organics or dust particles) that controls the properties of a deep convectively generated anvil cloud.
From Eqs. (2) and (3), we can see that Setting T to −38 • C, the intercept (i in Eq. 2) of the INP parameterisation can be calculated as follows:

The Hallett-Mossop process ice production efficiency (HM-eff)
The HM-eff in the model is varied from 1 to 1000 splinters produced per milligram of rimed liquid. The default efficiency of splinter production from the Hallett-Mossop process in MONC-CASIM is 350 mg −1 . This value is the best estimate of ice production based on a number of laboratory studies and has been used in previous modelling studies (Connolly et al., 2006;Hallett and Mossop, 1974;Mossop, 1985). However, other rates have been reported. An upper limit of 1000 mg −1 aligns with previous modelling studies where the efficiency of ice production by the Hallett-Mossop process was varied (Connolly et al., 2006). This upper limit also allows us to account somewhat for the possibility that the Hallett-Mossop process operating in real clouds is stronger than that observed in laboratory studies Takahashi et al., 1995).

Selection of uncertain input parameter combinations
MONC was run with combinations of values of N INP −38 , λ INP , and HM-eff from within the ranges shown in Table 2. Combinations of parameter values were defined using a maximin Latin hypercube design algorithm. Latin hypercube sampling is based on the Latin square and ensures optimum space filling (Johnson et al., 2015;Lee et al., 2011;Mckay et al., 2000). The maximin algorithm maximises the minimum distance between points in the cube (Lee et al., 2011). The use of a maximin Latin hypercube sampling design means that the parameter combinations cover the threedimensional parameter space in an optimum way. We can therefore evaluate the full effects of the parameters (individual and interacting) using traditional analysis on just the simulation data themselves, as well as employing statistical emulation (described in Sect. 2.4) to produce response surfaces of cloud properties.
The applied parameter values are shown in Fig  to compensate for this, and the parameter combinations of the additional simulations were selected by augmenting points into the largest gaps in the realistic section of the original Latin hypercube design.

Statistical emulation of the model output
Statistical emulation is a "process by which the computer model is replaced by a statistical surrogate model that can be run more efficiently" (Lee et al., 2011). This approach has previously been used to look at deep convective cloud microphysical properties in a 3D model (Johnson et al., 2015), hail formation , nocturnal stratocumulus (Glassmeier et al., 2019), and aerosol forcing from volcanic eruptions (Marshall et al., 2019). In this study, as well as using traditional methods of analysis, we explore the usefulness of statistical emulation as a tool to understand the interacting effects of mixed-phase ice production mechanisms.
Statistical emulation (as well as the applied Latin hypercube sampling methodology) has advantages over traditional one-at-a-time tests (where one variable is varied at predictable values from a control or base case while all other variables are held constant). Firstly, it allows the exploration of the effects of simultaneously perturbing multiple uncertain input parameters on output variables of interest across the entirety of reasonable parameter space for a much reduced number of complex simulations. Secondly, dense sampling via statistical emulation enables techniques such as variancebased sensitivity analysis to be applied, through which we can identify the input parameters that are contributing the most uncertainty to important output responses. This subsequently allows for the direction of resources towards quantifying and accurately representing those key parameters that contribute large amounts of uncertainty to output variables of interest.
We use a Gaussian process as the basis for the emulator (Johnson et al., 2015;Lee et al., 2011;Marshall et al., 2019). A more detailed description of the process used to construct emulators can be found in Johnson et al. (2015) and Lee et al. (2011). Separate Gaussian process emulators are built from the output of the training simulations (i.e. the emulation runs shown in Fig. 5) for each of the target output variables listed in Table 1 using the statistical software R (R Core Team, 2017) and the DiceKriging package (Roustant et al., 2012). These emulators are 3D maps of how the values of the target output variables change in the simulation output depending on the three uncertain input parameters (Sect. 2.2); i.e. they are surrogate statistical representations of the MONC-CASIM model. The emulators assume a linear mean function including all uncertain inputs and a Matérn covariance structure. The Matérn covariance structure allows slightly more roughness in the output response than a pure Gaussian function (Rasmussen and Williams, 2006).
An underlying assumption of the Gaussian process emulator is that the output of the cloud model varies smoothly and continuously. Based on this assumption, the emulator fits a smooth response surface that passes directly through each training point. Techniques to allow extra noise in the output response, such as a variance nugget, were explored but were not used in the final emulator design as they did not substantially improve the emulator performance. To test whether the emulator can accurately predict the output of the cloud model, it is necessary to validate the prediction against output from simulations that have not been used to train the emulator. The simulations used to train and validate the emulator are shown in Fig. 5a-c. Fifty-two simulations are used to train each emulator. This is well in excess of the thirty sim-ulations recommended by Loeppky et al. (2009), who states that 10 times the number of variable parameters is required. Eighteen simulations are used to validate the emulator. The model's output from these 18 simulations is compared with the mean and 95 % confidence interval predicted by the emulator at those combinations of the uncertain input parameters.
Variance-based sensitivity analysis is used to measure the sensitivity of the cloud model outputs to the three uncertain input parameters and their interaction effects (Johnson et al., 2015;Saltelli et al., 2000). The overall variance attributed to each input can be separated into the individual or main effect index of each input parameter and the total effect index which comprises the variation attributed to the input parameter in question itself and the variation due to interactions of that parameter with other input parameters (Saltelli et al., 2000). The main effect index of a parameter tells us the proportion of variance in the value of an output variable that could be minimised if the value of the given individual input parameter was known exactly. The difference between the total and main effect indices of a parameter tells us how much variance in the output variable is determined by the input parameter in question interacting with other input parameters (Johnson et al., 2015). In this work, the variancebased sensitivity analysis is carried out using the extended-FAST (Fourier amplitude sensitivity test) approach detailed in Saltelli et al. (1999).

Anvil cloud properties
We first examine the effect of variations in N INP −38 , λ INP , and HM-eff on anvil cloud properties. We focus on the anvil ice properties because the anvil cloud can persist in the atmosphere longer than the deep convective cloud that forms it (and beyond the simulation period presented here) and is therefore climatically important for cloud-radiation interactions. Tropical convectively produced cirrus can persist in the atmosphere for 1-2 d (Luo and Rossow, 2004) while the convective stage of the deep convective cloud simulated here has decayed after approximately 3 h. In Sect. 3.1.1 and 3.1.2, we examine the anvil ICNC and ice crystal size, respectively. An anvil with more numerous, smaller crystals will persist longer in the atmosphere than one with fewer, larger crystals. In Sect. 3.1.3, we examine the simulated anvil cloud fraction and the microphysical properties controlling it. The anvil region of the cloud is defined as the clouds occurring between 150 and 240 min in the simulations with a cloud base height greater than 9 km and an ice water path less than 0.04 kg m −2 . These quantities for specifying anvil cloud were based on qualitatively selecting the anvil region of the deep convective cloud by analysing a number of cloud properties. Other thresholds were tested, e.g. altitudes of 8-11 km, and did not change the results substantially or affect the qualitative findings in any meaningful way.

Anvil ice crystal number concentration
The column-integrated anvil ICNC from all simulations is shown in Fig. 6a-c. Figure 6d shows the associated mean anvil ICNC profile in each simulation. Anvil ICNCs are predominately controlled by the value of N INP −38 (Fig. 6a-d), with a higher N INP −38 causing lower anvil ICNCs at all altitudes ( Fig. 6b and d). This is because the higher the N INP −38 , the higher the rate of heterogeneous freezing at the top of the mixed-phase cloud (Fig. 6e-h), reducing droplet transport to the homogeneous freezing regime and therefore homogeneous freezing rates (Fig. 6i-l). The homogeneous and heterogeneous ice particle production rates shown in  Fig. 5a-d) with the corresponding emulator predictions for anvil ICNC and convective heterogeneous and homogeneous ice crystal number production, along with 95 % confidence intervals on the emulator predictions. All three outputs validate well with points close to or on the 1 : 1 line and small 95 % confidence intervals that overlap the 1 : 1 line most of the time. This indicates that the emulator can capture the variability in the idealised cloud model well for the output variables in question.
Figure 7d-f show the results of variance-based sensitivity analysis and indicate the relative importance of the uncertain input parameters in controlling the variance in the value of the output variable in question. As was inferred from Fig. 6, N INP −38 is the dominant input parameter controlling the variance of anvil ICNC and heterogeneous and homogeneous ice particle production rates, while λ INP and interaction effects contribute a non-negligible but secondary amount to the variance in anvil ICNC. Figure 7d indicates that N INP −38 contributes to over 60 % of this output's variance. This means that the uncertainty in the exact value of the anvil ICNC could be significantly reduced if the value of N INP −38 was to be known exactly. Similarly, this parameter is almost completely controlling the variance in the columnintegrated heterogeneous ice particle production (Fig. 7e), with no real contribution from the other parameters here. In- Figure 6. Anvil ICNC and ice particle production rates. Dependence of anvil ice crystal number concentration (a-d), ice particle production by heterogeneous freezing (e-h), and ice particle production by homogeneous freezing (i-l) on the three uncertain input parameters: λ INP (a, e, and i), N INP −38 (b, f, and j), and HM-eff (c, g, and k). In-cloud profiles of anvil ICNC (d), ice particle production by heterogeneous freezing (h), and ice particle production by homogeneous freezing ( and λ INP is a substantial factor in determining the total number of ice crystals in a convectively generated anvil and may therefore affect the cloud lifetime. Figure 8a, d, and g show the emulator surfaces for homogeneous (Fig. 8a) and heterogeneous (Fig. 8d) ice particle production and anvil ICNC (Fig. 8g) at a fixed HMeff of 350 mg −1 . We hold the HM-eff constant because it had a minimal effect on the variance in the output variables ( Fig. 7d-f); therefore, variations in its value do not alter the shape of the emulated surface substantially. It is important to note that the emulator response surface passes through each simulation point exactly, so it does not allow for noise on each point caused by internal variability of the cloud. As a result, the emulator surfaces should be interpreted by examining the general smoothly varying trends rather than individual bumps which may be an artefact of the emulator representing non-deterministic variations across the parameter space. Methods of smoothing the emulator surfaces could be explored in future studies (e.g. Marshall et al., 2019).
Ice particle production from homogeneous freezing is high and relatively constant between an N INP −38 of 10 −4 and 1 cm −3 before decreasing rapidly at higher N INP −38 values (Fig. 8a). Meanwhile, N INP −38 has the opposite effect on heterogeneous freezing, with heterogeneous ice particle production increasing relatively uniformly with increasing N INP −38 . This is because as more cloud liquid is consumed by heterogeneous freezing at mixed-phase levels due to droplet freezing and the associated increase in secondary ice production, Figure 7. Emulator validation and uncertain input contributions to output uncertainty. Validation of emulator results (a-c) and results of the variance-based sensitivity analysis (d-f) for anvil ICNC (a, b), ice particle production by heterogeneous freezing (c, d), and ice particle production by homogeneous freezing (e, f). In panels (a)-(c), the dots show the value of the validation run on the x axis and the corresponding emulator mean prediction on the y axis. The 95 % confidence intervals on the emulator predictions are also shown. An emulator that validates well will have dots close to the 1 : 1 line and small error bars. Panels riming, and deposition, fewer cloud droplets are available for homogeneous freezing.
Interaction between λ INP and N INP −38 freezing can be seen in the emulator response surfaces. At low N INP −38 , the heterogeneous ice particle production is highest for shallowλ INP values, while at high N INP −38 , the heterogeneous ice particle production rates are highest at steep λ INP values (Fig. 8d). This is because at low N INP −38 values, hetero-geneous freezing at warm temperatures does not limit the number of cloud droplets reaching the upper mixed-phase region. However, at high N INP −38 , a shallow λ INP inducing substantial freezing at warm temperatures will cause substantial cloud liquid consumption (by droplet freezing, secondary ice production, riming, and deposition) that will limit the availability of cloud droplets for heterogeneous freezing at colder mixed-phase levels. In Fig. 8g, we can see that when the N INP −38 is high, the highest anvil ICNCs occur when the λ INP is steep (between −0.3 and −0.5 • C −1 ). This is because at high-N INP −38 values, homogeneous freezing is very low and upper level mixed-phase heterogeneous freezing controls the anvil ICNC. At steep λ INP values the consumption of cloud liquid at warm temperatures is lowest, leading to higher overall rates of heterogeneous freezing. values of 1 and 100 as heterogeneous freezing approaches becoming, and subsequently becomes, the dominant mechanism for primary ice production. Homogeneous freezing is the dominant mechanism of ice crystal production at N INP −38 < 10 cm −3 at all λ INP values (Fig. 8b), above which heterogeneous freezing becomes the dominant mechanism of ice particle production (Fig. 8e). Homogeneous freezing is essentially completely shut off at N INP −38 > 100 cm −3 (Fig. 8b), meaning that at very high N INP −38 all primary ice crystals in the simulated deep convective cloud are formed via heterogeneous freezing. This is because heterogeneous freezing and subsequent processes in the mixed-phase region of the cloud significantly reduce the amount of cloud liquid reaching the homogeneous freezing altitude.
Anvil ICNC decreases sharply as INP concentration increases when heterogeneous freezing approaches becoming the dominant mechanism of primary cloud ice production (Fig. 8h) (> 1 cm −3 ), the highest anvil ICNCs occur at steep λ INP values, which allow more cloud droplets to reach the upper mixed-phase temperatures or the homogeneous freezing regime. Figure 8b, e, and h indicate that an INP number concentration of 1 cm −3 or more is enough to allow heterogeneous freezing to begin to compete with homogeneous freezing, while an INP number concentration of 100 cm −3 will shut off homogeneous freezing completely, regardless of λ INP .   (Fig. 8c). In particular, homogeneous ice particle production declines linearly with increasing λ INP at a N INP −38 of 10 cm −3 . At an N INP −38 of 100 cm −3 , homogeneous freezing is completely shut down, and therefore there is no sensitivity to λ INP evident in the emulator surface. Meanwhile, at a N INP −38 of 10 −3 cm −3 , heterogeneous ice particle production is insufficient at all λ INP values to affect homogeneous ice particle production, which remains uniformly high across the parameter space. Ice particle production by heterogeneous freezing is insensitive to changing λ INP values except at the extremes of the N INP −38 perturbations. There is a slight increase in heterogeneous ice particle production with increasing λ INP at low N INP −38 values and a slight decrease in heterogeneous ice particle production with decreasing λ INP at high N INP −38 (Fig. 8f). Anvil ICNC is sensitive to λ INP values at high N INP −38 (> 10 cm −3 ) where anvil ICNCs decline as λ INP becomes more shallow (Fig. 8i)  −38 values. The anvil ICNC is reduced substantially when the number of heterogeneously frozen ice crystals exceeds the number of homogeneously frozen ice crystals due to the efficient consumption of liquid at upper mixed-phase cloud levels before droplets can be frozen homogeneously. The emulator response surfaces shown in Fig. 8 highlight the complex interactions between heterogeneous and homogeneous freezing and between heterogeneous freezing at different mixedphase temperature levels (determined by the interaction between λ INP and N INP −38 ).

Anvil ice crystal size
We now examine how N INP −38 , λ INP , and HM-eff affect the anvil ice crystal size. The measure we use to quantify ice crystal size is the effective radius or the ratio of the third to the second moments of the ice crystal size distribution. A larger ice crystal size indicates that anvil ice particles will have lower scattering potential (much like the Twomey effect for cloud droplets). A larger ice crystal size also indicates higher fall speeds and lower lifetimes, theoretically reducing the lifetime of the anvil cloud and reducing its radiative effect. The simulated ice crystal effective radius in the anvil cloud region at 14 km can be seen in Fig. 9a-c. We used the effective radius at 14 km because 14 km is the altitude of peak ICNC shown in Fig. 6d.
Anvil ice crystal effective radius exhibits two distinct regimes depending on the value of λ INP , which can be seen in Fig. 9a and b. Simulations with a λ INP shallower (larger) than approximately −0.3 • C −1 (Fig. 9a) exhibit a large jump from under 25 µm and very little variation between simulations to over 27 µm, with a large spread in ice crystal effective radius between simulations. In simulations with a shallow λ INP and an effective radius greater than 27 µm, the value of the effective radius is dependent on the N INP −38 , with simulations with larger N INP −38 values having a larger ice crystal size (Fig. 9b). This indicates that while anvil ICNC was determined predominately by N INP −38 with λ INP having a secondary role, ice crystal size is determined predominately by λ INP with N INP −38 having a secondary role. This is because ice crystal size is more strongly affected than ICNC by the altitude of ice formation, the amount of liquid available for riming, and the time available for deposition growth. Therefore, ice crystal effective radius is predominately affected by INP number concentration at warm temperatures where liquid is available for ice crystal growth which is determined by λ INP .
The mechanism for the increased ice crystal size at shallow-λ INP and high-N INP −38 values is as follows: ice crystals in clouds with a shallower λ INP values have larger concentrations of heterogeneously frozen crystals at warm mixed-phase temperatures (Fig. 9d-f). This increase in heterogeneously frozen ice crystals in the Hallett-Mossop re-gion leads to an increase in ice particle production by the Hallett-Mossop process (Fig. 9g-i). We see a large increase of approximately 1 order of magnitude in ice particle production by the Hallett-Mossop process at shallow λ INP (Fig. 9g) and a bifurcation in cloud behaviour because of this enhancement. The output data are split into two populations, or regimes, based on the λ INP value, with ice effective radius in each regime having a linear dependence on HM-eff (Fig. 9i). Within the warmer temperature mixed-phase cloud region, liquid is still available when crystals are frozen for riming. Therefore, with more heterogeneously frozen ice crystals at lower cloud altitudes, there are higher riming rates (Fig. 9jl), more ice crystal growth, and overall larger ice crystal sizes. Figure 9a, g, and j illustrate a regime change at shallowλ INP values with large increases in anvil ice crystal size (Fig. 9a), Hallett-Mossop ice particle production (Fig. 9g), and accretion of water by ice (Fig. 9j) at values of λ INP above approximately −0.3 • C −1 . This regime change is further illustrated in Fig. 10, which shows the variation in anvil ice crystal effective radius (Fig. 10a), convective Hallett-Mossop ice particle production (Fig. 10b), and accretion of water by ice (Fig. 10c)  −38 greater than 10 cm −3 , the regime change occurs when λ INP is greater than −0.3 • C −1 . The regime change occurs in the same location of parameter space in all three variables (Fig. 10). Simulations in the shallow-λ INP regime with a HM-eff above 600 mg −1 are highlighted with a red outline, and the lack of distinction in colour between simulations with a high HM-eff in the low N INP −38 and steep λ INP regions indicates that a high HM-eff does not have the same effect in the cloud as a shallow λ INP ; i.e. simulations with a steep λ INP and high HM-eff cannot experience the same elevated Hallett-Mossop ice particle production, accretion rates, and resultant increase in the anvil ice crystal effective radius as a simulation with a shallow λ INP and a low HM-eff. However, simulations on the border of the regime transition seem more likely to have elevated ice effective radius and thus be in the shallow-λ INP regime if they have a high HM-eff.
Statistical emulation of anvil ice crystal effective radius at 14 km, Hallett-Mossop ice particle production, and accretion of water by ice crystals was attempted. Figure 11a-c show the validation of the emulator surface against the cloud model validation points. In all three cases the emulator does not validate as well as was seen in Fig. 7 with larger 95 % confidence intervals. Applying a nugget, a term to introduce noise, to allow the emulator to pass nearby to, rather than directly through, the training points (Johnson et al., 2011) was tested as a means to improve the validation. However, Figure 9. Anvil ice crystal size and driving processes. Dependence of anvil ice crystal effective radius (a-c), ice particle production by heterogeneous freezing between 5 and 7.5 km altitude (d-f), ice particle production by the Hallett-Mossop process (g-i), and the accretion of water by ice crystals (j-l) on the three uncertain input parameters: λ INP (a, d, g, j),  Figure 6 shows the total column-integrated heterogeneous ice particle production, while Fig. 9 (here) shows only the heterogeneous ice particle production occurring in the Hallett-Mossop region (5-7.5 km). because the poorer validation occurs mainly as a result of the emulator struggling with the sharp transitions at shallowλ INP values seen in Fig. 9a, g, and j, a nugget term did not change the results. Nevertheless, in most cases the points are relatively close to the 1 : 1 line, indicating that the emulator has some skill in predicting ice crystal size and the cloud development properties that control ice crystal size. Figure 11d-f show the results of variance-based sensitivity analysis and indicate that for all three output variables considered here, λ INP accounts for a large proportion of the variance with a main effect index of 30 % to 60 %. Interaction Figure 10. Regime change in anvil ice crystal effective radius and driving processes. Variation in anvil ice crystal effective radius (a), ice particle production by the Hallett-Mossop process (b), and the accretion of water by ice crystals (c) due to variation in λ INP and N INP −38 . Marker colours indicate the value of anvil ice crystal effective radius (a), ice particle production by the Hallett-Mossop process (b), and the accretion of water by ice crystals (c). Circular markers indicate an ice crystal effective radius above 25 µm (a), an ice particle production by Hallett-Mossop exceeding 2 × 10 4 m −2 s −1 (b), and a rate of water accretion by ice over 1 × 10 −5.5 kg m −2 s −1 (c). Simulations with a HM-eff above 600 splinters mg −1 are indicated with a red outline. Panel Appendix Fig. A1 shows emulator response surfaces for anvil ice crystal effective radius at 14 km (Fig. A1a), ice particle production by the Hallett-Mossop process (Fig. A1b), and accretion of water by ice crystals (Fig. A1c). In Fig. A1a and c, the Hallett-Mossop splinter production efficiency is held constant at 350 splinters produced per milligram of rimed liquid. In Fig. A1b, N INP −38 is held constant at 1 cm −3 . The emulator response surfaces are noisier, with more bumps than those shown in Fig. 8. This is expected due to the larger 95 % confidence intervals on the emulator pre- Figure 11. Emulator validation and uncertain input contributions to output uncertainty. Validation of emulator results (a-c) and results of the variance-based sensitivity analysis (d-f) for anvil effective radius at 14 km (a, b), ice particle production by the Hallett-Mossop process (c, d), and water accretion by ice (e, f). In panels (a)-(c) the dots show the value of the validation run on the x axis and the corresponding emulator mean prediction on the y axis. The 95 % confidence intervals of the emulator mean predictions are also shown. An emulator that validates well will have dots close to the 1 : 1 line and small error bars. Panels dictions shown in Fig. 11a-c. Emulation using a Gaussian process assumes that the uncertain input parameters cause changes in output variables that vary smoothly over the parameter space. This is not the case for the three variables emulated in Fig. 12. For example, ice particle production by the  Figure 12 shows the dependence of anvil cloud fraction (Fig. 12a-d), in-cloud updraught speed (Fig. 12e-h), and total cloud condensate amount (Fig. 12i-l) on the uncertain input parameters. The mean cloud fraction profile occurring between 180 and 240 min of the simulations is shown in Fig. 12d. The anvil cloud fraction values shown in Fig. 12ac, and used in all further analysis of the anvil cloud fraction, are the peak values of the profile shown in Fig. 12d. A similar regime shift at shallow-λ INP values as was seen in the anvil ice crystal size is seen in all three of the output variables considered here (Fig. 12a, e,  −38 region of parameter space (Fig. 13) as was anvil ice crystal size, Hallett-Mossop ice particle production, and ice accretion rates (Fig. 10).

Anvil cloud fraction
Anvil cloud fraction is enhanced at shallow-λ INP values due to an invigoration effect caused by enhanced heterogeneous ( Fig. 9d-f) and secondary freezing (Fig. 9g-i) and increased riming (Fig. 9j-l) in the mixed-phase cloud region, as well as the resultant enhancement in latent heat release, updraught speeds ( Fig. 12e and h), and vertical condensate mass transport ( Fig. 12i and l). Increased ice crystal sizes at shallow-λ INP values in the convectively generated anvil discussed in Sect. 3.1.2 (Figs. 9-12) would be expected to reduce anvil size, due to the associated increases in ice fall speed. Within the time period analysed here, the enhancement in convective strength (inferred from enhanced updraught speeds) and the resultant increase in anvil size at shallow-λ INP values compensate for the effect of increased anvil ice crystal size. The importance of the anvil ice crystal properties relative to the convective invigoration effect for anvil cloud fraction may change with a longer simulation period owing to the persistence of the anvil cloud after the decay of the convection that forms it, and this should be examined in future studies.
The small reduction of anvil cloud fraction within the shallow-λ INP regime with increasing N INP −38 (Fig. 12b) can be attributed to the changes in anvil ice crystal properties reported in Sect. 3.1.1 and 3.1.2

. At high-N INP
−38 val-ues, ICNC is reduced (Fig. 6b) and ice crystal size is increased (Fig. 9b). Fewer and larger crystals will sediment out faster and therefore will spread out over a smaller horizontal area, reducing anvil fraction in simulations with high-N INP −38 values. The chosen Hallett-Mossop splinter production efficiency has very little impact on anvil cloud fraction (Fig. 12c), updraught speeds (Fig. 12g), or cloud condensate amount (Fig. 12k).
Emulation of anvil cloud fraction was attempted but, the bifurcation of these output data into the two distinct regimes depending on the value of λ INP proved impossible to capture with our emulator approach, and validation of the emulation showed little predictive power (not shown). Recently developed emulator approaches that attempt to overcome the smoothness assumption of the Gaussian process emulator used here could be explored in future studies (Pope et al., 2021;Volodina and Williamson, 2020). This indicates that although emulation is a powerful tool to aid in our understanding of cloud processes, traditional methods of analysis are also still needed where there are complex sharp transitions such as those seen in Fig. 13. It is not clear why the emulation of some variables with two distinct regimes (such as ice crystal effective radius) worked relatively well and emulation of anvil cloud fraction did not.

The importance of the Hallett-Mossop process and its interaction with λ INP
One notable feature of the results presented so far is the apparent insensitivity of most cloud properties to the HMeff. For example, the results of the variance-based sensitivity analysis shown in Figs. 7 and 11 indicate that the HM-eff makes no significant contribution to the uncertainty in the value of anvil ICNC, heterogeneous or homogeneous freezing rates, anvil effective radius, or ice accretion of water. Ice particle production by the Hallett-Mossop process was the only output variable shown to have a notable dependence on the HM-eff, and up to 40 % of the uncertainty in its value was attributed to variation in the λ INP value owing to the role of λ INP in determining the regime shift evident in Figs. 9g and 10b. This regime shift induces an enhancement in the ice particle production by the Hallett-Mossop process of about 1 order of magnitude at shallow-λ INP values regardless of the value of HM-eff by increasing the number of primary ice crystals available to initiate the Hallett-Mossop process. In most simulations, over 99 % of ice crystals in the Hallett-Mossop region (5-7.5 km) are formed via the Hallett-Mossop process and not via heterogeneous ice formation (Fig. A2). Figure A2 shows that only 7 of 73 simulations have more than 10 % of the ice particle production between 5 and 7.5 km occurring via heterogeneous ice nucleation rather than via the Hallett-Mossop process. Many output variables, particularly those exhibiting a regime shift at shallow λ INP and high N INP −38 , show a strong correlation with ice particle production in the Hallett-Mossop re- Figure 14. Importance of ice production in the Hallett-Mossop ice production regime (5-7.5 km). Dependence of cloud properties on ice particle production in the Hallett-Mossop regime (5-7.5 km) in the convective stage of cloud development (60-180 min). Shown is ICNC at 7 km (a); column cloud droplet number concentration (b); accretion of water by ice (c); graupel mass (d); snow mass (e); in-cloud updraught speed at 7 km (f); cloud condensate from cloud droplets, rain, ice crystals, snow, and graupel (g); anvil ice crystal effective radius at 14 km (h); and anvil cloud fraction (i). The colour of the markers indicates λ INP and the marker size indicates the HM-eff. Panels gion of the cloud (Fig. 14). This is in spite of the apparent unimportance of the chosen HM-eff for the simulated cloud properties detailed in Sect. 3.1. This correlation indicates that the key role of the INP slope in determining cloud properties can be partly attributed to its role in enhancing Hallett-Mossop ice particle production rates (Fig. 9g-i), which dominate ice production in the Hallett-Mossop regime (Appendix Fig. A1). To avoid biasing the correlation analysis to simulations with unrealistically low concentrations of INP in the Hallett-Mossop regime which have very low variability between simulations, the simulations and correlation analysis shown Fig. 14 comprise only simulations from the realistic region of parameter space (Fig. 5).
Ice particle production by the Hallett-Mossop process is greatly enhanced at shallow-λ INP values due to both the larger availability of seed ice crystals and the enhanced riming that accompany these increased ICNCs. This indicates that INP particles can exert strong control over deep convective cloud properties even when heterogeneous freezing is not the dominant mechanism of ice production because they can alter the rate of ice production by SIP mechanisms (Figs. 9g, 14, and A1).

R. E. Hawker et al.: Effects of mixed-phase ice production on convective anvils
In particular, we note that high rates of ice particle production by the Hallett-Mossop process do not occur unless the λ INP is shallow. This is evident from the lack of distinction between simulations with a HM-eff above or below 600 splinters mg −1 in Figs. 10 and 14 (compare simulations shown with and without a red outline). This indicates that a steep λ INP and a high HM-eff cannot have the same effect on the cloud properties as a shallow λ INP regardless of the HM-eff. Furthermore, ICNCs at lower mixed-phase altitudes, regardless of the freezing mechanism in question, can be key determinants of deep convective cloud properties and the properties of the convectively generated anvil (Fig. 14). Figure 14 indicates that as ice production by the Hallett-Mossop process increases due to increased INP number concentrations at shallow-λ INP values, mixed-phase ICNCs (Fig. 14a) are increased and column cloud droplet number concentrations are reduced (Fig. 14b). Due to the enhancement in ICNC in the lower mixed-phase region with shallower λ INP and higher resultant Hallett-Mossop ice particle production, increases are seen in all mixed-phase freezing mechanisms including accretion of water by ice (Fig. 14c), as well as column graupel (Fig. 14d) and snow (Fig. 14e) mass concentrations due to the well-documented enhanced effectiveness of liquid collection by frozen hydrometeors relative to liquid ones (Johnson, 1987;Phillips et al., 2005). Enhanced latent heat release by the increased freezing events from multiple pathways leads to increased updraught speeds (Fig. 14f) and an overall increase in cloud formation (Fig. 14g). Enhanced riming in the mixed-phase region increases anvil ice crystal effective radius (Fig. 14h) as more anvil ice crystals are formed via heterogeneous freezing and are subject to riming than are formed via homogeneous freezing. The increased convective strength also leads to increased anvil cloud fraction (Fig. 14i).
Many studies have tried to establish a threshold concentration of INP where significant Hallett-Mossop ice production begins to occur and affect cloud properties (e.g. Crawford et al., 2012;Huang et al., 2017). Figure (Fig. 15a), the Hallett-Mossop ice particle production (Fig. 15b), the accretion of water by ice (Fig. 16c), the anvil cloud fraction (Fig. 15d), the updraught speed at 7 km (Fig.16e), and the total column condensate (Fig. 15f). We can see that there is little variation in the cloud properties shown below an N −5 INP of 10 −5 L −1 . Above this INP concentration, there is an increase in all properties shown. Simulations with HM-eff values above 600 splinters mg −1 generally show an enhancement in Hallett-Mossop ice particle production (Fig. 15b), accretion (Fig. 15c), invigoration (Fig. 15e), and cloud condensate (Fig. 15f)   a deep convective cloud particularly in the convective stage of development, though this has not translated to an obvious effect on the resultant cloud anvil effective radius (Fig. 15a) or cloud fraction (Fig. 15d). Figure 15b indicates the threshold INP number concentrations required to initiate or invigorate the Hallett-Mossop process. Hallett-Mossop ice particle production is enhanced to above 50 000 m −2 s −1 at INP number concentrations as low as 10 −4 L −1 when the HM-eff is above 600 splinters mg −1 . At HM-eff values below 600 splinters mg −1 , the Hallett-Mossop ice particle production is enhanced to above 50 000 m −2 s −1 only at INP number concentrations above 0.01 L −1 . This indicates that the threshold INP number concentration for initiating and enhancing the Hallett-Mossop process is dependent on the chosen HMeff. For a base case HM-eff of 350 splinters mg −1 , an INP number concentration of 0.01 L −1 may be enough to significantly enhance the Hallett-Mossop process (see simulations without a red outline in Fig. 15b), in agreement with previous studies (Crawford et al., 2012;Huang et al., 2017).

Limitations of this modelling study
The role of λ INP in determining the ice particle production by the Hallett-Mossop process highlights the importance of the interaction of INP number concentrations with SIP mechanisms. The Hallett-Mossop process is the only SIP mechanism included in simulations in this work, but other SIP mechanisms have been identified in convective clouds . We recommend that the effect of these other SIP mechanisms, including those occurring at temperatures below −10 • C such as droplet shattering Lauber et al., 2018) on deep convective clouds, be tested in similar studies in the future.
The simulated cloud is a single idealised case, and as such the results cannot be directly extrapolated to more realistic convective cloud cases, where less idealised triggering mechanisms are at play (Wellmann et al., , 2020 and different clouds in the population can interact (Hawker et al., 2021). It was not feasible to conduct the necessary number of simulations required to study the impact of three uncertain input parameters on a larger more complicated cloud field due to time and cost restrictions. However, the results presented here provide an interesting stepping stone to understanding the interacting effects of INP number concentrations, INP efficiency, and SIP on deep convective anvil properties. We recommend similar studies be undertaken with more realistic cases, including with less idealised triggering mechanisms, in the future.
The chosen uncertain input parameters are just three of a multitude of microphysical parameters that contribute uncertainty to convective cloud processes which should be considered in future work. For example, uncertainty in ice crystal number and mass concentrations were strongly affected by assumed ice crystal shape in simulations of a continental deep convective cloud simulated using the 3D MAC3 model (Johnson et al., 2015). Uncertainties in the riming, sedimentation, and aggregation rates of snow and graupel particles should also be addressed in the future. Uncertainty in environmental conditions that may affect the cloud properties, for example, the size and temperature perturbation value of the warm bubble initiating our deep convective cloud, or the potential temperature profile , have also not been addressed here. Uncertainties in the initial conditions of our simulations have also not been tested and should be explored in the future (Miltenberger et al., 2018a;Miltenberger and Field, 2020). Additionally, MONC-CASIM is configured to be a two-moment scheme in this work and uses multiple ice categories with fixed parameters for bulk physical properties. The representation of ice properties using a continuous spectrum of physical properties could be explored in the future.
In and steep λ INP on the cloud properties may be larger in reality than was found in this idealised study because of this design feature. This feature of our experiment design also means that combinations of INP number concentration and INP parameterisation slope causing unrealistically low INP concentrations at temperatures above −35 • C are very common in our sampling design. To compensate for this we conducted 22 additional simulations for use in the emulator design and 6 additional validation simulations in the realistic region of parameter space (red and black lines in Fig. 5d). Our variance-based sensitivity analysis is conducted over all simulated parameter space, including the unrealistic space shown in Fig. 5.
It should also be borne in mind that INP concentration spectra often do not follow a simple logarithmic relationship.  Sullivan et al., 2018). Hence, the INP concentration spectra in real clouds can be much more complex than those used in this model. The work presented here highlights the importance of improving our capabilities to represent the complexities of INP number concentrations at mixed-phase temperatures.
The ice crystal properties of the convectively generated anvil are analysed and the implications for anvil lifetime and radiative effect hypothesised. However, the short length of our simulations due to computational limitations means we do not examine the full life cycle of the generated anvil. Conducting similar simulations covering a longer time period would address this limitation and is recommended for the future. tive cloud. A schematic of the main effects identified in this study is shown in Fig. 16. Overall, we find that both λ INP and N INP −38 play a role in determining the anvil cloud properties, with the HM-eff being relatively unimportant in determining the anvil cloud properties. Despite this, we find that the interaction of λ INP with HM-eff is important for determining the resultant amount of ice particle production by the Hallett-Mossop process, which in turn has large effects on the cloud development.
Anvil ICNC is strongly reduced at high-N INP −38 values, with the reduction being more pronounced at shallowλ INP values. Conversely, anvil ice crystal size is increased at shallow-λ INP values, with the enhancement being more pronounced at high-N INP −38 values. This is because the lower the altitude of heterogeneous freezing, the more cloud droplets are consumed by riming and depositional growth and the lower the number of droplets that are available for either homogeneous or cold temperature heterogeneous freezing. A shallow λ INP reduces the number of cloud droplets reaching the top of the mixed-phase regime, and a high N INP −38 reduces the number of cloud droplets reaching the homogeneous freezing regime. Consequently, the anvil consists of a smaller number of large heterogeneously frozen crystals in a high-N INP −38 and shallow-λ INP scenario. The ice crystals transported from the heterogeneous freezing regime to the anvil in this case are larger than those that would be frozen heterogeneously at very cold temperatures or homogeneously. Anvil cloud fraction within the time period studies is enhanced at shallow-λ INP values, and this en-hancement is lower at high-N INP −38 values. This is because a shallow λ INP induces an invigoration effect (inferred from higher cloud updraught speeds) due to an increase in Hallett-Mossop ice production and enhanced glaciation at mixedphase temperatures, leading to larger condensate mass divergence in the upper troposphere and a more extensive anvil. The anvil is smaller in a high-N INP −38 scenario due to the reduced ICNC and increased ice crystal size discussed above, which serves to reduce cloud lifetime.
Statistical emulation and variance-based sensitivity analysis allow us to identify complex interdependencies between input and output variables. We find that the interaction between λ INP and N INP −38 account for up to 30 % of the variation in values of anvil ICNC and anvil ice crystal size. The emulator surfaces help us to see variation in the importance of input parameters depending on the value of other uncertain inputs that would not otherwise be apparent in one-ata-time tests. In particular, we identify several important interdependencies between different freezing mechanisms. For example, at high N INP −38 , λ INP is important in determining the anvil ICNC and homogeneous freezing rates because the shallower the λ INP , the fewer cloud droplets available to be frozen homogeneously or heterogeneously at the top of the mixed-phase regime. However, at lower N INP −38 values, λ INP is nearly inconsequential to homogeneous freezing rates and anvil ICNC because the N INP −38 is low enough that a shallow λ INP does not substantially affect the number of droplets reaching upper cloud levels. Furthermore, the dominant effect of λ INP on many cloud properties is in part attributed to the fact that high-λ INP values provide seed crystals for the Hallett-Mossop process, vastly increasing the number concentration of ice crystals between 5 and 7.5 km and subsequently the cloud riming and deposition rates.
The amount of Hallett-Mossop ice particle production is determined by both HM-eff and λ INP , with a λ INP above ∼ −0.3 • C −1 causing a jump of about an order of magnitude in Hallett-Mossop ice particle production while the effect of HM-eff is linear. A regime shift to a cloud with extensive glaciation at warm temperatures, stronger convective updraughts, larger condensate mass divergence in the upper troposphere, and a more extensive anvil occurs for λ INP values between −0.3 and −0.1 • C −1 , with the exact value of the transition depending on the N INP −38 values (the transition occurs at steeper λ INP values if the N INP −38 is high). This regime change is driven by a shallow λ INP (particularly when combined with a high N INP −38 ) forming more ice crystals in the Hallett-Mossop temperature regime and thus seeding a stronger ice particle production by the Hallett-Mossop process. We find a strong enhancement in Hallett-Mossop ice particle production occurs at an INP number concentration anywhere between 10 −4 to 1 L −1 depending on the HMeff. This indicates that the threshold INP number concentration needed to cause substantial SIP varies depending on the strength of the SIP mechanism in question. Whether the interaction of INP with other SIP mechanisms that operate over different temperature ranges has a similar effect on SIP production rates should be explored in the future.
Tropical cirrus can typically persist, and therefore affect radiation, in the atmosphere for days after the convective cloud that formed them has decayed (e.g. Luo and Rossow, 2004). While the simulations were not long enough to study the full life cycle of the anvil, we can infer possible implications of the uncertain input parameters on anvil lifetime based on the anvil ice crystal properties. An anvil with more numerous, smaller crystals will persist longer in the atmosphere than one with fewer, larger crystals. High-N INP −38 values lead to both fewer and larger anvil ice crystals and a slightly reduced cloud fraction, implying that deep convection generated in an environment with a high concentration of ice-nucleating aerosol, such as in high dust loading, may have a shorter lifetime than those generated in lower N INP −38 environments. At shallow-λ INP values, the anvil ice crystals are larger, indicating that a cloud anvil formed in the presence of an INP population with high efficiency at warmer temperatures, e.g. marine organics (Wilson et al., 2015), may have a shorter lifetime than that formed in the presence of an INP population with a steeper temperature dependence. However, this effect is compensated for by an invigoration effect (inferred from higher in-cloud vertical velocities) driven by higher rates of mixed-phase ice formation in simulations with shallow-λ INP values, which leads to larger condensate mass divergence in the upper troposphere and a larger anvil in the timescale simulated in this study. Future studies should cover the entire life cycle of the generated anvil cirrus in order to quantify the resultant lifetime of the convective anvil due to compensation between these two effects of λ INP and N INP −38 . We demonstrate with the present study that statistical emulation is a powerful tool for visualising and quantifying the relationships between cloud responses and different uncertain parameters (Fig. 8). However, statistical emulation struggles to accurately represent cloud responses where there is a significant regime shift at shallow λ INP (Figs. 10 and 14), for example, the anvil cloud fraction. We therefore suggest that emulation be used alongside traditional analysis methods for the further study of the complex processes occurring within deep convective clouds, particularly where sudden transitions or regime shifts are evident or likely. The use of Latin hypercube sampling to capture cloud responses over the full realistic parameter space to multiple uncertain input parameters is very effective, even without undertaking statistical emulation of the simulation data.
The microphysical effects of the variations in INP number concentrations and INP parameterisation slope detailed here build on the results of Hawker et al. (2021) and further our understanding of the role of these two uncertain inputs on deep convection. In both the complex cloud-field simulation of Hawker et al. (2021) and the idealised deep convective cloud presented here, INPs in the warm mixed-phase region enhance Hallett-Mossop ice particle production and increase snow and graupel formation, leading to an invigoration effect, more cloud condensate, and an increased cloud fraction at mixed-phase cloud levels. INPs in the mixed-phase region also reduce homogeneous ice production, leading to reduced overall column-integrated ICNCs in both studies. Conversely, in the Hawker et al. (2021) study, shallow-λ INP values led to a reduced cloud fraction above 9 km due to reduced ice particle production by homogeneous freezing. In the deep convective cloud simulated here, a shallow λ INP leads to an increased anvil cloud fraction due to an invigoration effect caused by enhanced ice formation and latent heat release in the mixed-phase cloud region. This indicates that the microphysical effects of INPs and the interaction of INPs with the Hallett-Mossop process are relatively consistent between realistic and idealised case studies. However, the consequences of microphysical changes due to INPs on the cloud macro-physical properties such as cloud fraction and outgoing radiation can be different depending on the specific conditions of the simulation and the type of cloud being perturbed. For example, the deep convective cloud simulated here was initiated with a relatively strong warm bubble perhaps predisposing the cloud fraction to be more sensitive to enhancements in an already strong convective updraught strength than the more realistic clouds in the Hawker et al. (2021) study. This study has enhanced our understanding of aspects of cloud microphysical behaviours that was not fully explained in Hawker et al. (2021). For example, in both studies a reduction in homogeneous freezing rates occurs where there are high INP number concentrations at low temperatures, and the results from this study indicate that it is a high N INP −38 rather than a steep λ INP that is the main driver of this effect, a distinction that we were unable to make in the previous study. The differences in resolution between the two studies (250 m (here) and 1 km; Hawker et al., 2021) may also cause divergences in the microphysical responses of clouds to perturbations (Varble et al., 2020).
This work highlights the complexity of interactions between mixed-phase ice processes and the challenge of representing them accurately in numerical weather prediction models. Our work indicates that the sensitivity of deep convective cloud properties to mixed-phase ice processes varies depending on ambient ice-nucleating aerosol concentrations (e.g. absolute dust concentrations) as well as the efficiency of the available ice-nucleating aerosol (e.g. whether the INP number concentration consists of dust or marine organic particles). The potential for ice particle production by INPs and SIP to impact anvil cirrus ice properties also presents a challenge for climate models. Climate models do not typically use INP number concentrations to determine ice water path and the resultant outgoing radiation, and this is an important area for future work (Baran et al., 2014;Waliser et al., 2009). The role of the temperature dependence of INP number concentration in determining the observed cloud properties indicates the importance of quantifying the concentration of INP at all mixed-phase temperatures, adding to work by, for ex-ample, Hawker et al. (2021), Liu et al. (2018), Shi and Liu (2019), and Takeishi and Storelvmo (2018). Furthermore, the temperature dependence of the INP parameterisation had a substantial effect on Hallett-Mossop ice particle production rates, indicating that heterogeneous freezing can be an important determinant of deep convective cloud properties by affecting SIP mechanisms, even when heterogeneous freezing is not the dominant mechanism of ice formation in the SIP region.  Simulations with a HM-eff above 600 splinters mg −1 are indicated with a red outline. Ice particle production rates are calculated from the convective stage of the simulation (between 60 and 180 min).
Data availability. The datasets generated and analysed in this study are available from the corresponding author on reasonable request.
Author contributions. REH, AKM, BJM, JSJ, PRF, and KSC contributed to the design, development, and direction of the study. AKM and REH set up a default deep convective simulation in the MONC-CASIM model. JSJ provided the base R code needed for the uncertain input parameter combination selection and to carry out the statistical emulation and uncertainty analysis, and they provided advice about statistical emulation and uncertainty analysis throughout. REH carried out the model development in the MONC-CASIM code to build the base case simulation and allow for the perturbations to the uncertain input parameters determined by the sampling design, used and modified the R code to select the uncertain input parameter combinations, ran all MONC-CASIM simulations presented here, conducted all analysis, and wrote the manuscript. JMW, AAH, and BJS built and maintained the Met Office MONC-CASIM model used to run the simulations. Everyone provided comments and edits to the manuscript.
Competing interests. At least one of the (co-)authors is a member of the editorial board of Atmospheric Chemistry and Physics. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Acknowledgements. We acknowledge the use of MONSooN, a collaborative high-performance computing facility funded by the Met Office and NERC. We acknowledge the use of JASMIN, the UK collaborative data analysis facility. The simulation initial conditions are based on the "Ice in Clouds Experiment -Dust" (ICE-D) field campaign, which was also funded by the Natural Environment Research Council (NERC, grant NE/M00340X/1). The ICE-D campaign used the BAe-146-301 atmospheric research aircraft, which is operated by Directflight Ltd (now Airtask) and managed by the Facility for Airborne Atmospheric Measurements (FAAM). At the time of the measurements FAAM was a joint entity of NERC and the UK Met Office. We thank all the people involved in the ICE-D campaign.
Financial support. This research has been supported by the European Research Council, H2020 (MarineIce, grant no. 648661), and the Natural Environment Research Council (grant no. NE/M00340X/1).
Review statement. This paper was edited by Jianzhong Ma and reviewed by Xiaohong Liu and one anonymous referee.