Aerosol indirect effects on summer precipitation in a regional climate model for the Euro-Mediterranean region

Aerosols affect atmospheric dynamics through their direct and semi-direct effects as well as through their effects on cloud microphysics (indirect effects). The present study investigates the indirect effects of aerosols on summer precipitation in the Euro-Mediterranean region, which is located at the crossroads of air masses carrying both natural and anthropogenic aerosols. While it is difficult to disentangle the indirect effects of aerosols from the direct and semi-direct effects in reality, a numerical sensitivity experiment is carried out using the Weather Research and Forecasting (WRF) model, which allows us to isolate indirect effects, all other effects being equal. The Mediterranean hydrological cycle has often been studied using regional climate model (RCM) simulations with parameterized convection, which is the approach we adopt in the present study. For this purpose, the Thompson aerosol-aware microphysics scheme is used in a pair of simulations run at 50 km resolution with extremely high and low aerosol concentrations. An additional pair of simulations has been performed at a convection-permitting resolution (3.3 km) to examine these effects without the use of parameterized convection.

While the reduced radiative flux due to the direct effects of the aerosols is already known to reduce precipitation amounts, there is still no general agreement on the sign and magnitude of the aerosol indirect forcing effect on precipitation, with various processes competing with each other. Although some processes tend to enhance precipitation amounts, some others tend to reduce them. In these simulations, increased aerosol loads lead to weaker precipitation in the parameterized (low-resolution) configuration. The fact that a similar result is obtained for a selected area in the convection-permitting (high-resolution) configuration allows for physical interpretations. By examining the key variables in the model outputs, we propose a causal chain that links the aerosol effects on microphysics to their simulated effect on precipitation, essentially through reduction of the radiative heating of the surface and corresponding reductions of surface temperature, resulting in increased atmospheric stability in the presence of high aerosol loads.

Keywords. Atmospheric composition and structure (aerosols and particles)

1 Introduction

The hydrological cycle in the Mediterranean is a key environmental and socioeconomic question for a wide region including southern Europe, northern Africa, and the Middle East (Drobinski et al., 2014). The Mediterranean climate is characterized by hot summers and mild winters, and by the fact that wintertime is much more rainy than summertime. Historically, it has been proposed to define the Mediterranean climate by the criterion that wintertime precipitation total exceeds 3 times the summertime precipitation total (Köppen and Geiger, 1936). Xoplaki et al. (2004) show that wintertime precipitation accounts for at least 50–60 % of the total annual rainfall in the western and northern parts of the Mediterranean basin and up to 70–90 % in the southern and eastern parts.

There is a need for better understanding precipitation variability, including intensity and frequency, in the Mediterranean region. This is one of the main objectives of the HyMeX (Hydrological cycle in the Mediterranean Experiment) international program (Drobinski et al., 2014). Regional climate models (RCMs) have frequently been used for that purpose (Lionello and Giorgi, 2007; Somot et al., 2008; Flaounas et al., 2013; Berthou et al., 2014; Drobinski et al., 2016, 2017) with horizontal scales ranging from 20 to 50 km and parameterized convection. While direct effects of aerosol have already been examined on the regional scale over the Mediterranean (Nabat et al., 2015, 2016), these simulations do not usually take into account the indirect effects of aerosols. Therefore, the present study aims at exploring the possible impacts of aerosol indirect effects in such models on precipitation over that region. For that purpose, we perform a pair of simulations using the WRF (weather research and forecasting model version 3.7.1; Skamarock et al., 2008) at a 50 km resolution which is within the range of the typical long-term regional climate simulations, with the Kain (2004) convection scheme and the Thompson and Eidhammer (2014) aerosol-aware microphysics. One simulation is forced with extremely strong concentrations of cloud condensation nuclei (CCN), and the other one with extremely low concentrations. The use of this configuration permits us to disconnect the direct effect of the aerosols from their indirect effect, thereby isolating the indirect effect of aerosols on precipitation in a RCM for the Mediterranean area. The processes governing precipitation occurrence markedly differ between wintertime and summertime: while large-scale stratiform precipitation largely dominates wintertime precipitation, the contribution of small-scale convective precipitation is significant in summertime. Therefore, restricting the focus of the study to summertime precipitation permits us to examine in a more detailed way the processes and causal chains by which the indirect effect of the aerosols affect precipitation, allowing to examine both large-scale and parameterized precipitation. As the representation of precipitation in RCMs is very sensitive to the choice of the convection scheme (Di Luca et al., 2014), and also because it has been found that the use of parameterized convection might hinder the correct representation of the microphysical effects of the aerosols (Khain et al., 2015), two additional simulations for a smaller domain at a convection-permitting resolution (3.3 km) have been performed (HR for high resolution). These two simulations allow us to check that the effects observed in the two low resolution (LR) simulations are not mere artifacts of the convection scheme, and also to examine how the aerosol indirect effect affects summertime precipitation in a convection-permitting model.

The Mediterranean area is an area in which it is particularly important to include all aerosol effects in RCM studies because it is located at the crossroads of air masses carrying strong concentrations of both natural and anthropogenic aerosols (Lelieveld et al., 2002; Rea et al., 2015). Natural sources are mineral dust from arid areas, secondary organic aerosols due to biogenic emissions, biomass burning, volcanic emissions, and sea-salt emissions. Anthropogenic sources are mainly due to fossil fuel burning by industrial facilities, power plants, vehicles, and ships, as well as agricultural processes (Querol et al., 2009).

