Feedback of observed interannual vegetation change: a regional climate model analysis for the West African monsoon

Two WRF ensemble experiments, driven by ERA-Interim reanalysis data (ERA-I) (Dee et al. 2011), are conducted with (1) the observed state of surface variables (ALB, VF, LAI), that was derived from remote-sensing (DYN) and (2) default climatological datasets with a fixed annual cycle (CLIM) for the rainy seasons April–September 2009 and 2010, including four months of spin-up. The model analysis focuses on August–September, when vegetation anomalies reach their maximum (cf. Fig. 3). The two rainy seasons were simulated separately instead of conducting a continuous simulation in order to have the same integration time after the first soil moisture initialization for both years and to receive eight independent simulations per ensemble.

\(\Delta {\mathrm {Y}}_{\mathrm {Srfc}}\) represents a “vegetation-induced modulation of the large-scale driven changes” between 2009 and 2010. For example, for a temperature change from 25 to \(30\,^\circ \)C in CLIM and a change from 25 to \(26\,^\circ \)C in DYN, we obtain \(-4\,^\circ \)C for the vegetation signal \(\Delta {\mathrm {Y}}_{\mathrm {Srfc}}\), which means that \(\Delta {\mathrm {Y}}_{\mathrm {Srfc}}\) can be negative although DYN exhibits an effective net warming.

The WRF model domain encompasses the whole WAM region as can be seen in Fig. 1. If not indicated otherwise, all analyses are carried out for the domain (\(9^{\circ }\)W–\(9^{\circ }\)E, 7–16\(^{\circ }\)N), focussing on the Sudanian zone (7–12\(^{\circ }\)N) and the Sahel (12–16\(^{\circ }\)N) to reduce effects from the ocean and the coastline.

The provided boundary conditions from ERA-I include sea surface temperatures and lateral atmospheric variables that are updated every 6 h. In addition, soil moisture and soil temperatures are initialized from ERA-I for all simulations. The WRF model output is written at an hourly basis.

Different model configurations and their impact on the representation of the WAM dynamics were already investigated by Klein et al. (2015). On this basis, we run the WRF model with the RRTM/Dudhia Long and Shortwave Radiation schemes (Dudhia 1989), the LIN Microphysics scheme (Lin et al. 1983) and the YSU planetary boundary layer scheme (Hong et al. 2006). The YSU scheme is a first-order closure scheme that also takes into account non-local vertical mixing in convective boundary layers by including a counter-gradient term in the vertical diffusion equation. The surface exchange coefficients are provided by the MM5 surface layer scheme based on the Monin–Obukhov similarity theory (Zhang and Anthes 1982), which is tied to the YSU scheme.

This parameterisation configuration reproduces the WAM without showing an extreme wet or dry behaviour. The horizontal resolution is 7 km without using a convective parameterisation to allow the model to explicitly resolve convection. An explicit treatment of convection is especially important for land–atmosphere interaction studies: the diurnal cycle of cloudiness and rainfall affecting the available energy at the surface (Marsham et al. 2013) and planetary boundary layer processes affecting atmospheric stability (Hohenegger et al. 2009) can realistically be represented with explicit convection but are not properly captured by current cumulus parameterisations.

The atmospheric part of the WRF model forces the Noah LSM with atmospheric variables such as incoming long- and short-wave radiation, near surface temperature, pressure, humidity and precipitation. The Noah LSM then calculates the outgoing radiation with respect to ALB and TS of the previous time step and thus provides the radiative forcing \(R_{net}\) at the ground (Eq. 2). Altogether, four variables are fed back to the atmosphere: (1) LH (evapotranspiration), (2) SH, (3) outgoing long-wave radiation (via calculation of TS) and (4) upward short-wave radiation (reflected portion given by surface albedo).

This summary of variable dependencies describes how state variables and surface fluxes directly react on changes in ALB, VF or LAI: ALB has a direct effect on \(\mathrm {R}_{\mathrm {net}}\) and therefore impacts the available energy for the turbulent surface fluxes. VF is the key parameter for the partitioning between bare soil evaporation and plant transpiration. The efficiency of the latter, and therefore the actual difference between bare soil evaporation and plant transpiration, is modified by the canopy resistance and thus by LAI. However, in this study the objective is not to disentangle the separate effects of these intrinsically connected surface parameters, but to investigate whether their “bulk effect” leaves an imprint on the atmosphere. As an approximation, we therefore use VF changes as a proxy for associated changes of all three variables in our feedback analyses.

The model simulations are compared to different satellite, observational and reanalysis data to evaluate their capability to capture the monsoon regime and interannual differences between 2009 and 2010.

For precipitation, the NASA Tropical Rainfall Measuring Mission (TRMM) 0.25\(^{\circ }\) resolution 3B42V7 rainfall estimates (Huffman et al. 1995, 1997), the African Rainfall Estimate version 2 (RFE) (NOAA CPC 2013) \(0.1^{\circ }\) combined satellite-gauge data product and the \(0.25^{\circ }\) resolution daily global CMORPH (CPC Morphing technique, Joyce et al. 2004) are used. All three datasets are available at daily resolution.

