Detection and attribution of aerosol-cloud interactions in large-domain large-eddy simulations with ICON

Clouds and aerosols contribute the largest uncertainty to current estimates and interpretations of the Earth’s changing energy budget. Here we use a new-generation large-domain large-eddy model, ICON-LEM, to simulate the response of clouds to realistic anthropogenic perturbations in aerosols serving as cloud condensation nuclei (CCN). The novelty compared to previous studies is that (i) the LEM is run in weather prediction mode and with fully interactive land surface over a large domain, and 5 (ii) a large range of data from various sources are used for the detection and attribution. The aerosol perturbation was chosen as peak-aerosol conditions over Europe in 1985, with more than five-fold more sulfate than in 2013. Observational data from various satellite and ground-based remote sensing instruments are used aiming at a detection and attribution of this response. The simulation was run for a selected day (2 May 2013) in which over the selected domain of central Europe a large variety of cloud regimes was present. 10 It first is demonstrated, using satellite aerosol optical depth retrievals available for both 1985 and 2013, that the aerosol fields for the reference conditions and also for the perturbed ones, as well as the difference between the two, were consistent in the model and the satellite retrievals. In comparison to retrievals from ground-based lidar for 2013, CCN profiles for the reference conditions were consistent with the observations, while the ones for the 1985 conditions were not. 1 https://doi.org/10.5194/acp-2019-850 Preprint. Discussion started: 21 October 2019 c © Author(s) 2019. CC BY 4.0 License.

. Temporal sequence of satellite images at (a) 9 UTC, (b) 13 UTC and (c) 17 UTC on 2 May 2013 as natural color composite of the 0.6-, 0.8-, and 1.6-µm channels and the high-resolution visible channel from the Spinning Enhanced Visible and Infrared Imager (SEVIRI) instrument on board the geostationary Meteosat satellite. Liquid clouds appear in shades of white and an increase in ice content shifts the color of the clouds towards cyan. Table 1. Summary of simulations for 2 May 2013. Control ("C") runs are with (low) 2013 aerosol concentrations, Perturbed ("P") runs, with (high) 1985 aerosol. The first set takes into account two moments of the cloud particles size distributions in radiation ("2R"), the specific mass and number, whereas the second one, only one moment ("1R"), which is the specific mass.