Apart from their direct radiative effect (Schulz et al., 2009; Myhre et al., 2009; Loeb and Su, 2010; Su et al., 2013; Myhre et al., 2013), interactions between aerosols and cloud water and ice occur through several processes. Aerosols act as nuclei for the formation of cloud water and ice, changing the number and size of the cloud droplets and ice particles, which affects the atmospheric physics and dynamics in several ways. Aerosols have been hypothesized to increase the liquid water content, the lifetime of clouds, and cloud height through a decrease in droplet radius and the corresponding decrease of the precipitation efficiency (Charlson et al., 1987; Radke et al., 1989; Pincus and Baker, 1994; Rosenfeld, 2000), an effect often called the Albrecht effect (Albrecht, 1989). However, recent studies suggest that the Albrecht effect on liquid water content could be compensated by enhanced evaporation and dry air entrainment (Small et al., 2009; Seifert et al., 2015). Such an effect may depend on cloud types and properties (Fan et al., 2016) and on the microphysics scheme (Zhou and Penner, 2017). Another process known as the Twomey effect refers to increased cloud optical depth (COD) with increased aerosol concentrations due to the diminished droplet radius and increased droplet number (Twomey, 1977). Under certain conditions, aerosols may also invigorate deep convective clouds by increasing the release of latent energy (Rosenfeld et al., 2008).

Indirect effects are not as well understood as direct effects because many different physical processes are involved and partly compensate each other (Stevens and Feingold, 2009). The buffered characteristic of clouds has been diagnosed in several aerosol indirect effect studies (Fan et al., 2013; Grabowski and Morrison, 2011; Morrison and Grabowski, 2011; Seifert et al., 2012; van den Heever et al., 2011). In numerical simulation, Seifert et al. (2012) observed compensations between several of these feedbacks, which explained a reduced overall aerosol indirect effect on summertime precipitation over Germany. Although similar to our HR simulations, they used day-by-day simulations and on a rather small domain, which potentially prevents the long-term effects described by Fan et al. (2013) from being captured. Our simulation domain and duration are large enough for such behaviors to be possibly observable.

Among them, there is the aerosol cooling effect on the surface caused by reduced solar energy reaching the surface with increased aerosol concentrations. This cooling can be due to the direct effect of the aerosols or to their indirect effects. Surface cooling due to aerosol effects and the consecutive reduction of precipitation have been mostly described on the global scale by studies such as Ramanathan et al. (2001), Salzmann et al. (2014), Bollasina et al. (2011), and by Lelieveld et al. (2002) for the Mediterranean area. On the cloud-resolving scale, the effect on precipitation of the surface cooling due to the indirect effect of aerosols has been studied by Seifert et al. (2012) and Fan et al. (2013). Both of these studies find a surface cooling due to the aerosol indirect effect, with weak effects on precipitation in spite of the stabilization of the atmospheric columns. Morrison and Grabowski (2011) examine the stabilization of the atmosphere due to upper-tropospheric warming as a consequence of the aerosol indirect effects, yielding a small reduction of precipitation. However, stabilization of the atmospheric column does not always cause a reduction of precipitation, since this effect is in competition with the invigoration effect discussed by Fan et al. (2013). This aerosol invigoration effect has been shown to potentially enhance precipitation in certain conditions but cannot act on the global overall amount of precipitation since it does not modify the evaporation flux.

Observational studies of these effects are hindered by the fact that observations can not easily separate the aerosol indirect effects from other effects (Denman et al., 2007); therefore, global or regional models can be used for this purpose. However, quite often, such studies are performed by changing aerosol concentrations simultaneously in the radiative and microphysical schemes of the model, which does not permit us to separate direct and indirect effects from each other. This is the case for studies performed with high-resolution models (Fan et al., 2012; Li et al., 2008; Lerach et al., 2008; Lynn et al., 2005; van den Heever and Cotton, 2007) or low-resolution climate models (Lohmann et al., 2007; Song and Zhang, 2011; Jacobson, 2001). In the present study, as discussed above, we will use the model configurations to modulate the aerosol indirect effects without affecting the aerosol direct effects.

Section 2 details the configuration of the WRF model used and the simulations that have been performed for this sensitivity analysis. Section 3 analyses the sensitivity experiment in both the HR and LR simulations, and compare their sensitivities to the change in CCN concentration. Section 4 discusses the results and proposes a synthesis of the various processes and interrelations involved in the aerosol indirect effects on the Mediterranean summer precipitation.

2 Model presentations and configurations

The model used in this study is version 3.7.1 of the Weather Research and Forecasting model (WRF). WRF is a non-hydrostatic limited-area model with mass coordinate vertical discretization designed to serve both operational and research needs (Skamarock et al., 2008). The WRF simulations presented here have been performed with a 50 km (LR) and a 3.3 km (HR) horizontal resolution and 33 vertical layers from the surface to 50 hPa, with about 8 levels in the lowest 1000 m of the atmosphere. The simulation domain is shown in Fig. 1. The Global Forecast System (GFS) model (National Centers for Environmental Prediction National Weather Service, 2000) provided the initial and lateral conditions updated every 6 h, with a 1^∘ resolution in longitude and latitude. A Newtonian-type nudging with relaxation towards GFS analysis data and a coefficient of $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{5}}$ s⁻¹ has been applied for temperature, humidity, geopotential, and velocity components, as recommended by Salameh et al. (2010) and Omrani et al. (2013, 2015).