Model surface temperatures are compared to the Global Historical Climatology Network (GHCN) \(5^{\circ }\) gridded temperature product Version 3 (Lawrimore et al. 2011), the University of Delaware 0.5\(^{\circ }\) station-analysis product (Legates and Willmott 1990) and the station-data based Climatic Research Unit (CRU) TS3.21 \(0.5^{\circ }\) product (Harris et al. 2013). All three datasets come at monthly resolution.

Due to the small amount of available flux tower data in West Africa, the information on surface fluxes is probably the most uncertain. The only available gridded product for latent and sensible heat fluxes is the FLUXNET Multi-Tree Ensemble (MTE) at \(0.5^{\circ }\) (Jung et al. 2011) on a monthly basis. Additionally, we use two reanalysis-based products ERA Interim/LAND (Balsamo et al. 2015; Albergel et al. 2013) and MERRA-Land (Reichle et al. 2011) at monthly resolution. Both products are land-only reanalysis products from off-line simulations with updated land surface models and a new merged precipitation forcing from gauge data and reanalyses. Monthly soil moisture estimates are taken from the \(0.5^\circ \) gridded CPC soil moisture (V2) dataset (Van Den Dool et al. 2003). The monsoon dynamics are compared to ERA-I data. If not indicated otherwise, all data, including WRF, is transformed (averaged) to temporally and spatially match the coarsest available resolution of the considered datasets for each individual analysis.

The \(\mathrm {LAI}_{\mathrm {Dyn}}\) and \(\mathrm {VF}_{\mathrm {Dyn}}\) time series were created by fusing LAI and VF datasets derived from data of the Spinning Enhanced Visible and Infrared Imager (SEVIRI; LSA SAF 2013) and Satellite Pour l’Observation de la Terre-VEGETATION (SPOT_VGT; Baret and Weiss 2010, 2014; Camacho and Cernicharo 2015). The SPOT_VGT-based datasets provide good quality information on LAI and VF in 10-day intervals at 1 km spatial resolution, but show frequent cloud gaps during the rainy season, while the SEVIRI-based daily datasets exhibit fewer cloud gaps but at a coarser spatial resolution of 3 km (Camacho and Cernicharo 2015; Gessner et al. 2013). Both datasets were aggregated to a spatio-temporal resolution of 1 km and 10-daily intervals. The SEVIRI datasets were matched to SPOT_VGT values based on slope and offset of three linear regressions fitted for gap unaffected areas and three seasons (Jan–Mar, Apr–Sep and Oct–Dec). Finally, the cloud gaps in the SPOT_VGT products were filled with the adapted SEVIRI-based LAI and VF values. Remaining outliers in the time series were removed and remaining gaps were filled by temporal linear interpolation.

The \(\mathrm {ALB}_{\mathrm {Dyn}}\) time series is based on the bi-hemispherical reflectance (white-sky albedo) for the shortwave spectral range of the MODIS MCD43B3 product (LP DAAC 2014). The MCD43B3 is a composite product based on 16-day intervals at 1km spatial resolution. Like for LAI and VF, the product features gaps, especially during the wet season. To generate a gap free input, we aggregated the products to monthly composites by calculating the mean of corresponding months. The remaining data gaps were filled with long-term monthly mean values (15 years: 2000–2014). Finally, the time series underwent temporal smoothing using the Savitzky-Golay filter (William et al. 1992).

The DYN datasets for ALB, VF and LAI only span the focus study region \(9^{\circ }\)W–\(9^{\circ }\)E, 4–16\(^{\circ }\)N. In the rest of the WRF domain, global datasets (GLOB, cf. Table 1) are used to provide dynamical surface information. Remaining gaps due to cloud cover in GLOB are filled with the GLOB long-term mean value for the respective pixel and month.

All datasets shown in Table 1 are linearly interpolated to 6-h timesteps. They act as boundary condition for the WRF model and are updated every 6 h. Spatially, the datasets are interpolated (coarser datasets) or regridded and averaged (finer datasets) to the model grid at 7 km resolution. While the averaging implies a loss of information for LAI, it preserves a proper grid-scale representation of the variables affected by VF and ALB (Kunstmann 2008).

Other necessary land surface parameters (e.g. surface roughness, root zone depth) rely on table values that are associated with a certain land use classification. Therefore, a regionally optimized land use classification (DLC) map at 250 m spatial resolution, based on the year 2006 (Gessner et al. 2015), is additionally implemented in the model in order to ensure consistency of the assigned table values with the new DYN datasets. Figure 1b shows the 7 km upscaled maps of DLC and the WRF default 1 km MODIS land classification map based on 2002 (MLC, Friedl et al. 2010). The 7 km upscaled DLC and MLC consider the dominating land class only.

