Impact of Holuhraun volcano aerosols on clouds in cloud-system resolving simulations

. Increased anthropogenic aerosols result in an enhancement in cloud droplet number concentration ( N d ), which consequently modiﬁes the cloud and precipitation process. It is unclear how exactly cloud liquid water path (LWP) and cloud fraction respond to aerosol perturbations. A volcanic eruption may help to better understand and quantify the cloud response to external perturbations, with a focus on the short-term cloud adjustments. The goal of the present study is to understand and quantify the response of clouds to a selected volcanic eruption and to thereby advance the fundamental understanding of the 5 cloud response to external forcing. In this study we used the ICON (ICOsahedral Non-hydrostatic) model in its numerical weather prediction setup at a cloud-system-resolving resolution of 2.5 km horizontally, to simulate the region around the Holuhraun volcano for one week (1 – 7 September 2014). A pair of simulations, with and without the volcanic aerosol plume, allowed us to assess the simulated effective radiative forcing and its mechanisms, as well as its impact on adjustments of LWP and cloud fraction to the perturbations of N d . In comparison to MODIS (Moderate Resolution Imaging Spectroradiometer) 10 satellite retrievals, a clear enhancement of N d due to the volcanic aerosol is detected and attributed. In contrast, no changes in either LWP or cloud fraction could be attributed. The on average almost unchanged LWP is a result of some LWP enhancement for thick, and a decrease for thin clouds.

numerical weather prediction (NWP) variant is used (ICON-NWP). The resolution corresponds to approximately 2.5 km in the horizontal (R2B10 triangular grid). In the vertical, 75 layers with top height at 30 km were chosen.
The physics package of ICON-NWP includes a comprehensive double moment cloud liquid and ice microphysical scheme 60 (Seifert and Beheng, 2006). Because of using a rather fine resolution, deep convection is considered to be resolved, whereas, for shallow convection, the parameterization scheme by Tiedtke (1989) with modifications by Bechtold et al. (2008) was used. To achieve a more realistic representation of the Twomey effect, we furthermore coupled the hydrometeor number concentrations from the double moment microphysical scheme to the radiation scheme as proposed in Kretzschmar et al. (2020).
Initial and boundary conditions were derived from the European Centre for Medium-Range Weather Forecast (ECMWF) 65 Integrated Forecasting System (IFS) operational analysis. The 2014 Holuhraun eruption was a fissure eruption that started on 20 August 2014 and ended on 25 February 2015. By 7 September 2014, the lava field had extended more than 11 km to the north (Kolzenburg et al., 2017). We choose the period from 1 to 7 September 2014 for our analysis because the lava field had developed sufficiently in this period and substantial amounts of SO 2 had been emitted into the atmosphere, while, at the same time, a well-defined plume is observable. An additional feature in simulations that must be mentioned, is the 70 implementation of a satellite simulator into the model. Satellites are essential tools to assess the character of clouds due to their global coverage and availability (Lai et al., 2019). Differences between model simulations and satellite retrievals stem in 3 https://doi.org/10.5194/acp-2022-38 Preprint. Discussion started: 24 January 2022 c Author(s) 2022. CC BY 4.0 License. part from a different definition of the respective quantities that are compared. Therefore, one approach to reduce inconstancy between model simulations and satellite retrievals is to use satellite simulators in models to mimic the observational processes (Roh et al., 2020). The COSP satellite simulator (Bodas-Salcedo et al., 2011) is an open source work package developed by 75 CFMIP (Cloud Feedback Model Intercomparison Project) to replicate active and passive satellite data using variables from the model as an input (Webb et al., 2017). In this study, just satellite simulator for MODIS (Moderate Resolution Imaging Spectroradiometer) observations (Pincus et al., 2012) was used. The COSP simulator uses several model variables as input such as temperature, pressure, cloud fraction and cloud water content (Kretzschmar et al., 2019) to generate what the MODIS retrievals would capture given the simulated clouds fields (Saponaro et al., 2020).

80
In the cloud-resolving simulation (each grid box is either fully cloudy or clear), the use of sub-grid variability, one of the features of COSP for application in general circulation models, was not necessary. In order to evaluate COSP related variables in our simulations, the collection 6.1 Level-2 MODIS-Aqua optical and physical cloud data product was used ; therefore, swaths with 1 km spatial resolution for r e , cloud optical thickness (τ c ) and LWP were used and remapped to the model resolution to have an accurate comparison. Furthermore, the planetary albedo at the top of the atmosphere is 85 analyzed as retrieved by the Clouds and the Earth's Radiant Energy System (CERES) instrument onboard the Aqua satellite (Su et al., 2015;Loeb et al., 2016).