The sea surface temperature (SST) is provided by GFS. The geographical data are from 5^′ resolution United States Geophysical Survey data. Soil type is based on a combination of the 10^′ resolution, 17 category United Nations Food and Agriculture Organization soil data and US State Soil Geographic 10^′ soil data. The set of parameterizations used for these simulations include the improved Mellor–Yamada closure scheme (MYNN; Nakanishi, 2000, 2001; Nakanishi and Niino, 2004, 2006) for the boundary layer. The surface layer is the revised Monin–Obukov scheme (Jiménez et al., 2012). The microphysics scheme is the most important for this study. The new Thompson and Eidhammer “aerosol-aware” formulation (Thompson and Eidhammer, 2014) has been chosen for its ability to explicitly represent cloud droplet nucleation and ice activation by aerosols. The scheme gives one-moment prediction (mass mixing ratio) for snow, and a hybrid graupel–hail category and two-moment predictions (mass mixing ratio and number concentration) for cloud water, cloud ice, rain, and aerosols. The latter are divided into two categories depending on their capacity to serve as CCN (“water friendly aerosol”, WFA) or ice nuclei (“ice friendly aerosol”, IFA). During model integration, the number of WFA (NWFA) and the number of IFA (NIFA) are advected and diffused exactly as other scalars. Aerosol number concentration is initialized and forced at domain boundaries by a climatology. For WFA, a virtual surface emission flux is computed from the horizontal grid spacing and the initial NWFA values to approximately balance the loss of WFA due to nucleation and scavenging, as described in Thompson and Eidhammer (2014). No surface emissions are applied for NIFA. The 2-D tendency field is added at each time step to the first model vertical level NWFA value. The radiation scheme is RRTMG (rapid radiative transfer model for general circulation models; Iacono et al., 2008), which uses the correlated-k approach to calculate longwave fluxes and heating rates efficiently and accurately for application to global climate models. We choose this scheme so that the microphysics scheme can communicate cloud water droplets, ice, and snow effective radii to the radiation scheme, which will then use these values to resolve the equations for radiative transfer in each model column, providing quantities such as shortwave (SW) and longwave (LW) radiative fluxes and the COD. It is important to note at this point that the Tegen et al. (1997) climatology provided to the RRTMG radiative scheme will be unchanged when we perform sensitivity experiments by modifying the microphysical NWFA and NIFA climatologies, which enables us to modulate the indirect aerosol effect without changing the direct aerosol effect.

The convection is parameterized using the Kain (2004) scheme. This scheme triggers convection whenever a given atmospheric layer is diagnosed as unstable, taking into account its temperature and moisture, the temperature gap between the adiabatically lifted atmospheric parcel and its environment at its lowest condensation level (LCL), and the large-scale vertical wind speed. If a given atmospheric layer (60 hPa thick) is diagnosed to be able to reach its LCL given the above-mentioned parameters, then it is released at this altitude with an initial vertical velocity of up to several meters per second, depending on the atmospheric conditions. A Lagrangian parcel method including entrainment, detrainment, and water loading is applied. If the atmospheric parcel is able to raise by at least 3–4 km, convection is triggered, and lasts until when the Convective Available Potential Energy (CAPE) is consumed. It is worth noting that at that point the microphysical effects of aerosols are not taken into account explicitly in this parameterization. However, they can modulate the occurrence and intensity of parameterized convection because they affect the background temperature and moisture profiles through the above-mentioned indirect effects, taken into account by the Thompson and Eidhammer (2014) scheme.

Figure 1 Simulation domain with gray-shaded area indicating a topography higher than 500 m. The EUR and MED boxes indicate the regions that will be used below. The Low Resolution domain is the entire map (LR) and the small box inside the EUR box indicate the High Resolution domain (HR).

The model was run with two different NWFA and NIFA microphysical forcings for 6 months covering spring and summer 2013. The two simulations start on 1 April 2013 (after 1 month spin-up) and end on 30 September 2013. We performed two “extreme” simulations in terms of NWFA and NIFA concentrations. In the first simulation, hereafter referred to as MAX or polluted simulation, a very high aerosol emission level is applied (1.75×10⁷ kg s⁻¹ for the whole domain), whereas for the other simulation, hereafter referred to as MIN or pristine simulation, a very low aerosol emission level is applied ($\mathrm{1.75}\times {\mathrm{10}}^{-\mathrm{4}}$ kg s⁻¹ for the whole domain). Although emission rates are extreme between MIN and MAX simulations, the application of the microphysics scheme constrains the range of variation in NWFA to be between ∼10 cm⁻³ and ∼10 000 cm⁻³ and of NIFA to be between 0.005 cm⁻³ and 10 000 cm⁻³, so that the extreme emission rates imposed in the MAX and MIN simulations only ensure that NIFA and NWFA concentrations in the MAX (resp. MIN) simulation stay close to their maximal (resp. minimal) permitted values, corresponding to a factor of ∼10³ for NWFA and $\sim \mathrm{2}\times {\mathrm{10}}^{\mathrm{6}}$ for NIFA between the MAX and MIN simulations. These extreme values ensure that aerosol indirect effects emerge from the “natural noise” between MIN and MAX simulations. Except this difference in the NWFA and NIFA concentrations between the MAX and MIN simulations, the configuration of both simulations is strictly the same, so that all the observed differences between the MAX and MIN simulations can be attributed to these different concentrations of NIFA and NWFA.

Figure 2 Accumulated precipitation (April–September 2013) from the MIN simulation (a) and difference between MAX and MIN simulations of accumulated precipitation (b).

Apart from these two simulations on the LR domain, two companion simulations have been performed at a convection-permitting resolution (3.3 km) on the HR domain, which is a smaller domain (to keep computational cost under control) shown on Fig. 1. The location of the domain was chosen far from the oceans (to avoid sea contamination) and from the edges of the LR domain (to avoid boundary condition contamination). One simulation on this HR domain has been performed with maximal concentrations of NWFA and NIFA, the other one with minimal concentrations. These simulations are forced at their boundaries by the LR MAX and MIN simulations, respectively, through one-way nesting, using an intermediate resolution domain at 16.6 km resolution, also shown in Fig. 1. Due to these successive increases in the resolution by a factor of 5 and then by a factor of 3, each grid cell of the LR domain corresponds to exactly 15×15 grid cells of the HR domain. The configuration of the low-resolution and intermediate resolution simulations is exactly the same. Regarding the high-resolution simulations, the only difference in model configuration is that the Kain (2004) convection scheme is turned off since model horizontal resolution becomes sufficient to explicitly resolve the convection processes.