There is a considerable albedo offset between DYN and CLIM, illustrated in Fig. 2, with a mean deviation of 2.4 % in the Sahel and 1.9 % in the Sudanian zone, which arises because \(\mathrm {ALB}_{\mathrm {Dyn}}\) is based on MODIS shortwave broadband white-sky albedo, while \(\mathrm {ALB}_{\mathrm {Clim}}\) represents AVHRR clear-sky albedo. This offset and its physical effects hinders a direct comparison of the DYN and CLIM simulations. Instead, we compare the interannual change of the respective simulation case (cf. Sect. 2.1), assuming that for \(\Delta \mathrm {Y}\), the effect of this offset is negligible.

Differences between \(\mathrm {VF}_{\mathrm {Dyn}}\), \(\mathrm {LAI}_{\mathrm {Dyn}}\) and the climatological dataset are smaller, but show pronounced higher maxima during the rainy season. As expected, the standard deviations are generally higher for the DYN datesets compared to the climatological annual cycles.

The DYN monthly anomalies for the Sahel in Fig. 3a show good agreement with each other as well as with observed precipitation (PCP) and soil moisture estimates (SM), giving confidence in their quality. \({\mathrm {LAI}}_{\mathrm {Dyn}}\) co-varies with the \({\mathrm {VF}}_{\mathrm {Dyn}}\) anomalies and is therefore not shown. In the Sahel, SM and \({\mathrm {VF}}_{\mathrm {Dyn}}\) anomalies show maximum correlations with PCP anomalies of 0.71 and 0.74 with a lag of one month, as shown in Fig. 4. Even after two months, the correlations do not drop below 0.55. In contrast, Lohou et al. (2014) identified an immediate increase of the evaporative fraction at flux station sites in Niger and Mali with a recovery time to its original state of only 1–3 days. This illustrates that in the Sahel, changes in vegetation cover after rainfall events are much slower but longer lasting than changes in surface soil moisture.

In the generally moister and more densely vegetated Sudanian zone, VF reacts with a smaller temporal lag to PCP anomalies shown by similar correlations for a lag of 0 and 1 month. However, corresponding to Camberlin et al. (2007), vegetation changes are generally less sensitive to PCP anomalies (r ~ 0.47) since water availability is a weaker constraint in this region than in the Sahel.

Nevertheless, this illustrates that vegetation may retain a long-term memory of preceding rainfall anomalies over the whole region. We hypothesize that surface patterns that were induced by precipitation in the early months of the rainy season may affect the patterns of rainfall in the late rainy season. We therefore focus on the period from August to September (Aug–Sep), when vegetation anomalies reach their maximum (cf. Fig. 3) and the error that is introduced by climatological datasets should be largest.

The years 2009 and 2010 show a transition from a negative (−4.6 %) to a positive (4 %) VF anomaly in the Sahel and are therefore suitable to assess the potential contribution of the changing land surface to this transition as compared to a climatological case. The regional patterns of \(\Delta {\mathrm {Y}}_{\mathrm {Dyn}}\) in Fig. 5 show that differences in VF locally reach over 10 % with more vegetation in the Sahel but lower values in the eastern parts of the Sudanian zone. In the Sahel, which is especially prone to interannual PCP variations and where the bare sand is considerably brighter than vegetation, the albedo locally decreases by more than 3 %. The Aug–Sep spatial correlation of \(\Delta \mathrm {Y}(\mathrm {VF})_{\mathrm {Dyn}}\) with the observed \(\Delta \mathrm {Y}(\mathrm {PCP})\) from June is still about 0.3 and confirms the local lagged co-variation of rainfall and vegetation.

In the following, we investigate the plausibility of the modelled effect of year-to-year vegetation changes for Aug–Sep domain-wide averages with respect to observations. As presented in Table 2, WRF correctly captures more PCP in 2010 than in 2009 in both the DYN and CLIM case. This illustrates that the general monsoon regime is already determined by remote drivers that are fed to the WRF model via the domain boundaries from ERA-I.

DYN shows higher LH because of higher maxima in \({\mathrm {VF}}_{\mathrm {Dyn}}\) (cf. Fig. 2) which reduces the bias for both surface fluxes with respect to the reference data and indicates a more realistic flux partitioning. For both years, the ratio of evapotranspiration to precipitation is 6 % higher than for CLIM and also 2 % higher than for the reference datasets. The sensible heat flux generally remains too high due to an overestimation of incoming solar radiation related to model uncertainties regarding the generation of clouds (not shown).

DYN also decreases the precipitation bias because of less absolute rainfall of about 12 mm month\(^{-1}\) for both years. This distinct change is most likely due to the higher values of \(\mathrm {ALB}_{\mathrm {Dyn}}\) (+1.8 % on average) in comparison to \({\mathrm {ALB}}_{\mathrm {Clim}}\) (cf. Fig. 2) leading to more reflected solar energy (countering the bias in incoming short-wave radiation), overall lower TS of about −0.6 K and smaller SH of 6 W m\(^{-2}\). However, these improvements are related to technical differences in the albedo dataset instead of realistic surface changes and are thus purely artificial.