CCN activation
The ICON-NWP version applied in this study does not contain an interactive aerosol model; therefore, in this section, we discuss how CCN are activated into clouds droplets in the default model setup and afterward we introduce a new method for 90 CCN activation in microphysical scheme, which had specifically been developed for this study. In the default setup of ICON, CCN activation uses a parameterization that employs a functional dependency of grid scale vertical velocity and pressure to derive the number of newly activated CCN (Hande et al., 2016). Hande et al. (2016) performed model simulations that considered a multi-modal interactive aerosol scheme to provide information on the formation and transport of aerosols in Europe and, by using the parameterization of Abdul-Razzak and Ghan (2000, ARG), derived CCN number concentrations for 95 different vertical velocities for a selected date (30 April 2013). This parameterization thus assumes a temporally and spatially constant profile of CCN which is representative for CCN background over Europe. For that reason, this parameterization alone can not provide information about CCN concentration within a plume of volcanic aerosol.
In order to more accurately represent the aerosol plume, we use look-up tables that contain the number of activated CCN as a function of pressure p and vertical velocity w as an input for the ICON simulation. The number of activated CCN is interpolated to obtain the information about the spatial-temporal distribution of the aerosol mass mixing ratio by aerosol species. The

105
CAMS reanalysis provides aerosol mass mixing ratios at 60 hybrid sigma/ pressure levels up to 0.1 hPa, and covers the 2003 to calculates the number of activated aerosols employing the Köhler theory (Köhler, 1936), we created look-up tables of activated CCN for our simulation.
In our study, the ARG-parameterization is employed offline, by running a box model setup and using aerosol mass mix-110 ing ratio from the CAMS reanalyses as an input for various vertical velocities. The ARG-parameterization has been used in microphysical schemes in a wide range of resolutions before, ranging from global climate models to cloud resolving models (Ghan et al., 2011;West et al., 2014;Luo et al., 2008). The ARG parameterization is based on the competition between aerosol particles for available water vapor which depends on aerosol particle composition, size distribution and most importantly the supersaturation forcing rate obtained by the updraft. We evaluate supersaturation S 0,i at ten specific values of vertical velocity 115 used in the look-up tables (see Costa-Surós et al., 2020). After calculating S max , the critical radius of activation for each aerosol mode is obtained in the box model. When the supersaturation for each aerosol mode is smaller than maximum supersaturation S max ≥ S 0,i , the environment has gained the needed supersaturation to activate the particles. Using this approach, an observations-tied spatially-temporally varying input number concentration of activated CCN for ten prescribed vertical velocity classes was produced. In the CAMS reanalyses data, the aerosols emitted from the Holuhraun volcano are not represented: 120 it is firstly not constrained by the data assimilation. CAMS assimilates MODIS aerosol optical depth (AOD) retrievals (Levy et al., 2013), but due to the presence of extensive clouds in the region of interest, MODIS was not able to capture sufficient information about AOD. Secondly, it is also not included in the model simulation, because in the emission source model of CAMS, no volcanic emissions are considered. Therefore, the CAMS data was used to obtain background spatial and temporal aerosols concentration and in order to implement aerosol concentrations inside the plume, the sulfate aerosol concentration 125 in CAMS was scaled based on the SO 2 emission monitored by Ozone Mapping and Profile Suite (OMPS) satellite retrievals which will be explained in more detail in the next session (Yang, 2017).

The volcanic-aerosol plume in the model simulations
Lava flows and emitted gases from volcanic eruptions are the most common features that remotely can be monitored globally and at different time scales. SO 2 is one of the most common gases emitted from volcanic eruptions and can be retrieved by 130 spaced-based sensors (Fioletov et al., 2020). In this study, the OMPS data product (Level 2) for SO 2 was used. This data set provides information about vertically integrated SO 2 (in Dobson units, DU). The SO 2 retrievals for 1 to 7 September 2014 for the lower troposphere are shown in Figure 2.
The SO 2 plume was detected on 1 September shortly after the beginning of the eruption and evolved over time mostly eastwards, towards Scandinavia. Former studies compared OMPS satellite retrievals with surface observations for the Holuhraun 135 eruption and showed that satellite retrievals are able to detect spatial and temporal evolution of the volcanic plume (Ialongo et al., 2015). In this study, we performed two simulations over the domain shown in Figure 1, one with background aerosol concentrations only, which is referred to as the no-volcano simulation, and one with scaling the sulfate concentrations in the CAMS reanalysis data within the plume as defined by the OMPS satellite retrievals, referred to as the volcano simulation in this article. As shown in Figure 2, grid-points with SO 2 concentrations in the lower troposphere exceeding 1 DU are considered to constitute the plume. For these grid-points, a scale factor field was computed by dividing the SO 2 concentrations retrieved within the plume by the mean SO 2 concentration for the entire domain outside the plume region. In the next step, the sulfate aerosol mass mixing ratio from the CAMS reanalyses was scaled inside of plume by these scaling factors before deriving a new CCN distribution that now considers the volcanic plume with the enhancement consistent with the OMPS satellite retrievals. location of the plume can smoothly be identified. This information lead us to perform two simulations one with a background activated CCN concentration (left panels) referred as no-volcano simulation, and one with scaled activated CCN concentration 150 (right panels) referred to as volcano simulation.