The analysis of the simulations will be conducted on two distinct sub-domains in Europe (EUR) and the Mediterranean region (MED; see Fig. 1). In the absence of an interactive ocean model, the SST remains prescribed and is the same in both simulations. Therefore, the SST is not in equilibrium with the heat fluxes at the air–sea interface, which is the reason why we will focus all our analyses on atmospheric columns above the continental surface in this study.

3 Sensitivity experiment analysis

3.1 Sensitivity of the LR simulation to the CCN load

Figure 2 displays the geographical pattern of total accumulated precipitation from the MIN simulation. Except for the northeastern part of the Atlantic Ocean, summer precipitation is mostly located over land, particularly over central and eastern Europe as well as over the main mountain ranges. Figure 2 also shows the total accumulated precipitation difference between MAX and MIN simulations (MAX-MIN). Differences are ranging around 5–10 %. Such values are large compared to those found in the study of Thompson and Eidhammer (2014) with the same model and microphysics scheme (the difference between their polluted and clean simulations for a winter cyclone in the USA is 1.1 %). This is as expected because in our study the ratio of aerosol concentrations between our clean and polluted numerical experiments is much larger (10³ vs. 10 in Thompson and Eidhammer, 2014). The pattern of accumulated precipitation differences displays neighboring dipoles of positive and negative differences, suggesting compensating effects or shifts in the precipitation maxima.

Let us first evaluate the first indirect aerosol effect and assess its consequences on precipitation. Figure 3 displays downward surface SW, LW, and total radiation differences as a function of COD differences over the continental regions of the EUR and MED domains. The COD is computed using the following simplified equation:

$$\begin{array}{}\text{(1)}& \mathit{\tau }=\int \mathrm{3}{\displaystyle \frac{\text{LWC}}{\mathrm{2}{\mathit{\rho }}_{\text{w}}{R}_e}}\text{d}z,\end{array}$$

where LWC is the liquid water content, ρ_w the water density, and R_e the droplet radius, and the integral is performed from the surface to the top of atmosphere (or, in our case, from the surface to the model top which is at 50 hPa). Only results for cloud liquid water are shown as results for ice are comparatively negligible. It shows that the effect of this increased COD is twofold, reducing the downward SW radiative flux at the surface (due to increased cloud albedo), but increasing the downward LW radiation due to LW emissions by the optically thicker clouds. However, the effect in the SW band largely dominates the effect in the LW band, leading in our case to a radiative cooling proportional to the change in COD between both simulations. The fact that the modifications of the SW flux dominate those of the LW flux in our study is in line with the findings of Fan et al. (2013) for the United States, China, and the tropical western Pacific. The fact that the modifications of the SW flux dominate those of the LW flux cannot be generalized to all conditions, since the comparison between these two terms depend on several parameters such as the solar zenith angle, surface albedo, cloud optical depth, and cloud altitude, as described by Randall et al. (1984), Stephens and J (1984), Ackerman et al. (1988), and Shupe and Intrieri (2004).

These results are a logical consequence of the Twomey effect, which tends to increase the COD in the presence of stronger concentrations of WFA. On the one side, the increase in the total cross section of cloud droplets leads to higher cloud albedo and therefore to lower downward surface solar radiation. On the other side, the larger the COD, the larger the LW radiation emission by clouds. As discussed before, the difference in surface radiative heating which accounts for both SW and LW radiations remains systematically negative in our study.

Figure 3 Mean downward surface shortwave (ΔR_SW, black), longwave (ΔR_LW, green), and total (ΔR_tot, red) radiation differences as a function of the mean COD difference (Δτ) over the EUR (a) and MED (b) domains.

Let us split total precipitation (ΔP_tot∕P_tot) into their explicit (i.e., large-scale; ΔP_expl∕P_expl) and convective (i.e. parameterized; ΔP_conv∕P_conv) components. Figure 4 shows the relative difference between MAX and MIN simulations of convective, explicit, and total precipitation over continent as a function of the fraction of convective precipitation (P_conv∕P_tot) in the EUR and MED domains (Fig. 1). At each grid point in each domain, the relative differences in accumulated convective, explicit, and total precipitation are paired with convective precipitation fraction. The pairs are placed in bins with an equal number of samples in each bin (50 samples) and therefore varying ranges of convective precipitation fraction for each bin. For each bin, the median convective precipitation fraction and median relative precipitation difference are computed. Summer convective precipitation are clearly weakened by the addition of aerosols, while it is the opposite for large-scale precipitation. Similar results are obtained for low precipitation rate (20th percentile) and high precipitation rate (80th percentile; not shown), which therefore confirm the robustness of the results. Figure 4 suggests a compensating effect between explicit and convective precipitation leading to a non-significant difference for total accumulated precipitation. To explain this, a similar pair of simulations have been performed without activating the convection scheme. The analysis of these two simulations also confirm an increase in large-scale precipitation when aerosol concentrations increase, thus discarding the compensating effect (not shown). On the contrary, it suggests an independent positive effect of aerosols on non-convective precipitation.

Figure 4 Relative difference between the MAX and MIN simulations of continental summer convective (ΔP_conv∕P_conv, dotted), not convective (ΔP_expl∕P_expl, dashed), and total (ΔP_tot∕P_tot, solid) accumulated precipitation as a function of convective precipitation fraction (P_conv∕P_tot) over Europe (EUR domain, Fig. 1) (a) and the Mediterranean region (MED domain, Fig. 1) (b).