To exclude the effect of such artificial offsets, it is helpful to compare the change in \(\varDelta Y\), the modelled interannual difference for any variable, instead of comparing the change for a specific year.

The absolute \(\Delta \mathrm {Y}\) in Table 2 is generally larger for DYN than for CLIM, as could be expected under the influence of the two opposing vegetation anomalies during 2009 and 2010.

To get an impression on the significance of the difference between \(\Delta {\mathrm {Y}}_{\mathrm {Dyn}}\) and \(\Delta {\mathrm {Y}}_{\mathrm {Clim}}\) (=\(\Delta {\mathrm {Y}}_{\mathrm {Srfc}}\), cf. Eq. 1), we compare \(\Delta {\mathrm {Y}}_{\mathrm {Srfc}}\) to the spread of \(\Delta \mathrm {Y}\) per ensemble in Fig. 6. Interestingly, the spread of \(\Delta \mathrm {Y}\) for the corresponding reference datasets (blue) is always smaller than the spread of the WRF ensembles, in this case attributing larger uncertainty to the model’s internal variability than to the choice of a certain reference dataset.

For PCP and SH, the WRF ensemble spreads for \(\Delta \mathrm {Y}\) strongly overlap. This means that, for the study region average, instead of changing the surface information, we could have changed the initial conditions of the WRF model to produce differences of the same or an even larger margin for these two variables, rendering \(\Delta {\mathrm {Y}}_{\mathrm {Srfc}}\) insignificant. The spread in PCP increases by 83 % from 12 mm month\(^{-1}\) for \(\Delta {\mathrm {Y}}_{\mathrm {Clim}}\) to 22.3 mm month\(^{-1}\) for \(\Delta {\mathrm {Y}}_{\mathrm {Dyn}}\), suggesting that the dynamical surface considerably increases the internal variability of precipitation generation in the model.

Other than for PCP, the \(\Delta {\mathrm {Y}}_{\mathrm {Srfc}}\) cooling signal of −0.1 K for TS and the related LH (+3.3 W m\(^{-2}\)) is significant (\({P}\le 0.01\)), estimated by a two-tailed t test. In accordance to the observations, DYN shows a positive \(\Delta \mathrm {Y}\)(LH), although overestimates. CLIM does not produce a clear difference in LH between the two years suggesting that the dynamical vegetation improves the models ability to represent the interannual difference of LH in this case. Associated with LH, the decrease in TS is too strong for DYN, indicating a too large sensitivity of surface fluxes on changing vegetation conditions in the WRF model that leads to a larger \(\Delta \mathrm {Y}\) than observed. This overestimated sensitivity of \(\Delta \mathrm {Y}(\mathrm {LH})_{\mathrm {Dyn}}\) to vegetation change might partly be related to an overestimation of the LAI due to the simple upscaling via averaging.

The clear impact of vegetation changes on TS raises the question whether the monsoon dynamics, such as the African Easterly Jet (AEJ), could also be affected. Being a thermal wind, the AEJ follows the largest temperature gradient over the continent that can be modulated by vegetation patterns. Disturbances of the AEJ may trigger or support mesoscale convective systems (MCS) in its proximity and a shift of the jet may therefore change the meridional distribution of precipitation (Cook 1999).

The normalized meridional cross-sections in Fig. 7a reveal a broadened but weaker maximum of the Sahelian temperature gradient between 1 and 15 \(^\circ \)N in 2010 compared to CLIM. The resulting northward shift of the AEJ leads to a better representation of the difference between the 2010/2009 AEJ positions with respect to ERA-I. Accordingly, the precipitation peak in the Sahel for \(\Delta {\mathrm {Y}}_{\mathrm {Dyn}}\) is broadened to the north (Fig. 7b) and directly coincides with the largest \(\Delta {\mathrm {Y}}_{\mathrm {Dyn}}\) for VF and LH (cf. Fig. 7c). This PCP increase is however not visible in the reference datasets. Instead, WRF generally overestimates \(\Delta \mathrm {Y}\)(PCP) in the Sahel and DYN further increases this positive bias via a positive feedback between LH and PCP.

The vegetation-induced change in \(\Delta \mathrm {Y}\)(LH) is directly proportional to \(\Delta \mathrm {Y}\)(VF), as the overlapping curves for \(\Delta {\mathrm {Y}}_{\mathrm {Srfc}}\) in Fig. 7c illustrate. However, such a straightforward relationship does not exist for \(\Delta \mathrm {Y}(\mathrm {PCP})_{\mathrm {Srfc}}\): Although \(\Delta \mathrm {Y}(\mathrm {VF})_{\mathrm {Srfc}}\) is generally positive, there is a negative effect on PCP between 10 and 12 \(^\circ \)N (Fig. 7b, \(\Delta {\mathrm {Y}}_{\mathrm {Srfc}}\)) that we only find in CMORPH amongst the observations. A further examination of \(\Delta {\mathrm {Y}}_{\mathrm {Srfc}}\) in Fig. 7a reveals that in this region, the temperature (and pressure) gradient is weakened by the stronger vegetation gradient, suggesting a relatively lower near-surface moisture convergence that leads to reduced rainfall. Correspondingly, the areas with a positive feedback on PCP are marked by a stronger temperature (pressure) gradient (positive \(\Delta {\mathrm {Y}}_{\mathrm {Srfc}}\)).