Results
The present study aims at a detection and attribution approach, assessing the differences in cloud properties within and outside the volcanic plume by comparing a factual and a counterfactual simulation with satellite observations. This aims to evaluate how cloud microphysical properties (N d and LWP) behave differently in and outside the volcano plume. (left column in Figure 3). Nevertheless, the grid points that are located inside of the volcano plume are compared to the ones outside the plume to assess differences due to different meteorological conditions. 160 N d is the first microphysical variable we assess. N d is not directly retrieved by the operational MODIS satellite retrievals.
Instead, r e and τ c are retrieved using the method described by Nakajima and King (1990). On the basis of such retrievals, assuming clouds that behave like adiabatic ones, N d can be computed as follows (Grosvenor et al., 2018): In this relation, γ depends mainly on the adiabatic condensation rate and can be approximated as 1.37·10 −5 m − 1 2 (Quaas et al., . In order to obtain N d by Equation 1 in our analyses both in simulations and MODIS, r e less than 4 µm and τ c less than 4 were excluded from data set because they are less reliable (Nakajima and King, 1990). For consistency,N d is derived from the  and an increase in the probability of thicker clouds (with higher LWP) compared to the no-volcano simulation. The MODIS 195 distribution for LWP inside the plume indicates that the probability for shallower clouds is less than what the simulations show, but the probability for thicker clouds is higher than in the no-volcano run, albeit also less than in the volcano run. In terms of the mean values for LWP (Table 1)   average, however, seems to come about by a decrease in LWP for the clouds with low LWP, and an enhancement of LWP for large LWP values ( Figure 5. This is qualitatively consistent with the results of the ICON-NWP model. The model, however, exaggerates the increase in large LWP values, leading to the exaggerated mean increase. The question is now what is the underlying process leading to an increase in LWP in the volcano simulation? One reason is 205 the suppression of precipitation (e.g., Seifert et al., 2012). Therefore, the distribution of rain water path (RWP) was analyzed to investigate the alteration of precipitation inside and outside the volcano plume in both, the volcano and no-volcano simulations.
The comparison is shown in Figure 6. Since the precipitation information is not available from MODIS or other satellite retrievals, RWP is only depicted for the simulations. Inside the volcano plume, there is a decrease in light rain and an increase in heavy rain for the volcano simulation, compared to the no-volcano simulation. In terms of mean values for RWP (Table 1), there 210 is a decrease in the volcano run by 15 % on average, while the precipitation profile for outside of plume is quite similar which is in the agreement of the fact that LWP for outside of plume didn't alter significantly. Moreover, suppression in precipitation can also lead to enhancement in cloud horizontal extent (cloud fraction). Therefore, the modification in cloud fraction was examined in simulations and MODIS. The analyses for mean values of total cloud fraction in Table 1 demonstrates that, in the volcano simulation, cloud fraction is enhanced in the plume compared to outside the plume by 40 %, while the enhancement is 215 only 32 % in the no-volcano simulation. However, even in the no-volcano simulation, cloud fraction inside of plume is higher than outside of plume by 32 % due to the different weather conditions, and this is consistent with what MODIS shows (29 %).
Finally, the effect on radiation (indicative of the effective radiative forcing due to the modification of cloud properties by the volcanic aerosol) is examined. Therefore the TOA albedo was analyzed inside and outside of plume in simulations and CERES level-2 footprint data (Su et al., 2015). For the comparison, the simulation output was remapped to 20 km horizontal resolution to be consistent with the resolution of the CERES footprint. In Figure 7 TOA albedo for the cloudy sky is depicted for inside and outside the volcano plume for both simulations and CERES data. Clear sky was excluded because, in the model, no aerosol-radiation interactions are considered, but in the CERES this effect is in the data and would bias the analysis for clear sky. An additional important aspect that should be considered, is that the TOA albedo distribution is considered here for liquid clouds with τ c more than 4 because in obtaining N d the data with τ c less than 4 were excluded as well. Considering the 225 TOA albedo distribution inside the plume, it is seen that in the volcano simulation, there is a higher probability for TOA albedo larger than 0.6 compared to the no-volcano simulation. In the CERES data, there is a peak at TOA albedo between 0.4 and 0.6 that is not as pronounced in either simulation. In turn, the probability for TOA albedo larger than 0.7 is smaller in the data than in both simulations. This bias, however, is clear outside the plume but much less so inside the plume -possibly indicative of the albedo enhancement due to the volcanic aerosol.

230
For the mean values (Table 1) volcano simulation is 15 % larger compared to no-volcano simulation. In CERES data there is an 18 % enhancement inside the volcano plume compared to outside the plume. When compared to the difference between inside and outside the plume in the no-volcano simulation (27 %), it is difficult to conclude that there is a signal of alteration in TOA albedo in CERES 235 data. We also analyzed cloudy sky TOA albedo mean values in simulations and CERES. The values in Table 1 demonstrate an enhancement of 9 % in CERES and 7 % in volcano simulation while no changes were obtained in no-volcano simulation. The daily mean incoming solar radiation was obtained 260 W m −2 ; therefore, effective radiative forcing except cloud cover effect can be estimated as 10 W m −2 in CERES and 8 W m −2 in volcano simulation.