Since 50–85 % of the radiative heating is balanced by evaporation (Kiehl and Trenberth, 1997), a correlated modification of low-level water vapor mixing ratio should be observed. Figure 5 displays the relationship between the difference between MAX and MIN simulations of convective precipitation (ΔP_conv) and of surface water mixing ratio (Δq). For each grid points, the differences in hourly convective precipitation (ΔP_conv) are paired with the difference of surface water mixing ratio (Δq). The pairs are placed in bins with an equal number of samples in each bin (20 000 samples). For each bin, the median values are computed as well as the 20th and 80th percentiles and are displayed in Fig. 5. Figure 5 shows that a decrease (increase) in surface water mixing ratio is associated almost systematically with a decrease (increase) in convective precipitation. Despite less local evaporation, some convective events are characterized by higher available surface water vapor with an associated increase in convective precipitation. The origin may be a “contamination” by the sea. Indeed, since the SST is imposed and is not in equilibrium with the reduction of downward solar radiation and because evaporation flux at the surface of the sea is in part driven by the temperature difference between the air and the sea surface (e.g., Mahrt and Ek, 1984; on which the Unified NOAH Land Surface Model – which is used in this study – is based), erroneous evaporation fluxes are expected to be simulated over the sea.

In the coastal land of the Mediterranean region in summertime, sea breezes advect marine moist air over a few tens of kilometers inland (e.g., Bastin et al., 2006, 2007; Drobinski et al., 2006, 2017), which can cause inland precipitation. The absence of air–sea feedbacks in this simulation therefore affects inland humidity advection and convective precipitation in the coastal areas (e.g., Lelieveld et al., 2002; Lebeaupin Brossier et al., 2013).

In these simulations, aerosols reduce total precipitation through their parameterized part. In order to check if this effect is not only a parameterization effect, we performed a 3.3 km resolution simulation without activating the convective scheme. The following part aims at determining if the same precipitation reduction effect is also present in the convection-permitting simulation.

Figure 5 Difference between MAX and MIN simulations of convective precipitation (ΔP_conv) and of surface water mixing ratio (Δq) for EUR (a) and MED (b) domains. The thick black line corresponds to the ensemble average. The lower and upper dotted lines delimiting the shaded area are the 20th and 80th percentiles, respectively.

Table 1 Comparison of precipitation statistics between LR and HR simulations. Space-averaged temporal correlation for the hourly, daily, and weekly precipitation amounts, respectively (first row); same for the time-averaged spatial correlation (second row); accumulated precipitation (mm) averaged over the HR domain from the whole LR and HR simulation time series (third and fourth rows).

3.2 Sensitivity of the HR simulations and impact of the resolution

At this point, it is suitable to check how the LR simulation behaves in terms of precipitation amount, localization, and timing compared to the convection-permitting HR simulation. Figure 6 shows daily MIN time series of precipitation for a LR grid point that is inside the HR domain. HR precipitation time series were made considering the average of the 15×15 grid cell that are inside the LR grid cell. This time series shows that the LR simulation is able to reproduce the main features of the HR simulation, including the succession of long-lasting rainy and dry sequences, as well as the occurrence of precipitation maxima for individual days, even though some of these peaks are stronger in the LR simulation than in the HR simulation or conversely. Table 1 shows a comparison of temporal and spatial statistics of precipitation between the LR and HR simulations. In the first line we can see the spatial average of the hourly, daily, and weekly correlations of precipitation amounts between LR and HR simulations, performed by first calculating the correlation of the hourly (resp. daily, weekly) precipitation amount between each of the cells of the LR model vs. the corresponding precipitation amount for the 15×15 matching cells of the HR domain, and then averaging all these local time correlations over the whole intersection of the HR and LR domains.

One can note a correlation around 0.3 in both MIN and MAX simulations on the hourly timescale between the LR and the HR simulations. The correlation coefficient increases for longer timescales, reaching 0.57 on daily timescale and around 0.71 on weekly timescale for both MAX and MIN simulations, which shows that even though the LR simulation does not succeed in having the precipitation events with exactly the same timing and intensity as in the HR simulation, these differences tend to be strongly reduced on the weekly timescale, which is satisfying since RCMs are usually used for climate studies.

Figure 6 Daily MIN precipitation time series for the LR (blue) and HR (red) simulations in one LR grid cell located close to the center of the HR domain.

Let us now evaluate the spatial distribution of precipitation in the LR simulation. Figure 7 displays the MIN hourly cumulated precipitation for the day of 12 May 2013 in the HR domain, for the HR simulation (top left), the LR simulation (top right), and for the HR precipitation aggregated on the LR grid by averaging on each LR grid cell the 15×15 matching HR grid cells (bottom right). We can see that the spatial distributions of the HR and LR simulations are rather similar on a large scale, even though the precipitation tends to be spread over larger areas in the LR simulation than in the HR simulation, even after aggregation of the latter on the LR grid, with a corresponding reduction of the maximal values of precipitation rates in the LR simulation.

As seen in the second line of Table 1, averaged spatial correlation increases from around 0.30 on an hourly timescale to around 0.65 on a weekly timescale for both MIN and MAX simulations. Similarly, averaged spatial correlations increase when we cluster grid points (not shown). It shows that even if precipitation is not exactly located at the same grid point and at the right time, spatial patterns are globally respected as seen in our example (Fig. 7). Therefore, as for the temporal correlations, the examination of spatial correlations shows that the spatial agreement between the HR and LR becomes satisfactory when the analysis is focused on daily or weekly timescales, which once again is in line with the typical uses of long-term RCM simulations.