This confirms that large-scale differential heating due to vegetation heterogeneities can affect the dynamics and related precipitation even on continental scale. However, in this case, DYN reduces the correspondence of modelled LH and PCP change with the average of the observational datasets, shown by higher mean absolute deviations (MAD) and lower r\(^2\). Again, the change in LH with vegetation and the resulting lower surface temperatures seem to be overestimated, especially in the Sahel, leading to an exaggerated effect on the meridional gradients.

In summary, CLIM is able to capture the domain wide interannual differences of variable averages and meridional gradients as good or better than DYN. While the dynamical vegetation even exerts a clear influence on the modelled monsoon dynamics, the magnitude seems to be overestimated. We therefore conclude that large-scale drivers dominate the observed interannual changes and that a dynamical land surface description does not add clear value at this seasonally averaged spatio-temporal scale within the uncertainty range of observations. It remains to be answered whether the changing vegetation patterns might exhibit a regional or local effect that cannot be captured by CLIM.

The meridional differences in the surface response to vegetation changes are an expected result if we consider two factors that modify the impact of vegetation on spatial LH variability: (1) the evaporative regime and thus the variable that limits evapotranspiration, (2) the actual difference between transpiration and bare soil evaporation.

Regarding the first factor (1), Seneviratne et al. (2010) distinguish between moisture-limited regimes, where soil moisture (SM) changes have a maximum effect on the turbulent surface fluxes, and energy-limited regimes where SM is plenty and these fluxes are limited by incoming solar radiation. Only in the moisture-limited regime, vegetation (root zone soil moisture) can act as a first-order constraint for LH changes, analogous to moist bare soil.

In the moist region of the monsoon rainband however, the frequent cloud cover limits the available energy for evapotranspiration and LH patterns are therefore predominantly controlled (limited) by incoming short-wave radiation, reducing the secondary effect of spatial variations in vegetation. For the investigated time period, the transition zone of the two regimes lies approximately between 12 and 14\(^\circ \)N, as depicted in Fig. 9a. In the Sahel, positive daily correlations between \(\Delta \mathrm {Y}(\mathrm {LH})_{\mathrm {Srfc}}\) and \(\Delta \mathrm {Y}(\mathrm {SM})_{\mathrm {Srfc}}\) illustrate the moisture-limited region, while throughout the Sudanian zone, \(\Delta \mathrm {Y}(\mathrm {LH})_{\mathrm {Srfc}}\) variability is controlled by incoming short-wave radiation.

The second factor (2) that affects the ratio of \(\Delta \mathrm {Y}(\mathrm {LH})_{\mathrm {Srfc}}\) and \(\Delta \mathrm {Y}(\mathrm {VF})_{\mathrm {Srfc}}\) is related to the fact that vegetation exhibits an evaporative advantage compared to bare soil. The advantage exists because the volume of water in the soil that is available for vegetation transpiration is larger than for soil evaporation. This is especially important in the Sahel, where the surface soil layer dries out quickly and where vegetation provides a profound longer-lasting moisture supply to the atmosphere, resulting in stronger LH changes in response to VF changes. This effect is less important in the Sudanian zone, where surface soil moisture is replenished more frequently.

Nevertheless, we also identify an increased sensitivity to \(\Delta \mathrm {Y}(\mathrm {VF})_{\mathrm {Srfc}}\) at the southern border of the domain in Fig. 8, where we find an increase of \(\Delta \mathrm {Y}(\mathrm {LH})_{\mathrm {Srfc}}\) of 1 W m\(^{-2}\) (0.04 mm day\(^{-1}\)) per unit \(\Delta \mathrm {Y}(\mathrm {VF})_{\mathrm {Srfc}}\), comparable to the moisture-limited Sahel. This is because the evaporative advantage also increases when vegetation density increases due to a larger evaporative surface as compared to the bare ground. The canopy is much more dense in the Sudanian zone as illustrated in Fig.  2 with a more than two times higher LAI than in the Sahel. Consequently, a high amount of available root zone soil moisture in dry regimes as well as a high canopy density in wet regimes contribute to a stronger increase of \(\Delta \mathrm {Y}(\mathrm {LH})_{\mathrm {Srfc}}\) per increase of unit \(\Delta \mathrm {Y}(\mathrm {VF})_{\mathrm {Srfc}}\).