CCN for 2013 and 1985
The anthropogenic aerosol emissions in the Northern Hemisphere exhibited a strong increase since the industrialization starting in the early 18th century and accelerating from the 19th century onwards. Over Europe and North America, they reached a maximum in the 1980s, and declined since then (Smith et al., 2011;Cherian et al., 2014). Over Europe, the peak was at about 40 Tg sulfur per year. After that, they decreased to early 1900s levels because of the introduction of air quality policies, and by the year 1985) and the mid-2010s (represented here by the year 2013), a very substantial change in aerosol concentrations occurred (Smith et al., 2011).
A prerequisite to realistically simulate the aerosol-cloud interactions is a realistic representation of the aerosol concentrations and their capacity to serve as CCN. Compared to the earlier model version (Heinze et al., 2017;Hande et al., 2016), new time-varying 4D distributions of CCN concentrations were generated from the emissions valid for 2013 and for the peak-5 aerosol conditions around 1985. The simulation imposing these latter CCN concentrations is hereafter called "perturbed", in comparison to the "control" run which employs CCN concentrations for 2013. The difference between peak-aerosol in 1985 and 2013, rather than, e.g. a comparison of 2013 conditions to pre-industrial aerosol, has been chosen for two reasons: (i) the perturbation is much larger, since the emissions in 2013 were much closer to pre-industrial than to 1985 levels, and (ii) some observations for 1985 are available to assess the aerosol fields. 10 The CCN distributions were created with the regional coupled model system COSMO-MUSCAT, which consists of two online-coupled codes: the chemistry transport model Multi-Scale Chemistry Aerosol Transport (MUSCAT) (Wolke et al., 2004(Wolke et al., , 2012 and the Consortium for Small-scale Modelling (COSMO) model, which is the the operational weather forecast model of the German Meteorological Service (Deutscher Wetterdienst, DWD) and other national meteorological services in Europe (Baldauf et al., 2011;Schättler et al., 2014). The COSMO model, which is a non-hydrostatic meteorological model and solves 15 the governing equations on the basis of a terrain-following grid (Baldauf et al., 2011), provides MUSCAT with all required meteorological fields. Based on these fields, MUSCAT calculates the transport and transformation processes which include advection, turbulent diffusion, and physico-chemical conversion of particles and trace gases in the air, as well as sink processes (sedimentation, dry and wet deposition) (Knoth and Wolke, 1998;Wolke et al., 2012). The emissions of anthropogenic primary particles and precursors of secondary aerosols is prescribed using emission fields from the European Monitoring and Evaluation 20 Programme (EMEP, 2009; http://www.emep.int/ ). The emission fluxes of natural primary aerosols (e.g. desert dust, primary marine particles), are computed online within the model depending on meteorological fields (surface wind speed, precipitation) (Heinold et al., 2011).
For the present study, two periods of the year 2013 were simulated with COSMO-MUSCAT, coinciding with the measurement campaigns of the HOPE experiment (Macke et al., 2017). It includes measurements around JOYCE (Löhnert et al.,  UTC) were analyzed. The ICON-COSP simulations were temporally and spatially matched to the satellite observations as well as regridded to the MODIS data resolution of 1 km. Only cloudy satellite pixels with assigned liquid cloud phase as well as good quality and solar zenith angles below 50 • were considered, to exclude uncertain/problematic cloudy retrievals.
For liquid clouds, cloud droplet number concentration (N d ) is derived from cloud effective radius (r e ) and cloud optical thickness (τ c ) as in Quaas et al. (2006); where α = 1.37·10 -5 m -0.5 : an approach that assumes an adiabatic growth of clouds (Grosvenor et al., 2018). CCN number concentration, r e , liquid water content (q l ), and N d , which are compared with the model output (Section 3.4). The CCN concentration profiles of the lidar measurements were calculated using the method described in Mamouri and Ansmann (2016). Profiles of liquid cloud microphysical properties were derived from ground-based remote sensing using the recently established dual-field-of-view (DFOV) lidar techniques (Grosvenor et al., 2018). Such observations are available at Leipzig, Germany (51.3 • N, 12.4 • E), since 2013 and provide information about aerosol-cloud-interaction processes (Schmidt et al.,5 2014). Originally, the observations were based on the DFOV Raman lidar technique (Schmidt et al., 2013), which can only be applied to nighttime lidar observations in order to reduce effects of the solar background on the measurements of Ramanscattering lidar returns from nitrogen molecules. Progress that was made recently in the accuracy of polarization measurements with lidar (Jimenez et al., 2019a) allows to apply an alternative technique for profiling of liquid-cloud microphysical properties even during daytime (Jimenez et al., 2017(Jimenez et al., , 2019b. In this novel DFOV-polarization approach, the liquid water content, q l , 10 profile is assumed to increase adiabatically with height from the cloud bottom, while the N d remains constant. An inversion scheme exploits a non-ambiguous relation between r e and the extinction coefficient, both dependent on N d and q l , with the single-FOV depolarization ratio and the relative depolarization, quantities that a DFOV polarization lidar can measure (Jimenez et al., 2017(Jimenez et al., , 2019b.

Ground-based
The cloud radar (8.6 mm wavelength) also based in the MOL-RAO, is well suited for the study of thin, low-reflectivity clouds 15 such as non-drizzling and drizzling stratocumulus clouds, due to its high sensitivity. The cloud radar transmits linear polarized Only liquid clouds are selected following Cloudnet criteria (Illingworth et al., 2007). In section 3.5 the ICON-LEM forward simulated reflectivities have been compared to the radar observations for four different categories of drizzle-cloud droplets-rain drops.

Evaluation of the perturbed and unperturbed CCN and aerosol distributions
A comparison of the AOD over the North and Baltic seas as retrieved from AVHRR for 1985 and 2013 and AOD simulated by COSMO-MUSCAT is shown in Fig. 2 and summarized in Table 2. Note that since AVHRR retrieves AOD only in cloudless cases over sea, only a rather small fraction of the time a valid retrieval is available (Table 2) Table 3). In the following sections (3.3 and 3.4) it is explored to which extent this perturbation is also seen by satellites and from ground-based remote sensing in terms of N d retrievals. The total-water 10 mass difference at about 2 km altitude (vertically integrated relative increase by 8.8 %) is mainly due to decreased rain in the perturbed simulation (-12.3 %). In the following sections, the changes in the LWP and liquid water content (q l ) are investigated.
A slightly higher homogeneous cloud droplet freezing for the perturbed simulation is triggered by upward transport of cloud droplets. In turn, graupel number and mass concentrations are higher in the cleaner environment in low to mid altitudes (3 -   (Table 4). A reason could be the MODIS instrument sensitivity, since optically very thin clouds are not observed or give problematic retrievals (Grosvenor et al., 2018), and for optically very thick clouds the measurements can go into saturation (note that N d is computed from the cloud optical thickness and cloud effective radius, see Eq. 1). The simulated N d distribution for the perturbed simulation (P2R) is shifted to significantly higher values. A factor of about 2 between ICON-control and ICON-perturbed is reflected in the percentile values.

5
However, even if the ICON-LEM output is sampled along the MODIS swath, it does not show this distinct behaviour, but a smooth PDF. The difference between LWP from the ICON reference simulation (C2R) and ICON perturbed simulation (P2R) is small compared to the LWP variability, and small in comparison to the model bias with respect to the MODIS retrievals. It is nevertheless a systematic increase in LWP that is simulated, even if it is small as also expected from recent investigations of satellite data (Malavelle et al., 2017;Toll et al., 2017;Gryspeerdt et al., 2019;Toll et al., 2019). An exception is at large 10 LWP (larger than 200 g m −2 ), where the control simulation is much smaller than the perturbed one, and closer to the satellite retrievals. This is firstly consistent with the expectation that and increase in N d leads to an invigoration of convective clouds and, hence, deeper clouds with higher LWP in the tail of the LWP-distribution, where the convective cloud cores can be found.
It, secondly, reflects the fact that the thick, precipitating clouds most strongly respond to precipitation delay in response to the CCN perturbation.

15
In conclusion, the influence of the perturbed aerosol on N d clearly can be detected and attributed in comparison to satellite retrievals, but the systematic increase in LWP is too small in comparison to natural variability and model bias to be detected and attributed, except for large LWP (> 200 g m −2 ).

Liquid-cloud microphysics in comparison to ground-based remote sensing
Profiles of several cloud microphysical variables (N d , r e , and q l ) were retrieved from ground-based remote sensing as explained in section 2.3.2. In order to derive comparable profiles from ICON-LEM that are best suited for evaluation against the lidar observations, the temporarily high-resolved (9 s) meteogram output was used. All profiles from 21 of the 36 stations for which meteogram output of ICON-LEM was available (Table A1), were in a first step screened for the occurrence of the presence 5 of hydrometeors. The remaining 15 stations were not considered in the analysis because they were too close (approx. within 20 km) to an already considered station, which would lead to an unwanted weighting of the statistics towards a certain region of the ICON-LEM domain. This was specifically the case for the region of the HOPE campaign (Macke et al., 2017), for which output of 13 stations is available. Vertically continuous sequences of hydrometeors were classified as cloud layers to which cloud base height and cloud top height were assigned. Each identified cloud layer was in a subsequent processing step filtered 10 in such a way, that (i) precipitation was absent, (ii) ice was absent, (iii) LWP > 150 g m −2 , (iv) cloud depth < 500 m, and  3.5 Effects on precipitation 10 As discussed earlier, the model simulates large reductions in rain mass and number concentrations, and increases in cloud liquid mass and number concentrations in the perturbed vs. control simulation (Fig. 4). To enquire into the role of the autoconversion process and find out to which extent these changes might be observable, forward simulations of control and perturbed meteogram profiles at the site of the Meteorological Observatory Lindenberg -Richard Aßmann-Observatorium (MOL-RAO) have been performed using the PAMTRA tool . After that, the drizzle detection criterion described in Ac-15 quistapace et al. (2019) has been applied to the forward simulated Doppler radar moments that were generated with PAMTRA.
This approach is similar to the one presented in Rémillard et al. (2017). An ice cloud mask has been applied to the model output in order to filter out the ice pixels and therefore applying the drizzle detection criterion only to the pixels corresponding to liquid in the cloud mask, in both observations and model output.
The distributions for the different classes of drizzle development are shown in Fig. 7 and summarized in Table 5. The postprocessed model results shown in Fig. 7 are consistent with the raw results summarized in Fig. 4: the perturbed simulation (not shown). Mean Doppler velocities are too high (due to too large drops compared to reality) and this produces also broader spectra, giving larger spectrum widths compared to what observed. A closer look into the radar signal suggests that the small reflectivity values for precipitation observations are due to insects detected by the radar. the effort to make model and data comparable, and despite the rather strong signal in the model, this is not yet a useful tool for detection and attribution of differences between control (C2R) and perturbed simulations (P2R) in this case. This tool is at a preliminary stage and can be used to evaluate microphysical schemes to observations (Acquistapace, 2017).

Cloud macrophysics: cloud boundaries, cloud cover and cloud persistence
ICON-LEM cloud base height (CBH), as well as the calculated cloud cover (CC) and cloud persistence (CP) based on the 20 CBH measurements, are compared with high-resolution ceilometer measurements (15 s temporal resolution) from the DWD ceilometer network. In Table 6, mean CBH and CC have been calculated over the 51 stations for 2 May 2013. On average, the model produces less clouds than observed and the cloud base heights are too low in comparison to the ones measured by the ceilometer network. However, the cloud variability is very large, and the model output still is consistent with the data to within the uncertainty range. For both variables, the differences between control (C2R) and perturbed (P2R) simulations are so small, in comparison to the observations, that no significant deviations can be detected given the simulation and observation uncertainties. The absolute (relative) differences between perturbed and control simulations are for the mean CBH -4 m (-0.37 %) and for the CC 0.4 % (0.70 %). The same tendency is found in the all-domain simulation means (Table 7) where the perturbed ICON-LEM shows on average higher cloud tops and bases in comparison to the control simulations. The difference 5 between mean cloud top pressures is -263 Pa (0.35 %) and between mean cloud base pressures, -140 Pa (0.17 %). The perturbed simulation also shows higher total cloud cover by 0.16% (relative difference: 0.20 %) in the domain average (Table 7). Cloud fraction can also be assessed from MODIS satellite data and compared to the COSP-processed ICON-LEM output (as in Section 3.3 for LWP and N d ). The domain-average cloud fraction for the four MODIS overpasses is 0.84, compared to 0.49 and 0.50 for the control and perturbed ICON-LEM simulations, respectively. Despite the very different observational approaches 10 and spatiotemporal sampling, thus, the general conclusion is confirmed: ICON simulates less clouds than observed (a result that also has been noted by Heinze et al., 2017), and shows a positive response of cloud fraction to more aerosol. Possibly the ceilometer network is not able to detect the latter reliably. Further analysis of 30 min periods exhibit that both 15 simulations overestimate "clear sky" cases (0 -5 %) and underestimate "overcast" skies (87 -100 %). "Few" (5 -25 %) and "scattered" (25 -50 %) skies are also overestimated, and "broken clouds" (50 -87 %) are slightly underestimated (not shown). Cloud persistence analysis for the group of few, scattered, and broken cloud conditions (i.e. cloud cover between 5 % to 87 %) and with bases below 3000 m shows also small differences between control and perturbed distributions (Fig. 8). None of the simulations are able to capture the observed distributions for the short-lived clouds (less than 5 min.), while the longer-lived clouds are overestimated.  The comparison of CBH, CC and CP between model and ceilometer measurements shows systematic differences between 5 either model simulation and the observations. Detection and attribution of differences between control (C2R) and perturbed simulations (P2R) is not feasible in this case. This is mostly because the effect of the aerosol perturbation on these quantities is small compared to the model bias.

Sensitivity to cloud regimes
The International Satellite Cloud Climatology Project (ISCCP) regime classification (Rossow and Schiffer, 1991)

Radiative implications
The effective radiative forcing due to aerosol-cloud interactions (ERFaci) has been estimated by subtracting the control simulation (C2R) domain averages from the perturbed (P2R) ones (Table 7). For the simulated case, the ERFaci is -2.62±1.80 W m -2 15 in the TOA net solar radiation (R s toa ), and 0.21±0.40 W m -2 for the TOA net thermal radiation (R t toa ). Consistent with the expectation, the negative forcing in the solar spectrum is slightly reduced by a positive forcing in the terrestrial spectrum (Heyn et al., 2017).
Thanks to the extra simulations that didn't use the number concentration in the radiation transfer computations (C1R and P1R), which also have two different 4D CCN concentration distributions (for 1985 and 2013), the adjustments to the RFaci, 20 as far as they operate via cloud-and precipitation microphysical and dynamical changes, were quantified. In these simulations only the adjustments associated to aerosol forcing are responsible for the radiation changes. The results are noisy signals: the average changes in cloud fraction and LWP are not different from zero to within the temporal variability (not shown). On average, cloud fraction is simulated to decrease slightly (-0.17 %±0.40 %). This result is surprising and different from what is seen in satellite statistics (Gryspeerdt et al., 2019). Further analysis is ongoing. The consequence of the decreasing cloud 25 cover, which is more important radiatively than the increase in LWP, is a positive radiative effect of +0.23 ± 1.24 W m −2 .
As difference between the ERFaci and the adjustments, RFaci, or the cloud albedo effect (Twomey, 1974), is obtained as -2.85 W m -2 .
In order to put this number into context, we assessed the aerosol ERF from four different models from the 5th Coupled -4.0 W m −2 . This latter is larger than the global-annual average present-day vs. pre-industrial, since it is over a region with a large local aerosol perturbation which is even much larger than present-day minus pre-industrial, and the solar zenith angle in May is larger than in the annual mean. In total, a scaling factor of 3.4 is obtained. The ICON-LEM thus implies a global annual mean ERF due to aerosol-cloud interactions in 2000 of about -0.8 W m −2 .

5
This study used a new type of large-eddy simulation which was carried out over a very large domain and driven as a numerical weather prediction with realistic initial and boundary conditions, including an interactive land surface. A large set of observational data from various sources is used aiming for detection and attribution. Four simulations with ICON-LEM model over Table 7. Domain mean differences between perturbed and control simulation pair (P2R -C2R). The absolute changes are given along with the temporal standard deviation of the domain mean changes. The numbers in brackets are the relative changes.