Figure 7 Precipitation amount (in mm) between 13:00 and 14:00 UTC on 12 May 2013 for the MIN HR simulation (a), MIN LR simulation (b) and the HR simulation aggregated on the LR grid by averaging on each LR grid cell the 15×15 matching HR grid cells (c).

A comparison of the accumulated precipitation over the whole HR domain for the entire 6 months of simulation between the HR and LR simulations (Table 1) reveals that the accumulated precipitation in the LR simulation is significantly stronger than in the HR simulation, both for the MIN and MAX experiments. However, the indirect effect of aerosols on the precipitation is similar in the LR simulation than in the HR simulation (−11 % in the HR simulation and −7 % in the LR simulation). This shows that even though the LR simulation with the Kain–Fritsch convection scheme simulates stronger precipitation than the convection-permitting HR simulation, the microphysical effects of aerosols on precipitation in the LR simulation are of the same sign and magnitude than in the HR simulation, which tends to strengthen the results presented for the LR simulation in Sect. 3.1.

4 Discussion

Table 2 displays various key variables averaged over all of the simulation period and over the EUR and MED domains. By examining the sensitivity of these key variables to the CCN load in both the HR and LR simulations, we will try to discuss the plausible processes involved in the simulated reduction of precipitation due to the aerosol indirect effects, which we sum up in Fig. 8 – a summary of most of the effects discussed in the present study.

Table 2 Key variables for the MIN and MAX simulations: number of cloud condensation nuclei (NCCN), number of liquid cloud drops (N_drops), liquid water content (LWC), cloud drops effective radii (R_eff), cloud optical depth (COD), shortwave downward radiations (SWD), longwave downward radiations (LWD), latent heating (LH), surface temperature (T_surf), surface water vapor mixing ratio (Q_surf), number of hours per grid point with strictly positive MuCAPE (N_hour MuCAPE), MuCAPE, MuCIN, convective precipitation (P_conv) and total precipitation (P_tot). All the values are averaged in space and time over the indicated domain from April to September 2013, except P_tot and P_conv (precipitation totals averaged over space only) and N_hour MuCAPE (number of hours averaged over space only).

The ratio of CCN between the MIN and MAX simulations is about 1000, with no major difference between the EUR and MED domains as the forcing climatology is uniform over the whole domain. In detail, the NCCN is larger in the MED region in spite of the identical forcing emissions since a smaller NCCN is used for condensation, and wet scavenging of CCN is less active due to the weak amount of precipitation (Fig. 2a). As expected, the number of liquid cloud droplets (Ndrops) increases massively with increasing aerosol concentration (this number is multiplied by 40 to 50 between the MIN and MAX experiment in the three examined simulation pairs). One can note that the relative difference in Ndrops is smaller than that of NCCN as only a fraction of aerosols is used for condensation, especially in the MAX simulation in which the number of aerosols does not limit the condensation of the available water. An increase in the number of droplets due to increased availability of CCN is the first step in the process relating the increase in aerosol concentration to the reduction of convective precipitation as shown in Fig. 8. This larger number of droplets is associated with smaller droplet radius (Reff) as a given amount of water is distributed over more droplets (Table 2 and Fig. 8). The relative difference between the MAX and MIN simulations of the droplet radii is significant and on average equal to −60 %.

Table 2 also shows that the available amount of water is sensitive to the aerosol load, with an increase in the liquid water content (LWC) in the MAX experiment by about 78 to 190 % in the three simulation pairs. An increase in the cloud lifetime in the MAX simulation might be the main explanation (Warner, 1968; Albrecht, 1989; Rosenfeld, 1999). However, on the one hand, the aerosol cloud lifetime effect has been challenged in recent studies (Small et al., 2009; Seifert et al., 2015), and on the other hand, even if this effect stands in reality, it is not sure that our model is able to correctly simulate it since the representation of entrainment and droplet evaporation has been showen to be poor in parameterized models (Zhou and Penner, 2017).

Equation (1) shows that the COD is proportional to the LWC and inversely proportional to the droplet radius. Therefore, the increase in the COD between the MIN and MAX experiments (COD is multiplied by 4.9 up to 6.4 in the three simulation pairs) is in part due to the increase in the LWC, and in part due to the diminished droplet radius, these two effects being of similar magnitude in our simulations. As a consequence of a massive increase in the COD, the surface downward shortwave radiation (SWD) is reduced by about −1 to −6 % (Table 2 and Fig. 8).

Figure 8 Schematic summary of the proposed aerosol causal chain sequence for the indirect effect of aerosols on convective precipitation

The decrease in SWD implies a decrease in the surface temperature (Table 2) with the same order of magnitude (−0.2 to −0.6 K) as in Seifert et al. (2012) in Germany. The combination of SWD decrease and surface temperature decrease sum up to cause a decrease in surface evaporation (LH) by about 3 to 6 % due to the smaller quantity of energy available for evaporation and to the reduced capacity of the colder air to transport water vapor (Clausius–Clapeyron law). A consequence is the reduction of the surface water vapor mixing ratio. On average, it is the case over the EUR domain (−0.8 %) but not over the MED domain (+0.2 %) due to the “sea contamination” (prescribed sea surface temperature; Table 2). A reduced surface water vapor mixing ratio is associated with reduced convective precipitation (Fig. 5; with the exception of the contaminated coastal convective precipitation). The reduction of precipitable water (increase in LWC does not compensate the decrease in water vapor) available for convection partly inhibits convective precipitation (Fig. 8). The terms Most unstable Convective Available Potential Energy (MuCAPE) and Most unstable Convective Inhibition (MuCIN) have been chosen to characterize the stability of the atmosphere. MuCAPE represents the total amount of potential energy available to the most unstable parcel of the atmospheric column while being lifted to its level of free convection, and MuCIN represents the energy barrier needed for this parcel to reach its level of free convection. By construction, the MuCAPE parameter is always positive. Therefore, the difference in MuCAPE between both simulations can only be non-zero when there is a potentially unstable layer in at least one of the simulations. MuCAPE differences between both simulations at a given time and place can reveal very diverse situations. If, on one hand, MuCAPE_MIN>0 and MuCAPE_MAX=0, there is a potentially unstable atmospheric layer in the MIN simulation but not in the MAX simulation (and conversely). If, on the other hand, MuCAPE_MIN>0 and MuCAPE_MAX>0, then the atmospheric column is potentially unstable in both simulations, and MuCAPE_MIN−MuCAPE_MAX describes the difference in the potential intensity of convection between both simulations.