Interestingly, \(\Delta \mathrm {Y}(\mathrm {TS})_{\mathrm {Srfc}}\) does not reproduce the behaviour of \(\Delta \mathrm {Y}(\mathrm {LH})_{\mathrm {Srfc}}\) in the Sahel (Fig. 8b). While the decrease in \(\Delta \mathrm {Y}(\mathrm {TS})_{\mathrm {Srfc}}\) per unit \(\Delta \mathrm {Y}(\mathrm {VF})_{\mathrm {Srfc}}\) is indeed somewhat stronger (−0.09 K) in the southern Sudanian zone, the relationship varies around −0.05 K in the rest of the domain. One reason for this is that on smaller time-steps, a net-warming in the vegetated areas may occur when soil moisture is not sufficient to supply transpiration. In these cases, we find a daytime heating effect irrespective of a positive or negative \(\Delta \mathrm {Y}(\mathrm {VF})_{\mathrm {Srfc}}\), consequently weakening the time-averaged overall cooling by a positive \(\Delta \mathrm {Y}(\mathrm {VF})_{\mathrm {Srfc}}\).

In agreement with \(\Delta \mathrm {Y}(\mathrm {LH})_{\mathrm {Srfc}}\), the temporal variability of \(\Delta \mathrm {Y}(\mathrm {TS})_{\mathrm {Srfc}}\) is affected by the predominant evaporative regime. Figure  9b reveals that only north of the transition zone at about \(12^\circ \)N, changes in the \(\Delta {\mathrm {Y}}_{\mathrm {Srfc}}\) of temperature and turbulent heat flux partitioning are predominantly controlled by changes in \(\Delta \mathrm {Y}(\mathrm {LH})_{\mathrm {Srfc}}\) (negatively correlated), depending on the availability of soil moisture.

Opposed to that, the positive correlation between \(\Delta \mathrm {Y}(\mathrm {LH})_{\mathrm {Srfc}}\) and \(\Delta \mathrm {Y}(\mathrm {TS})_{\mathrm {Srfc}}\) in the energy-limited Sudanian zone is a spurious correlation that is actually linked to the dependence of both variables on solar radiation which also applies to the relationship of \(\Delta \mathrm {Y}(\mathrm {LH})_{\mathrm {Srfc}}\) and \(\Delta \mathrm {Y}(\mathrm {SH})_{\mathrm {Srfc}}\) (Fig.  9d). The \(\Delta {\mathrm {Y}}_{\mathrm {Srfc}}\) of surface temperatures and both turbulent fluxes increase simultaneously (positive correlation) with a higher \(\Delta {\mathrm {Y}}_{\mathrm {Srfc}}\) of the incoming solar radiation, explaining the weaker response of \(\Delta \mathrm {Y}(\mathrm {EF})_{\mathrm {Srfc}}\) to \(\Delta \mathrm {Y}(\mathrm {VF})_{\mathrm {Srfc}}\) changes in Fig. 8d.

These two different process pathways of dominant temperature control via soil-moisture anomalies (terrestial control) or via radiation anomalies (atmospheric control) were also identified by Berg et al. (2015) in a comparison of a set of climate simulations with and without interactive soil moisture. They conclude that land–atmosphere feedbacks in energy-limited regions can only play a minor role, since solar radiation determines temporal surface temperature variability. This drives latent heat flux variability which then drives soil moisture variability. This finding is supported by Dirmeyer (2011), who points out that the identification of strong surface–atmosphere coupling (e.g. a high correlation of the surface state and the atmospheric response) is only valid if there is surface (e.g. soil moisture) variability in time.

However, this purely temporal definition of land–atmosphere coupling does not take into account that persistent spatial surface characteristics might also have a considerable effect on the atmosphere. In the case of vegetation, we see from Fig. 8 that, although solar radiation controls the temporal changes of the surface variables in the Sudanian zone, the long-lasting vegetation patterns lead to aggregated spatial characteristics of TS and LH. We assume that such persistent surface gradients may force the atmosphere to reoccurring states just like, for example, a cold lake surface may impose frequent subsidence on the overlying air masses such that precipitation is suppressed. To be able to trace a land–atmosphere feedback for a temporally slowly changing feature like vegetation, we therefore define the vegetation–atmosphere feedback as the spatial correlation of atmospheric characteristics with vegetation patterns over an aggregated time period as opposed to the temporal correlation approach.

We now want to investigate whether the vegetation-driven change in patterns of the surface variables have a local effect on interannual changes in the atmosphere. As explained before, the vegetation patterns are relatively static with respect to typical atmospheric time scales. A land–atmosphere feedback (if any) should therefore manifest itself in a spatial accumulation of certain atmospheric states if the interannual vegetation change is strong enough to force the atmosphere.

The significant spatial correlations between \(\Delta \mathrm {Y}(\mathrm {VF})_{\mathrm {Srfc}}\) and the \(\Delta {\mathrm {Y}}_{\mathrm {Srfc}}\) of the different atmospheric variables (ranging from 0.46 to −0.76, \(P \le 0.01\)) in Fig.  10a–d indeed show that the surface characteristics modify the conditions in the PBL. The correlation coefficients are generally lower than for the surface variables, not only because the atmospheric sensitivity is ultimately determined by large-scale conditions, but also because our method can only capture pixelwise coinciding interactions between VF and the atmosphere. Therefore, we cannot relate any downstream changes in the atmosphere to their origin at the surface, for example.

Another factor for lower spatial correlations is the small-scale spatial heterogeneity of VF as compared to the well-mixed PBL conditions. To remedy this scale discrepancy, VF was upscaled to a length scale of 15 km, which is supposedly at the lower end at which PBL anomalies might persist and impact convection (Clark and Taylor 2004).

In accordance to \(\Delta \mathrm {Y}(\mathrm {TS})_{\mathrm {Srfc}}\), the \(\Delta {\mathrm {Y}}_{\mathrm {Srfc}}\) of the height of the PBL is anti-correlated with the interannual change in VF, which is most pronounced in regions where \(\Delta \mathrm {Y}(\mathrm {VF})_{\mathrm {Srfc}}\) surpasses 10 % (cf. Fig. 8a). In regions with a positive vegetation change \(\Delta \mathrm {Y}(\mathrm {VF})_{\mathrm {Srfc}}\) and corresponding negative \(\Delta {\mathrm {Y}}_{\mathrm {Srfc}}\) for PBL heights, we find an increase of daytime hours in which the lifted condensation level (LCL) lies within the PBL (Fig.  10b). This is because more moisture (higher LH) trapped in a lower PBL increases the relative humidity which ultimately lowers the LCL and leads to cloud formation. Note that this cloud formation process via PBL moistening is less frequent in the Sahel (local increases of 20–40 h), where the air is usually far from saturated due to less moisture advection (see also Findell and Eltahir 2003 for a description of the surface-driven mechanisms of cloud formation via moistening or warming of the PBL).

Although there is a tendency to more frequent cloud cover over regions with positive \(\Delta \mathrm {Y}(\mathrm {VF})_{\mathrm {Srfc}}\) during the day, this does not necessarily imply an increase in precipitation. For example, in the South of Mali, \(\Delta \mathrm {Y}(\mathrm {VF})_{\mathrm {Srfc}}\) is strongly positive and locally causes over +60 h in which the LCL lies within the PBL. However, there is no significant increase in the \(\Delta {\mathrm {Y}}_{\mathrm {Srfc}}\) of the convective available potential energy (CAPE) (Fig. 10a) suggesting that there is no additional destabilization of the atmosphere.

On the contrary, the number of daytime rainy hours is decreased over this region (Fig. 10e), indicating that higher VF fosters the formation of shallow clouds but decreases the potential for initiation of deep convection in the afternoon. Hohenegger et al. (2009) suggest that the development of shallow clouds is an important ingredient for the negative vegetation–precipitation feedback due to higher stability of the air above the PBL related to the longwave radiative cooling at the cloud tops. Garcia-Carreras et al. (2011) analysed convective cloud distributions over crop land and forests with large eddy simulations and attribute the suppressed initiation of convection over forests to a stabilizing capping layer of warm air due to subsidence above the mixed layer. They found this subsidence to be part of mesoscale flows initiated by temperature gradients between the two vegetation classes.

Such mesoscale horizontal flows are also visible in our simulations between cooler and warmer surface patches in Fig. 10c, d, where cooler near surface temperatures and shallower mixing layers correspond to a \(\Delta {\mathrm {Y}}_{\mathrm {Srfc}}\) signal of diverging winds that induce a negative VF-PCP feedback. The converging wind vectors over warmer regions indicate increased moisture convergence from the surroundings that favours the initiation of convection and ultimately increases the number of rainy hours during the day (Fig. 10e).

During the night however, there is a distinct increase of rainy hours over the positive \(\Delta \mathrm {Y}(\mathrm {VF})_{\mathrm {Srfc}}\) in southern Mali (Fig. 10f). This signal is weaker in the rest of the domain but the correlation coefficients in Fig. 10e, f indeed indicate an overall negative (positive) VF-PCP feedback during the day (night). Overall, for 62 % of the pixels (n \(=\) 3990) that show a shift of the diurnal maximum PCP from night (day) to day (night) between 2009 and 2010, this shift coincides with negative (positive) \(\Delta \mathrm {Y}(\mathrm {VF})_{\mathrm {Srfc}}\), which is 12 % more likely than what we would expect by chance (Fig. 10g). The absolute \(\Delta {\mathrm {Y}}_{\mathrm {Srfc}}\) of rainy hours (positive and negative) in the Sahel is 10 % higher during the night than during the day, indicating a slightly stronger effect of the positive VF-PCP feedback on nocturnal rainfall. In the Sudanian zone, the absolute \(\Delta {\mathrm {Y}}_{\mathrm {Srfc}}\) of rainy hours are similar for night and day.