To only analyze the potential intensity of the convection, we averaged values of MuCAPE_MIN−MuCAPE_MAX when both energies are non-zero (which makes a sample of about 2.4 million hours spread over 1499 grid points in the EUR domain and about 1.1 million hours spread over 1606 grid points in the MED domain).

On the other hand, MuCIN is defined only when a level of free convection exists, i.e., only when MuCAPE is strictly positive. For all these reasons, we averaged MuCAPE and MuCIN considering only the events with strictly positive MuCAPE in both simulations. Because of increased relative humidity, MuCIN is reduced by about 4 to 6 % in the LR simulations, and up to 16 % in the HR simulation. MuCAPE is also weaker by 6 to 8 % in the LR domain, and up to 19 % for the HR simulation, mostly due to the lack of water vapor and to lower surface temperatures. When convection is actually triggered, this reduced MuCAPE induces weaker convective updraft and therefore weaker convective precipitation. This reduction in convective precipitation (Fig. 8) is of 16 % in the EUR domain, and 12 % in the MED domain, and causes a reduction in total precipitation of 5 to 8 % in the EUR and MED domains of the LR simulation, and 11 % in the HR simulation (Table 2).

The effect of the stabilization of the atmospheric columns due to the indirect effect of the aerosols has already been examined by Morrison and Grabowski (2011), Seifert et al. (2012), and Fan et al. (2013). While Morrison and Grabowski (2011) find small reductions in the precipitation amounts due to this stabilization, Seifert et al. (2012) and Fan et al. (2013) find a slight increase in precipitation in a polluted atmosphere, due to the convective invigoration linked to increased ice formation in the upper troposphere, an effect that cannot be represented in our LR simulation since the Kain–Fritsch convective scheme we use here does not explicitly take into account the effect of aerosols and nucleation. While it can be represented by our HR simulation, it does not appear to play a dominant role in the conditions of our study.

5 Conclusion

Aerosols affect atmospheric dynamics and precipitation through their direct and semi-direct radiative effects as well as through their indirect effects. In this study, we investigated the indirect effects of aerosols on summer precipitation in the Euro-Mediterranean region. While it is difficult to disentangle the indirect effects of aerosols from the direct and semi-direct effects in reality, we have carried out a numerical sensitivity experiment using the WRF model, which allows us to isolate indirect effects, all other things being equal, while introducing uncertainties linked to the model formulation. Very large and unrealistic perturbations were chosen in order to allow the indirect effect of aerosols to emerge from these uncertainties. For this purpose, the Thompson aerosol-aware microphysics scheme has been used in a pair of simulations run at 50 km resolution for the spring and summer of 2013 (April–September) with extremely high and low aerosol concentrations, respectively. The outputs of these simulations have been analyzed in two areas, one over continental northern Europe and the other for the continental areas of the Mediterranean basin (Fig. 1). In order to strengthen the results obtained on precipitation with this first set of simulations, we performed an additional pair of simulations at a convection-permitting resolution (3.3 km) for a continental subdomain.

Analysis of the LR simulation outputs has revealed two opposite responses to this increase in aerosols concentrations. In the polluted simulation, the amount of convective (parameterized) precipitation is reduced in comparison to the pristine simulation (Fig. 4). The total precipitation is also reduced, even though an increase in stratiform precipitation partly compensates the reduction of convective precipitation (Table 2). An examination of the differences in several relevant variables between the pristine and polluted simulation reveals that the polluted simulation has more water droplets than the pristine simulation, and that these water droplets tend to be smaller (Table 2). Additionally, the liquid water content of the atmosphere is also generally increased in the polluted simulation (Table 2), probably in part because the smaller droplets tend to stay longer in the atmosphere before precipitating. Both the quantity and the size effect contribute to an increase in the COD, and a subsequent reduction in the SWD flux at the surface, which is only partly compensated by a much smaller increase in the net longwave radiative flux (Fig. 3 and Table 2). The reduction of radiative heating at the surface has consequences on the atmospheric stability through a reduced evaporation and lower surface temperature in the polluted simulation (Table 2). The reduced evaporation leads to a lower surface water vapor mixing ratio in the polluted simulation (Fig. 5 and Table 2) that accounts not only for the reduction of potential energy for convection (Table 2) but also the reduction of precipitable water. The increase in stability and the reduction of surface water vapor content contributes to the reduction of convective precipitation in the polluted simulation (Table 2). A visual summary of this causal chain has been presented in Fig. 8. In these simulations with parameterized convection, the stratiform (explicitly resolved) precipitation increase partly compensates the decrease in convective (parameterized) precipitation. Since it could be envisioned that this increase in stratiform precipitation could just be a consequence of the reduction of convective precipitation, for example through a greater availability of water for stratiform precipitation due to reduced convective rainfall, we performed an additional pair of pristine and polluted simulations at the same resolution but without activating the convection scheme. Analysis of these simulations has shown that an increase in stratiform precipitation also occurs in this additional pair of simulations where parameterized convection is turned off, thereby excluding the hypothesis that the increase in stratiform precipitation is just an effect of decreased convective precipitation. Therefore, the reasons for this increase in large-scale precipitation in the polluted case remains an open question.

In order to confirm and strengthen this result of a net reduction of total precipitation due to aerosol secondary effect in a regional climate model, at least in the study area, and to show that this result is not just a consequence of the choice of a particular scheme for parameterized convection (in our case the Kain-Fritsch scheme), we analyzed the same variables as in the low-resolution simulation in the pair of convection-permitting simulations. Even though differences exist between the behavior of the convection-permitting simulation pair compared to its large-scale counterpart, the main features are similar to those of the large-scale simulations: a strong increase in the number of cloud droplets and their liquid water content in the polluted simulation compared to the pristine simulation, a strong decrease in their size, a decrease in the SWD radiative flux, a decrease of the surface temperature and surface humidity, and a marked decrease in the MuCAPE, yielding an overall reduction of total precipitation during the simulation period (Table 2). Convective invigoration, which numerically can not occur in LR simulations because the cumulus scheme is insensitive to the aerosol concentration, actually may not be important in the convection-permitting simulation either. This result was expected since it has been shown that this invigoration effect is stronger in weak wind shear and warm cloud-bases (> 15 ^∘C) conditions (Li et al., 2011; Fan et al., 2009), which may not be the situations that occur most often in mid-latitude regions. Rosenfeld et al. (2008) evaluated the amplitude of the invigoration effect as a function of the CCN concentration. They claimed that the weight of the parcel condensate increases with CCN concentration since aerosols tend to suppress precipitation by the Albrecht effect. The increasing weight of parcels with CCN load tends to counteract the release of more freezing latent heating (due to invigoration effect) with increasing CCN. They stated the existence of an optimal CCN concentration at 1200 cm⁻³ for the strongest invigoration effect. Beyond this optimal CCN concentration value, the invigoration effect tends to be weakened by the increasing weight of parcels. The effect is even negative when reaching extreme values such as 10 000 cm⁻³, which is the typical order of magnitude of our MAX simulation (see Table 2). This is another possible explanation for the weak invigoration effect observed in our simulations.

Even though the reduction of precipitation occurs through their parameterized part, the fact that the reduction of the total precipitation amount persists in a high-resolution convection-permitting simulation, at least in the central Europe domain used, suggests that this impact of indirect aerosol effects on precipitation is not just a consequence of the formulation of the convection scheme but reflects underlying physical processes for which we propose a summary under the form of a causal chain (Fig. 8).

The magnitude of the simulated impact of the secondary aerosol effects on precipitation is estimated as a reduction of total precipitation by 5 to 11 % in our simulations, but since these simulations have been performed with extremely low or extremely high concentrations of CCN, these values need to be interpreted as an upper bound of the magnitude of this effect rather than as a realistic estimate. Our result can also be seen as a buffered response of atmosphere to the aerosol indirect forcing as mentioned in Stevens and Feingold (2009). Indeed, analysis of the LR simulation shows that stratiform clouds act in two opposite ways for surface precipitation in a polluted environment: in a direct way they produce more stratiform rain, and in an indirect way they decrease convective rain through their stabilization effect. This compensation explains why extreme changes in aerosol concentrations produce only small changes in total precipitation. If all of the mentioned aerosol indirect effects have already been discussed in past studies, in particular those of Seifert et al. (2012) and Fan et al. (2013), the present study makes it possible to analyze which effects are dominant for summertime in the Euro-Mediterranean area with a regional simulation, and comparing a low-resolution simulation with parameterized convection to a cloud-resolving simulation. Compared to studies such as the ones of Morrison and Grabowski (2011), Seifert et al. (2012), and Fan et al. (2013), our study permits us to observe this effect in the same region and period with these two different model configurations, showing that a model configuration with parameterized convection can represent precipitation reduction due to stabilization related to the indirect effect of aerosols in a way that is very similar to a cloud-resolving simulation. Therefore, the present study suggests that, for the mid-latitudes in summertime, global or regional climate models with parameterized convection may be able to capture the precipitation changes due to the effect of radiative cooling by the aerosol indirect effect in a realistic way, at least when compared to similar models in a convection-permitting configuration over the same area. Using such configurations with parameterized convections also permits us to clearly separate the processes affecting stratiform precipitation from the ones affecting convective precipitation, as we have done above.

Because of the time of computation, our high-resolution domain is rather small, which would limit the physical interpretation to that selected area. However, the typical scale of the precipitation difference patterns in Fig. 2 is less than half the size of the HR domain, so that the differences observed in the HR domain do not reflect only a particular case of positive or negative effect on precipitation but a reasonable sample of continental Europe. Of course, for future studies, it would be interesting to have such a high-resolution study for an entire continental-scale domain.

To strengthen the analysis of these processes and evaluate in a more realistic way their magnitude, it will be necessary to perform simulations with models that are able to take into account realistic aerosol levels as well as their variations in space and time, which will be possible in particular through the use of online-coupled models that treat aerosol emission and transport in a realistic way, such as described in Baklanov et al. (2014). Also, our study was not able to examine these effects over oceanic surfaces since the forcing by fixed SSTs that are not in equilibrium with the perturbed atmospheric forcings does not permit to obtain meaningful results over the oceans; this could be overcome by using models which are coupled with an oceanic model. Finally, to resolve the full lifetime of clouds and water vapor in the atmosphere, it would also be useful to use a larger domain (global or hemispheric).