Our findings are in line with Taylor et al. (2010), who analysed a MCS that was observed during the AMMA special observing period (Redelsperger et al. 2006) and report deepest convection over wet surface patches during its mature stage during the night (2130 UTC) while all areas of new convection emerged over drier soils. While the afternoon initiation of new convective cells is favoured over drier soil patches where thermals are more intense, nighttime PCP mostly falls from already existing convective systems that are enhanced by moist surfaces (Gantner and Kalthoff 2010). Wet soils have a stronger effect on the lifetime and strength of the MCS in the Sahel where moisture is limited, which might explain the enhanced positive feedback of vegetation in that region.

Due to the existence of these positive and negative VF-PCP feedbacks, especially in the Sudanian zone, there is no correlation for \(\Delta \mathrm {Y}(\mathrm {PCP})_{\mathrm {Srfc}}\) in Fig. 10h. On average, the significant local changes reach ±2.2 mm day\(^{-1}\), for which the sign largely depends on whether the daytime or the nighttime VF-PCP feedback dominates.

The previous section showed that an atmospheric signal that is related to the interannual vegetation changes can be extracted not only for variables that characterize the state of the PBL, but also for rainfall. The question remains whether the simulated vegetation-driven adjustments are realistic and if so, whether the signal is detectable in the observations for which the dominating large-scale signal cannot be removed i.e. it is not possible to determine a \(\Delta {\mathrm {Y}}_{\mathrm {Srfc}}\).

Nevertheless, the daily PCP statistics in Fig. 11 for 2009 and 2010 show a small but consistent increase of the Heidke skill score (HSS) for DYN for thresholds below 40 mm day\(^{-1}\) during both years. The HSS measures the improvement of the forecast skill over a random forecast and is computed for different thresholds over the whole study domain for every grid cell and every day. The HSS for both WRF cases is relatively low, but it should be kept in mind that WRF is only forced at the domain boundaries (cf. Sect. 2.1) and is otherwise allowed to create its own weather. Reasons why these free running simulations can perform better in capturing precipitating systems in space and time than a random process (\(\hbox {HSS} > 0\)) must therefore be related to strong atmospheric boundary conditions (e.g. atmospheric waves travelling through the model domain) or, to a lesser extent and only in the case of DYN, by the introduced land surface information. Since there is no substantial difference between DYN and CLIM in the frequency distribution or in the ability to capture daily mean PCP, the slightly higher HSS for DYN most likely stems from a surface-driven improvement in the spatio-temporal localisation of precipitating systems.

To test the impact of changes in VF, Fig.  12 shows the strength of the relationship between \(\Delta \mathrm {Y}(\mathrm {VF})_{\mathrm {Dyn}}\) and the \(\Delta \mathrm {Y}(\mathrm {PCP)}\) of DYN, CLIM and the reference datasets. We are only including pixels where \(\Delta \mathrm {Y}(\mathrm {PCP})_{\mathrm {Srfc}}\) is significant (cf. Fig. 10h) presuming that these regions have the potential for considerable land–atmosphere feedbacks (sufficient surface change) and to therefore strengthen the surface signal.

The difficulty here is to extract a vegetation–precipitation feedback from the obvious precipitation–vegetation forcing. We rely on the idea that the correlation between \(\Delta \mathrm {Y}(\mathrm {PCP})_{\mathrm {Clim}}\) and \(\Delta \mathrm {Y}(\mathrm {VF})_{\mathrm {Dyn}}\) represents the baseline relationship between large-scale driven PCP changes and resulting VF patterns. Any correlation surplus (strengthened relationship) for DYN or the reference datasets in comparison to CLIM should be an indicator for an effect of the surface on the PCP patterns.

This approach works well in the Sahel (Fig.  12a) where we find a \(\Delta \mathrm {Y}(\mathrm {VF})_{\mathrm {Dyn}}\) − \(\Delta \mathrm {Y}(\mathrm {PCP})_{\mathrm {Dyn}}\) correlation of 0.23 corresponding to correlations between 0.16 and 0.23 for the reference datasets. This relationship cannot be explained by large-scale interannual variability alone as illustrated by a very small correlation of 0.08 for \(\Delta \mathrm {Y}(\mathrm {PCP})_{\mathrm {Clim}}\) which confirms an improved spatial PCP distribution for DYN.

The correlations of the reference datasets show that \(\Delta \mathrm {Y}(\mathrm {VF})_{\mathrm {Dyn}}\) explains at most 5 % spatial variance (r\(^2\) = 0.23\(^2\) = 0.05 for RFE) of the average PCP changes between Aug–Sep 2009 and 2010. Note that this is valid only for the regions where a potentially strong feedback was detected.

In the Sudanian zone, there is no correlation for the WRF simulations but small positive correlations for the reference datasets most likely due to the precipitation–vegetation forcing (Fig. 12b). WRF has difficulties to capture the observed patterns of large-scale driven PCP differences in that region, which might explain the complete lack of a correlation with \(\Delta \mathrm {Y}(\mathrm {VF})_{\mathrm {Dyn}}\). Ultimately, the strong predominance of the large-scale monsoon dynamics in determining interannual precipitation differences in the Sudanian zone seems to inhibit any signal detection directly from \(\Delta \mathrm {Y}\).