The role of transient eddies and diabatic heating in the maintenance of European heat waves: a nonlinear quasi-stationary wave perspective

Abstract

European heat waves result from large-scale stationary waves and have major impacts on the economy and mortality. However, the dynamical processes leading to and maintaining heat waves are still not well understood. Here we use a nonlinear stationary wave model (NSWM) to examine the role played by anomalous stationary waves and how they are forced during heat waves. For our study, we use the Japanese Reanalysis (JRA-55) data for the period 1958 through 2017. We show that the NSWM can successfully reproduce the main features of the observed anomalous stationary waves in the upper troposphere. Our results indicate that the dynamics of heat waves are nonlinear, and transient momentum fluxes are the primary drivers of the observed anomalous stationary waves. The contribution from orographic forcing is moderate and mainly through nonlinear interactions with diabatic heating. Further decomposition of the transients indicates that the high-frequency transient vorticity fluxes make dominant contributions. Furthermore, our results reveal that the response to heating located in the tropical Indian Ocean and the west Pacific region is primarily responsible for maintaining the observed anomalous stationary waves linked to European heat waves. This is confirmed by exploring the relationship between heat waves and the Indian Ocean Dipole strength. The heating in the mid-latitude and tropical Atlantic region plays a secondary role. Our results suggest that European heat waves are potentially predictable by considering the nonlinear effects involved in anomalous stationary waves and the heating sources in the nearby and remote tropical region.

1 Introduction

Mid-latitude heat waves are expected to become more frequent and intense due to anthropogenic global warming (Perkins 2015; Kornhuber et al. 2019). Over the last few decades, heat waves had significant impacts on human health, terrestrial ecosystems, and global crop production (Horton et al. 2016; Deng et al. 2018). For instance, the 2003 European and the 2010 Russian heat waves produced significant socio-economic damages (Russo et al. 2015). In the context of anthropogenic global warming, heat waves in Europe will become commonplace by the end of the twenty-first century (Fragkoulidis et al. 2018; Carril et al. 2008). Hence, a better understanding of the physical and dynamical processes governing European heat waves is essential for improved predictability and adaptation measures.

In general, European heat waves are associated with persistent large-scale circulation anomalies, i.e., anti-cyclones or blocking highs (Stefanon et al. 2012; Woollings et al. 2018), which interrupt the prevailing westerlies from the North Atlantic and are associated with fair weather and clear skies, bringing an increase of radiative forcing, and then leading to temperature rises near the surface. For example, Zschenderlein et al. (2019) investigated the processes determining European heat waves using Lagrangian analysis and found that co-located anti-cyclones are quite important in establishing high temperatures near the surface due to the subsidence of air. Wolf et al. (2018) demonstrated that European heat waves have a clear connection with localized blocking anti-cyclones, which are also zonally elongated. However, the dynamical mechanisms involved in persistent anti-cyclonic circulations are still unclear.

The identification of the drivers of persistent anti-cyclones is complicated due to the complexity of the dynamical and physical processes involved. Blocking events can be attributed to internal atmospheric dynamics, but also can be triggered or forced by remote anomalous heating over tropical regions (Dole et al. 2011; Park and Lee 2019). Recent global climate change also contributes to the difficulty in examining the mechanisms, owing to the strongly nonlinear feedbacks in the thermodynamic and physical components of the climate variability (Nakamura and Huang 2018; Schär et al. 2004; Risbey et al. 2018).

Over Europe, persistent blocking anti-cyclones are part of slow-moving planetary waves or quasi-stationary Rossby waves (Kornhuber et al. 2017). These quasi-stationary waves are the primary drivers of the mid-latitude weather variations on weekly or monthly timescales, and their activities further impact on the position of the jet streams and act as guides for transient perturbations such as storms and blocking events (Sellevold et al. 2016; Held et al. 2002). Petoukhov et al. (2016) and Kornhuber et al. (2019) illustrate that persistent stationary waves can be explained by a phase-locked quasi-resonant amplification. This quasi-resonant amplification leads to persistent extreme circulation events and temperature anomalies through circumglobal wave patterns as well. Large-scale teleconnection patterns are also associated with stationary waves (Feldstein and Franzke 2017; Hannachi et al. 2017; Franzke and Feldstein 2005). For example, Luo et al. (2007) points out that the decay of the positive phase of North Atlantic Oscillation (NAO) events can result in blockings over Europe via enhanced downstream energy dispersion. Li et al. (2020) further demonstrates that the atmospheric circulation during European blocking and positive NAO events leads to a favorable circulation pattern for heat waves over northern and western Europe.

Furthermore, tropical diabatic heating is a major Rossby wave source (Teng and Branstator 2019). Heating from a tropical region has a significant influence on atmospheric blocking conditions which are necessary to cause anomalous summer heat waves over Europe (Della-Marta et al. 2007). The important tropical regions have been identified as the tropical Atlantic (Cassou et al. 2005), the tropical Indian Ocean and the northwest Pacific regions (Behera et al. 2013). Numerical studies by Bader and Latif (2005) also showed that anomalous Indian Ocean sea surface temperatures have an impact on the NAO and Europe based on a coupled GCM, which has a strong connection to the quasi-stationary waves over Europe (Wolf et al. 2018). However, there is no consensus regarding the locations where heating can force persistent blocking systems over Europe. Building on these findings, our study aims to investigate the dynamical mechanism of maintaining anomalous stationary waves and identify the key region of tropical heating related to European heat waves.

Based on linear stationary Rossby wave theory, stationary waves are defined here as departures from the climatological—and zonal—mean state, which arise from three forcing components: topography, diabatic heating, and transient eddy fluxes (Simpson et al. 2016; Ting 1994) and their interactions. Changes in any of the forcing terms or their nonlinear interactions between each other may modify the stationary wave patterns, depending on the magnitude and location of the perturbation of the forcing (Sobolowski et al. 2011). Models of intermediate complexity have been used to successfully explain and understand the forcing and maintenance of stationary waves (Held et al. 2002). In this study we use a fully nonlinear stationary wave model (NSWM) developed by Ting and Yu (1998). While this model has some limitations (Held et al. 2002), it enables us to decompose the response to given forcings to help identify causal relationships among the three forcing terms (Simpson et al. 2016). As a nonlinear diagnostic tool, NSWM has been used widely in investigating the maintenance mechanisms of stationary waves in seasonal or climatological settings (Ting 1994; Simpson et al. 2016), and has also been used in reproducing some basic features of observed waves for specific circulation anomalies (Liu et al. 1998; Schubert et al. 2014). An advantage of the NSWM is that it allows us to diagnose the nonlinear interactions among different forcings in maintaining the stationary waves. To our knowledge, there are no previous studies that provide a dynamic analysis of the nonlinear effects of European heat waves.

Here, we conduct a systematic study of the dynamic processes driving the circulation anomalies during persistent European heat waves using the NSWM. While most previous studies focused on individual heat wave events, here we will examine a large set of European heat wave events. Building upon previous studies, we hypothesize that the meteorological aspects of European heat waves are caused by anomalous stationary waves (blocking highs) in the upper troposphere. Specifically, we seek to address the following questions: (1) What are the significant forcing terms of the observed anomalous stationary waves? (2) What are the relative contributions from the different nonlinear interactions to the observed stationary waves? and (3) Which tropical region contributes most in forcing anomalous stationary waves associated with European heat waves? While Ghosh et al. (2017, 2019) found that the heating of North Atlantic mid-latitudes related to the positive Atlantic Multi-decadal Oscillation sea surface temperature is important for European heat waves, in this present paper we find that the heating in the tropical Indian Ocean and the west Pacific region is also important for maintaining a blocking anticyclone over Europe linked to European heat waves. The remainder of this paper is organized as follows: Sect. 2 describes the data, the NSWM and the experimental design, and the decomposition method of the nonlinear interactions of the NSWM. The results of the identification of persistent European heat wave events and our diagnostic analysis are presented in Sect. 3. The conclusions are given in Sect. 4.

2 Data, methods, and model experiments

2.1 Data source

For our study we obtain the reanalysis data from the 6-hourly JRA-55 (Japanese 55-year Reanalysis) data-set (Kobayashi et al. 2015), with a horizontal resolution of ${1.25{^\circ}} \times {1.25{^\circ}}$. We use the JRA-55 data for the period from 1958 to 2017, i.e., 60 years in total. Daily mean fields are obtained from these 6-hourly fields. We limit the analysis to the extended boreal summer period from June through August (JJA). The heating rates provided by JRA-55 can be used to calculate diabatic heating rates directly, by adding of the convective, large-scale condensation, solar and long-wave radiation, and vertical diffusion heating rates (Shi et al. 2020). The diabatic heating rate is a key parameter in the NSWM, and the resulting total diabatic heating is consistent with other reanalysis products (Zhang et al. 2017).

We use 2 m daily mean air temperature to identify European heat wave events. To examine large-scale circulation anomalies during heat waves, anomalous stationary waves are obtained in the upper troposphere using 200 hPa streamfunction. For each grid point, the daily climatological mean of the 200 hPa streamfunction field (${\overline{\Psi }}$) is determined by calculating the mean of the 60-year JJA data. The daily zonal mean of ${\overline{\Psi }}$ is subtracted at each grid point to obtain the daily climatological stationary wave $\overline{\Psi ^*}$. Then, the daily stationary wave $\Psi ^*$ is calculated at each grid point by subtracting the zonal mean for each day. Finally, the daily anomalous stationary wave is derived by subtracting $\overline{\Psi ^*}$ from $\Psi ^*$ for the corresponding day of the year. We average the daily anomalous stationary wave during the respective heat wave event period as the observed anomalous stationary wave for that heat wave event. As for the model simulations, we use the three-dimensional wind, temperature, relative vorticity, surface pressure, diabatic heating, and realistic orographic data.

Moreover, to explore the links between European heat waves and large-scale teleconnection patterns and corresponding sea surface temperature anomalies, the daily North Atlantic Oscillation (NAO) index, and the monthly index representing the strength of the Indian Ocean Dipole (IOD) and the NINO3.4 index which describes El Niño and La Niña events are obtained from the National Center for Environmental Prediction/Climate Prediction Center https://www.cpc.ncep.noaa.gov/.

2.2 Definition of heat waves

No universal definition exists for a heat wave event due to the fact that various societal sectors are impacted differently by heat waves (Perkins 2015). Hence, the definition of a heat wave is usually based on the research question at hand (Horton et al. 2016). In general, one can define heat wave events as a period of consecutive days when the daily temperature is above a certain threshold or percentile of the long-term temperature distribution (Stefanon et al. 2012; Deng et al. 2018; Russo et al. 2015).

In this study, we focus on heat wave events in the European region ($35^\circ$ - $60^\circ$ N, $10^\circ$W - $40^\circ$E) over the period 1958– 2017. A heat wave event here has to satisfy two key criteria: (1) the temperature has to be above a threshold, and (2) this threshold has to be crossed for a certain duration. We define that abnormally high temperatures occur when:

$$\begin{aligned} {\widehat{TG}} > {\widehat{TG}}_{90} \end{aligned}$$

(1)

where the hat represents an area-average over the European region, ${\widehat{TG}}$ is the area-averaged daily 2 m air temperature, ${\widehat{TG}}_{90}$ is the 90th percentile of the time series of ${\widehat{TG}}$ over the period 1958–2017.

We identify 80 heat wave events using the above definition (Eq. 1). We change the duration criterion from 3 to 10 days in steps of 1 day. Nearly half of the 80 heat wave events (38 events) can be identified by increasing the threshold to 5 days (Fig. 1). These 38 events are sufficient to represent the dynamical processes behind the 80 events. Also, as stated by Guo et al. (2017), the relaxation of the threshold of the heat wave duration will not modify the heat wave impacts on mortality. Thus, our threshold of 5 days is a good compromise between the duration, to satisfy our assumption of a stationary wave, and the number of heat wave events to analyze. Therefore, we conduct model simulations only focused on these 38 events to facilitate computation. Sensitivity experiments for heat wave events with a longer duration threshold revealed similar results. From this, we derive a heat wave list in Sect. 3.1, with corresponding start dates, end dates, and durations for each heat wave event.

Moreover, we measure the averaged magnitude of heat wave events by using the daily heat wave magnitude index (HWMId) from Russo et al. (2015). This index can measure the magnitude of a heat wave event at each grid point based on the corresponding inter quartile range of TG (Fragkoulidis et al. 2018). We averaged this index for each identified heat wave over the European region for easy comparison due to the different durations among the identified events.

Fig. 1

Number of heat wave events identified by using different duration thresholds

2.3 Model and experimental design

We focus on the maintenance mechanisms of anomalous stationary waves associated with European heat waves. The nonlinear stationary wave model (NSWM) by Ting and Yu (1998) is used to simulate the observed anomalous stationary waves and diagnose the various nonlinear effects in maintaining anomalous stationary waves. This idealized model is a time-dependent baroclinic model, which solves the three-dimensional nonlinear primitive equations for deviations from a prescribed zonal mean basic state in response to zonally asymmetric imposed forcings in $\sigma$-coordinates ($P/P_s$), where P denotes the pressure level and $P_s$ the surface pressure. The horizontal resolution is rhomboidal truncation R30 (about $2.25^\circ \times 3.75^\circ$) with 14 unevenly spaced vertical $\sigma$ levels. The spectral form of the primitive equations has four prognostic equations for vorticity ($\zeta$), divergence (D), temperature (T), the logarithm of surface pressure ($\ln P_s$ ) and two diagnostic equations for the geopotential ($\Phi$) and vertical velocity (${\dot{\sigma }}$):

$$\begin{aligned} \frac{\partial \zeta }{\partial t}&= -\nabla \cdot \left\{ (f+\zeta ) \cdot {\mathbf {V}} \right\} -{\mathbf {k}} \cdot \nabla \nonumber \\&\quad \times \Big \{RT \nabla \ln P_s + {\dot{\sigma }}\frac{\partial {{\mathbf {V}}}}{\partial {\sigma }}\Big \} - \varepsilon \zeta - \nu \nabla ^4\zeta \end{aligned}$$

(2a)

$$\begin{aligned} \frac{\partial D}{\partial t}&= {\mathbf {k}} \cdot \nabla \times \left\{ (f+\zeta ) \cdot {\mathbf {V}} \right\} \nonumber \\&\quad - \nabla \cdot \Big \{RT \nabla \ln P_s + {\dot{\sigma }}\frac{\partial {{\mathbf {V}}}}{\partial {\sigma }}\Big \}\nonumber \\&\quad - \nabla ^2 \Big \{\frac{1}{2}(u^2 + v^2)+\Phi \Big \} \nonumber \\&\quad - \varepsilon D - \nu \nabla ^4 D \end{aligned}$$

(2b)

$$\begin{aligned} \frac{\partial {T}}{\partial {t}}&= - {\mathbf {V}} \cdot \nabla T - {\dot{\sigma }}\frac{\partial {T}}{\partial {\sigma }} + \frac{\kappa T}{\sigma }{\dot{\sigma }} \nonumber \\&\quad +\kappa T \left\{ ({\mathbf {V}} - \tilde{{\mathbf {V}}}) \cdot \nabla \ln P_s - {\tilde{D}} \right\} - \varepsilon T - \nu \nabla ^4 T\end{aligned}$$

(2c)

$$\begin{aligned} \frac{\partial {\ln P_s}}{\partial {t}}&= - \tilde{{\mathbf {V}}} \cdot \nabla \ln P_s - {\tilde{D}} \end{aligned}$$

(2d)

$$\begin{aligned} \frac{\partial \Phi }{\partial \sigma }&= -\frac{RT}{\sigma } \end{aligned}$$

(2e)

$$\begin{aligned} \frac{\partial {\dot{\sigma }}}{\partial \sigma }&= -({\mathbf {V}} - \tilde{{\mathbf {V}}})\cdot \nabla \ln P_s - (D - {\tilde{D}}) \end{aligned}$$

(2f)

where the tilde denotes the vertically averaged quantities, $f = 2\Omega \sin \phi$ the Coriolis parameter, $\Omega$ the angular velocity of the Earth, $\phi$ the latitude, ${\mathbf {V}}$ the 2-D horizontal wind field, and R the specific gas constant.

To describe the evolution of the anomalous variables from a climatological state, the NSWM decomposes all model variables into the climatological state (represented by an over-bar) and the deviation (represented by a prime), and then removes the climatological state equation from Eq. 2. The resulting anomaly equations are (the primes have been omitted in the follow equations for brevity):

$$\begin{aligned} \frac{\partial \zeta }{\partial t}&= -\nabla \cdot \left\{ (f+{\bar{\zeta }}) \cdot {\mathbf {V}} + \zeta \bar{{\mathbf {V}}} + \zeta {\mathbf {V}} \right\} \nonumber \\&\quad - {\mathbf {k}} \cdot \nabla \times \Big \{RT \nabla \overline{\ln P_s} + R {\bar{T}} \nabla \ln P_s + RT\nabla \ln P_s\nonumber \\&\quad + {\dot{\sigma }}\frac{\partial {\bar{{\mathbf {V}}}}}{\partial {\sigma }} +\bar{{\dot{\sigma }}}\frac{\partial {{\mathbf {V}}}}{\partial {\sigma }} + {\dot{\sigma }}\frac{\partial {{\mathbf {V}}}}{\partial {\sigma }} \Big \} + TF_{V} - \varepsilon \zeta - \nu \nabla ^4\zeta \end{aligned}$$

(3a)

$$\begin{aligned} \frac{\partial D}{\partial t}&= {\mathbf {k}} \cdot \nabla \times \left\{ (f+{\bar{\zeta }}){\mathbf {V}} + \zeta \bar{{\mathbf {V}}} + \zeta {\mathbf {V}} \right\} \nonumber \\&\quad -\nabla \cdot \Big \{RT \nabla \overline{\ln P_s} + R {\bar{T}} \nabla \ln P_s + RT\nabla \ln P_s \nonumber \\&\quad +{\dot{\sigma }}\frac{\partial {\bar{{\mathbf {V}}}}}{\partial {\sigma }} +\bar{{\dot{\sigma }}}\frac{\partial {{\mathbf {V}}}}{\partial {\sigma }} + {\dot{\sigma }}\frac{\partial {{\mathbf {V}}}}{\partial {\sigma }} \Big \} \nonumber \\&\quad - \nabla ^2 \Big \{\frac{1}{2}(u^2 + v^2)+ u {\bar{u}} + v {\bar{v}} + \Phi \Big \} \\&\quad + TF_{D} - \varepsilon D - \nu \nabla ^4 D \nonumber \end{aligned}$$

(3b)

$$\begin{aligned} \frac{\partial {T}}{\partial {t}}&= - \bar{{\mathbf {V}}} \cdot \nabla T - {\mathbf {V}} \cdot \nabla {\bar{T}} - {\mathbf {V}} \cdot T - \bar{{\dot{\sigma }}}\frac{\partial {T}}{\partial {\sigma }} \nonumber \\&\quad -{\dot{\sigma }}\frac{\partial {{\bar{T}}}}{\partial {\sigma }} - {\dot{\sigma }}\frac{\partial {T}}{\partial {\sigma }} + \frac{\kappa T}{\sigma }\bar{{\dot{\sigma }}} + \frac{\kappa {\bar{T}}}{\sigma }{\dot{\sigma }} + \frac{\kappa T}{\sigma }{\dot{\sigma }}\nonumber \\&\quad + \kappa T ({\mathbf {V}} - \tilde{{\mathbf {V}}}) \cdot \nabla \ln P_s + \kappa {\bar{T}} ({\mathbf {V}} - \tilde{{\mathbf {V}}}) \nonumber \\&\quad \cdot \nabla \ln P_s + \kappa T ({\mathbf {V}} - \tilde{{\mathbf {V}}}) \cdot \nabla \overline{\ln P_s} \nonumber \\&\quad + \kappa {\bar{T}} (\bar{{\mathbf {V}}} - \bar{\tilde{{\mathbf {V}}}}) \cdot \nabla \ln P_s + \kappa {\bar{T}} ({\mathbf {V}} - \tilde{{\mathbf {V}}}) \nonumber \\&\quad \cdot \nabla \overline{\ln P_s} + \kappa T (\bar{{\mathbf {V}}} - \bar{\tilde{{\mathbf {V}}}}) \cdot \nabla \overline{\ln P_s} \nonumber \\&\quad + \kappa T (\bar{{\mathbf {V}}} - \bar{\tilde{{\mathbf {V}}}}) \cdot \nabla \ln P_s - \kappa {\bar{T}} {\tilde{D}} - \kappa T \bar{{\tilde{D}}}\nonumber \\&\quad - \kappa T {\tilde{D}}+ TF_{TEMP} - \varepsilon T - \nu \nabla ^4 T \end{aligned}$$

(3c)

$$\begin{aligned} \frac{\partial {\ln P_s}}{\partial {t}}&= - \bar{\tilde{{\mathbf {V}}}} \cdot \nabla \ln P_s - \tilde{{\mathbf {V}}} \cdot \nabla \overline{\ln P_s}\nonumber \\&\quad -\tilde{{\mathbf {V}}} \cdot \nabla \ln P_s - {\tilde{D}} - TF_{P_s} \end{aligned}$$

(3d)

$$\begin{aligned} \frac{\partial \Phi }{\partial \sigma }&= -\frac{RT}{\sigma } \end{aligned}$$

(3e)

$$\begin{aligned} \frac{\partial {\dot{\sigma }}}{\partial \sigma }&= -(\bar{{\mathbf {V}}} - \bar{\tilde{{\mathbf {V}}}})\cdot \nabla \ln P_s - ({\mathbf {V}} - \tilde{{\mathbf {V}}})\cdot \nabla \overline{\ln P_s}\nonumber \\&\quad -({\mathbf {V}} - \tilde{{\mathbf {V}}})\cdot \nabla \ln P_s - (D - {\tilde{D}}) \end{aligned}$$

(3f)

where $TF_{TEMP}$ denotes the transient heat flux convergence terms, $TF_{Ps}$ the surface pressure term, $TF_{V}$ the transient vorticity forcing term, and $TF_{D}$ the transient divergence forcing term. The orographic forcing term is in the hydrostatic equation in Eq. 3e as the lower boundary condition (Ting 1994), and also involved in the transient divergence equation in Eq. 3b through geopotential height. Furthermore, various damping terms including Rayleigh friction and Newtonian cooling ($\varepsilon$), and bi-harmonic diffusion ($\nu$) are used in the NSWM to prevent baroclinic instability and eliminate small-scale noise. The model coefficients used in this study are identical with those used in Held et al. (2002).

For the simulations in our study, all the NSWM inputs are from JRA-55. The basic state is a 60-year-average (1958–2017) zonal mean flow consisting of u, v, T, and $P_s$, first interpolated onto the model $\sigma$ levels before taking the zonal mean. The period of the basic state is event-dependent; that means that for each heat wave event we use the 60-year climatological mean state corresponding for that event period. For example, if an identified heat wave lasts from 5th to 15th June 2003, the basic state will be chosen as 5th–15th June from the 60-year-averaged zonal mean flow. Therefore, the NSWM simulates only the response of the zonally asymmetric components, which means the zonal-mean anomalies generated by zonal-mean forcing terms are excluded (Liu et al. 1998). We also tested other definitions of the mean state such as using the climatological JJA mean or the JJA mean in the respective year of the heat wave and the two week period of the respective heat wave. Using those mean states produced very similar results (not shown). This suggests that our results are robust concerning the particular choice of the mean state.

There are three different forcing components for the NSWM in simulating anomalous stationary waves: diabatic heating (Q), transient momentum fluxes (M), and orography (O). The $TF_{TEMP}$ is included in the diabatic heating forcing following Simpson et al. (2016), given the close linkage between the transient heat fluxes and the tendency for the transients to act diffusively on temperature gradients induced by other diabatic sources. We regard the $TF_{V}$ and $TF_{D}$ as transient momentum fluxes forcing, which are computed as in Ting and Yu (1998):

$$\begin{aligned} TF_{V} = -\nabla \cdot (\overline{V'\zeta '}) \end{aligned}$$

(4)

and

$$\begin{aligned} TF_{D} = {\mathbf {k}} \cdot \nabla \times (\overline{V'\zeta '}) - \frac{1}{2}{\nabla{^2}} \overline{(V' \cdot V')} \end{aligned}$$

(5)

where V is the horizontal wind. And the realistic orographic data is from JRA-55. These forcings are then linearly interpolated onto the $\sigma$-coordinates of the NSWM.

We simulate the observed anomalous stationary wave using a time step of 30 min. The NSWM reaches a (quasi-) steady-state after being integrated for around 20 days (Held et al. 2002). Therefore, we take the day 31–50 average for each simulation to ensure that our results are insensitive to the choice of the averaging period. To easily compare with reanalysis data, we also interpolate model results from the $\sigma$-coordinate to pressure coordinates.

As discussed in Sect. 2.2, we first simulate the observed anomalous stationary wave for each of the 38 heat wave events, and the NSWM is driven by all the forcing terms for each simulation. These simulations for each heat wave event refer to the model sensitivity analysis (see Sect. 3.2.1). Then, the composite analysis method (Boschat et al. 2016) is used to determine the robust mechanism of the maintenance of the anomalous stationary wave. We composite the observed anomalous stationary wave for the 38 events, and also the corresponding forcing terms to run the NSWM as the composite experiment.

2.4 Decomposition of the forcing terms in the NSWM

The NSWM enables us to perform simulations with just a subset of the forcings to systematically identify the relevant forcings and their interactions. Following Ting et al. (2001), the total nonlinear response of the NSWM to the three forcings could simply be due to the nonlinear effect of an individual forcing or could be caused by more complex interactions due to various combinations of the three forcing terms and their nonlinear interactions. Hence, we need to quantify these complex mechanisms arising from these forcings in the NSWM. To this end, we use the factor separation method (Stein and Alpert 1993; Cleveland et al. 2020) to calculate the NSWM response to individual forcings and the nonlinear interactions among the different forcings. Ting et al. (2001) and Sobolowski et al. (2011) have already used this method to decompose the results of the NSWM.

Based on the factor separation method, we can decompose the total nonlinear effect (TNE) between two forcings into individual nonlinear effects (INE) and total effects of nonlinear interactions (TENI):

$$\begin{aligned} \text {TNE}_{f_i, f_j} = \text {INE}_{f_i} + \text {INE}_{f_j} + \text {TENI}_{f_i,f_j} \end{aligned}$$

(6)

where the subscripts $f_i$ and $f_j$ represent two different forcings, $\text {TNE}_{f_i, f_j}$ is the total nonlinear effect of the forcings $f_i$ and $f_j$, $\text {INE}_{f_i}$ and $\text {INE}_{f_j}$ are the individual nonlinear effects of forcings $f_i$ and $f_j$ respectively, and $\text {TENI}_{f_i,f_j}$ is the total effect of nonlinear interactions between forcing $f_i$ and $f_j$. Using Eq. 6, we can quantify the total effect of nonlinear interactions between two forcings ($\text {TENI}_{f_i,f_j}$). Then, the nonlinear interactions between three forcings can be derived as:

$$\begin{aligned} \text {TENI}_{f_i,f_j,f_k}= & {} \text {TNE}_{f_i,f_j,f_k} - (\text {INE}_{f_i} + \text {INE}_{f_j} + {{\text {INE}}_{f_k}} \nonumber \\&+ \text {TENI}_{f_i,f_j} + \text {TENI}_{f_i,f_k} +\text {TENI}_{f_j,f_k}) \end{aligned}$$

(7)

And if we substitute Eq. 6 into Eq. 7, we can get the individual total nonlinear effect ($\text {ITNE}$) of a given forcing, for example, $f_k$. It can be obtained as

$$\begin{aligned} \text {ITNE}_{f_k} = \text {TNE}_{f_i, f_j, f_k} - \text {TNE}_{f_i, f_j} \end{aligned}$$

(8)

In this way, various mechanisms of the nonlinear responses in the NSWM can be qualitatively examined. For our simulations, the forcing input in the NSWM are Q, M, and O. Thus, to get the total linear nonlinear effect of three (TNE_QMO) or two forcings (TNE_QM, TNE_QO, and TNE_MO), the model is run with three or two forcing components. To get the individual nonlinear effects (INE_Q, INE_M, and INE_O), the model is driven by only one forcing. All the total effects of nonlinear interactions (TENI_QM, TENI_QO, TENI_MO, and TENI_QMO) are decomposed following Eq. 6 and Eq. 7, as also for the individual total nonlinear effects (ITNE_Q, ITNE_M, and ITNE_O).

Using reanalysis data for our stationary wave analysis has the advantage that our model results can directly be compared with the reanalysis data. However, the reanalysis data might have some dynamical inconsistencies with the NSWM dynamics (Wang and Ting 1999). A comprehensive comparison of the numerical model and the reanalysis, therefore, is difficult because of the brevity of the reanalysis data set (Teng et al. 2013). Here, an area-weighted spatial pattern correlation (hereafter referred to as pattern correlation) (Ting et al. 2001) is used to determine how well the NSWM simulates the spatial structure of the observed anomalous stationary wave. This pattern correlation can further be used to examine the relative contributions of the various mechanisms to maintain the observed anomalous stationary wave while we analyze the decomposed results. The weights are the square of the cosine of latitudes. Note that the pattern correlation measures the similarity between observed and simulated spatial patterns of anomalous stationary wave without considering the magnitude of each pattern. Besides, the statistical significance test is conducted for the pattern correlation following Walsh (2008).

3 Results and discussion

3.1 Identification of heat waves

We devote our attention to the 38 heat wave events which last at least 5 days. Table 1 describes the characteristics of these heat wave events, including dates, duration, region-averaged temperature anomalies ($T_{ano}$), and heat wave magnitudes (HWMId). These identified heat waves are in general consistent with recent studies (Petoukhov et al. 2013; Schubert et al. 2014). Specifically, 29 out of 38 (76.3%) events occurred after 2000, which is likely due to anthropogenic global warming. Furthermore, the frequency of heat waves occurring in one year is increasing, for example, 3 events in 2003, and 4 events in 2012 and 2015 respectively. Similarly, the increasing trend is also evident for characteristics of duration, $T_{ano}$, and HWMId. These increasing trends over Europe may be attributed to regional manifestations of global warming or changes in land-use and associated reductions in soil moisture (Behera et al. 2013). According to $T_{ano}$ and HWMId, the heat waves in 2007, 2010 and 2016 were the most extreme in the last 60 years. Research also shows that these heat wave events were the most extreme over the last 500 years (Casty et al. 2005; Barriopedro et al. 2011).

The composite of surface temperature anomalies (Fig. 2) for the 38 heat wave events clearly shows that the positive $T_{ano}$ over our defined European region (green box in Fig. 2a) has a magnitude of up to 4 K. This is the composite mean, individual events can reach higher anomalous temperatures. The composite of HWMId exhibits identical spatial patterns with the temperature anomaly as well. By examining the corresponding large-scale circulation conditions, we composite the observed anomalous stationary wave for each event in Fig. 2b. An upper-level anti-cyclonic system is in phase with the positive temperature anomalies over Europe. This system favors the occurrence and persistence of regional heat waves (Deng et al. 2018). The persistent blocking location over Europe blocks the weather disturbances from the Atlantic which could alleviate local summer high temperatures (Behera et al. 2013), and strengthens the advection of warm air masses from northern Africa and the Mediterranean basin (Cassou et al. 2005). Figure 2b is also consistent with Li et al. (2020), in which they find that European heat waves are corresponding more to a typical European blocking during positive NAO events. We confirm that European heat waves are dominated by an anti-cyclonic circulation aloft. Thus, we conduct a further numerical investigation for the maintenance of this dynamic linkage in the following.

Large-scale teleconnection patterns play an important role in the modulation of European temperatures, especially the NAO, and the atmospheric response to El Niño and IOD events (Della-Marta et al. 2007; Cassou et al. 2005; Schneidereit et al. 2012). Table 1 shows the values of the corresponding indices of the large-scale teleconnection patterns during each heat wave period. For the NAO index, we averaged the daily NAO index for the duration of the heat wave events; The monthly value of NINO3.4 and IOD is shown for the month when the heat wave occurred. Generally, European heat waves tend to occur during the positive IOD phase (30 out of 38 events), negative NAO period (23 out of 38), and positive NINO3.4 period.

We also examine this relationship by two-dimensional probability distribution function (PDF) plots of $T_{ano}$ and each corresponding indices for the summer season from 1958 to 2017 and the identified heat wave periods, respectively (Fig. 2c–e). We estimate the PDFs in a non-parametric way by using a kernel density estimator (Silverman 1986). There is a clear positive linear relationship between $T_{ano}$ and IOD when all days are considered. But for a given value of IOD, $T_{ano}$ is much stronger on heat wave days than it normally is, indicating that factors other than IOD must also contribute to persistent heat waves (Teng et al. 2013). There is a slight negative relation between $T_{ano}$ and NINO3.4 for all summer months. There is no evident relationship between NAO and $T_{ano}$ but high-temperature anomalies have a slight tendency for negative NAO events.

Our results are consistent with Behera et al. (2013) who also explored the relationship between European heat waves with the IOD and NINO3.4 indices. Their research demonstrated that IOD events have a significant influence on summer heat waves over western Europe. They also indicated El Niño and La Niña events play a role in the extreme Eastern European events. The role of El Niño and La Niña in European heat waves is also examined in Trenberth and Fasullo (2012) and Schneidereit et al. (2012). They found that the 2010 Russian heat wave had a strong link to La Niña conditions which is also shown in Table 1. The linkage of IOD, El Niño and La Niña events reveals the influence of the tropics on European summer heat waves. This motivates us to further examine the role of different tropical regions in forcing and modulating heat waves in the following section.

However, our results indicate that heat waves have no clear relationship with the NAO, which is inconsistent with previous studies like Li et al. (2020), Behera et al. (2013) and Folland et al. (2009). However, there are distinct differences between our study and the studies mentioned above. Our area of heat waves is in central Europe which has significant positive correlations with the summer NAO in the north and significant negative correlations in the south (not shown). These opposite correlations between north and south in our heat wave area likely causes the weak relationship with the NAO. The different role of the NAO for different areas is also found in Behera et al. (2013), where they investigated heat waves for western and eastern Europe and found the role of the NAO is not obvious for eastern Europe while it has a weaker influence on western European summers. The whole of western Europe and part of eastern Europe in their research are contained in our defined region, which could also result in this disagreement. Besides, the composite spatial pattern of the observed anomalous stationary wave (Fig. 2b) shows resemblance to the combination of the positive NAO and the Atlantic low regime rather than solely the positive NAO as shown in Cassou et al. (2005).

Table 1 The selected heat wave events and their features over Europe (35$^\circ N$–60$^\circ N$ and 10$^\circ W$–40$^\circ E$) with corresponding indices of the NAO, NINO3.4, and IOD indices. Dur. is the duration of a heat wave event, $T_{ano}$ and HWMId are the regional averaged surface temperature anomaly and heat wave magnitude over land-surface, respectively

Fig. 2

a Composite spatial pattern of temperature anomalies $T_{ano}$ (shading) and the heat wave magnitude HWMId (contour); the green dashed box depicts the used European region of our study. b The corresponding composites of the 200 hPa streamfunction anomaly; contour interval is $1 \times 10^6\,\,\mathrm{{m}}^2\,\mathrm{{s}}^{-1}$. Probability density functions of daily Europe $T_{ano}$ versus NAO c, NINO3.4 d and IOD e index for JJA (blue) and the heat wave events stated in Table 1 (red)

3.2 NSWM responses to heat waves

3.2.1 Sensitivity analysis of NSWM

The identified 38 heat wave events have different durations (Table 1). Correspondingly, the performance of the NSWM may vary depending on the duration. As stated in section 2.3, we perform a composite experiment using the NSWM with composite forcing terms for all 38 events. To test the robustness of our results, we carry out a sensitivity analysis for the performance of the NSWM. We first run the NSWM for each event with all the anomalous forcings, that is, global diabatic heating (Q), transient momentum fluxes (M), and orography (O). The pattern correlation is used to examine the reproducibility of the NSWM in simulating the observed anomalous stationary wave and is calculated for each simulation over Europe and in the mid-latitudes ranging from $25^\circ$N to $75^\circ$ N. After that, we divide these 38 events into different groups based on their duration. The duration for each group is extended from 5 to 10 days. We then compare the spread, mean and median value of the pattern correlations for these groups.

For the pattern correlations over Europe, variations of the pattern correlations tend to be rather stable for all duration thresholds (Fig. 3). As regards the mid-latitudes, the mean and median are slightly lower than those over Europe, but the values remain as stable as for over Europe. Most of the pattern correlations are statistically significant. Our results suggest that the performance of the NSWM for these 38 events is independent of the duration threshold in terms of the mean and median. The pattern correlation for the composite experiment is about 0.50 over Europe and larger than the significance threshold. Therefore, results from the composite experiment are robust for interpreting and analyzing the dynamic mechanisms in the NSWM.

Fig. 3

Sensitivity analysis using pattern correlations of the NSWM based on different heat wave thresholds. The NSWM is forced by all the forcings: global diabatic heating (Q), transient momentum fluxes (M), and orography (O). The red (blue) dashed line represents the threshold of the statistically significant positive (negative) pattern correlation. Above (below) the red (blue) dashed line indicates statistically significant values. The filled circles represent the mean value for each pattern correlations over all events

3.2.2 NSWM responses to total forcings

We use the NSWM for simulating the observed anomalous stationary wave (Fig. 2b) and diagnosing various nonlinear effects involved in maintaining the anomalous stationary waves. Before decomposing the different contributions, we have to verify that the NSWM is capable of reproducing the observed anomalous stationary wave (Fig. 2b). As stated above, the pattern correlation (0.50), between the simulated anomalous stationary wave from the NSWM forced by all forcing terms and the observed anomalous stationary wave, indicates that the NSWM faithfully reproduces the observed upper-level anomalous stationary wave in the composite experiment. Some differences are to be expected given the simplifications of the nonlinear model (Sobolowski et al. 2011), and also the composite forcing terms might miss some individual forcing information.

To corroborate the consistency between the NSWM and the reanalysis, we compare spatial patterns of the simulated and observed anomalous stationary wave (Figs. 4a and 2b). An obvious anti-cyclonic system is captured over Europe with a relatively lower amplitude compared with the observed anomalous stationary wave. This indicates that the streamfunction response of the NSWM subjected to all three forcings reproduces the robust features of the observed anomalous stationary wave during heat wave events well over Europe by visual inspection. The simulated result in Fig. 4a presents also European blocking with a weaker positive NAO pattern over Europe, in agreement with Li et al. (2020). However, clear discrepancies exist in simulated responses over the North Atlantic, for example, the NSWM fails to reproduce the cyclonic system located in the North Atlantic in Fig. 2b. These discrepancies in simulating the observed anomalous stationary wave may be attributed to the inaccuracy of the dissipation parameterization or other missing physical processes, such as land-ocean-atmosphere interactions (Liu et al. 1998; Ting et al. 2001; Behera et al. 2013). The deficiency of the NSWM in accurately modeling NAO events could also be responsible for these discrepancies (not shown). Overall, the NSWM can provide us with plausible results to simulate the anti-cyclonic system over Europe. To better understand the contributions from the different components, we will further decompose the total stationary nonlinear effects based on the method described in Sect. 2.4.

Fig. 4

The simulated anomalous stationary waves streamfunction fields of the NSWM. a NWSM forced by all the forcings as depicted in Fig. 3. b-d The NSWM is forced by different combinations of two terms in Q, M, and O. The contour interval is $1 \times 10^6\,\,\mathrm{{m}}^2\,\mathrm{{s}}^{-1}$

3.3 Decomposing the processes in NSWM

3.3.1 Combined forcing effects

We first decompose the NSWM responses to two combined forcings as stated in Sect. 2.4, i.e., diabatic heating and transient momentum fluxes (TNE_QM), diabatic heating and orography (TNE_QO), and transient momentum fluxes and orography (TNE_MO). This means we input two forcing terms in each simulation, thus the response of the NSWM not only exhibits the response to individual forcings but also includes their nonlinear interactions. The NSWM streamfunction responses due to the combined forcing mechanisms are shown in Fig. 4b–d. Generally, all the nonlinear responses generated from the combined forcings display a similar spatial pattern to the simulated anomalous stationary wave and all of them show a clear anti-cyclonic system over Europe as the observed anomalous stationary wave (Fig. 2b). However, the response to transient momentum fluxes and orography exhibits the largest response (Fig. 4d), in terms of magnitude, with smaller contributions coming from diabatic heating and orography (Fig. 4c).

Moreover, the center of the anti-cyclonic system is slightly shifted towards the Iberian Peninsula for the response to diabatic heating and orography, but for the other two, the locations of the center are nearly identical and closely resemble the observed and simulated anomalous stationary wave. These differences indicate that transient momentum fluxes and its nonlinear interaction with other forcings make major contributions to the simulation of the observed anomalous stationary wave, as its combinations with the other two forcings have a propensity for generating larger responses and capturing similar centers for the anti-cyclonic system over Europe. Similarly, diabatic heating might play a damping role in the simulations.

The pattern correlations between the NSWM responses to the combined forcing and the observed anomalous stationary wave can provide extra information that examination of the spatial maps might miss (Sobolowski et al. 2011). Results of the correlations confirmed that the contribution from the combined diabatic heating and transient momentum fluxes (0.59) and transient momentum fluxes and orography (0.48) are the two major contributors to the total NSWM response. The two correlation values have passed the significance test. The contribution from the combined diabatic heating and orography plays a smaller role with an insignificant pattern correlation of 0.005. We then investigate which individual mechanisms play an important role in the combined forcings.

3.3.2 Direct nonlinear effects

The NSWM streamfunction responses due to the forcing by the individual mechanisms are shown in the left panel in Fig. 5. These responses represent the direct contributions from individual nonlinear effects to simulate the observed anomalous stationary waves; these are INE_Q, INE_M, and INE_O. Our results reveal that the NSWM response to transient momentum fluxes and orographic forcing have a larger response. In particular, the response of the transient momentum fluxes can generate a consistent anti-cyclonic system compared with the observed anomalous stationary wave (compare Figs. 5b and 2b). As for the responses generated by orography, they show an extended anti-cyclonic system from the North Atlantic to Europe (Fig. 5c). This extended anti-cyclonic system constitutes a downstream flow emanating from the Tibetan Plateau. This is in agreement with Ting et al. (2001) and Held et al. (2002). Orographic forcing effects tend to strengthen the cyclonic flows downstream of the major northern hemispheric mountain chains.

However, the responses of the NSWM to diabatic heating are more complex over Europe and has the smallest magnitude (Fig. 5a). Cyclonic systems are visible over the UK and southern Europe, while an extended anti-cyclonic system reaches and dominates over the Iberian Peninsula. Generally, the role of diabatic heating is somehow different from Ting (1994), where she pointed out that the diabatic heating played a dominant role in maintaining the climatological summer stationary waves. The main reason for this disagreement lies in the fact that the climatological stationary waves during summertime have been removed in our study since we only focus on the anomalous stationary waves. This suggests that diabatic heating plays a major role in driving the zonal mean response.

We then calculate the pattern correlations between the nonlinear responses due to the individual forcing terms and the observed anomalous stationary wave over Europe. Consistent with those described above, transient momentum fluxes are the key contributor to simulate the observed anomalous stationary wave with a significant pattern correlation of 0.52. The pattern correlation for the response to orography is insignificant (0.08), although the responses have the same magnitude with transient momentum fluxes. The pattern correlation reveals that the diabatic heating makes a negative contribution to the total responses of the NSWM due to the significant negative pattern correlation ($-0.18$).

Interestingly, our results suggest that nonlinear interactions between transient momentum fluxes and diabatic heating exist, and should make a significantly positive contribution to the simulation. This is because the pattern correlation for combined transient momentum fluxes and diabatic heating with the observed anomalous stationary wave is 0.59, which is larger than their individual contributions, especially the damping role of diabatic. We will investigate the role of various nonlinear interactions in the following section.

Fig. 5

Same as Fig. 4. (a)–c (left panel) The NWSM is forced by the individual Q, M, and O forcings. d–f (right panel) The NWSM responses to the total nonlinear effects of Q, M, and O

3.3.3 Individual total nonlinear effects

Before analyzing the effect of the nonlinear interactions between the three forcings, we first explore the individual total nonlinear effects to demonstrate whether they play an important role in our simulations. These are expressed by Eq. 8, i.e. ITNE_Q, ITNE_M, and ITNE_O. This individual total nonlinear effect represents the nonlinear response due to the particular forcing and all its interactions with the other forcing terms.

The right panel in Fig. 5 shows the nonlinear response of the NSWM due to the total effects of each forcing term. They should be compared with the corresponding direct nonlinear effects in Sect. 3.3.2. By including the nonlinear interaction effects, the magnitude of both positive and negative responses is increased to some extend. Particularly in the total nonlinear effects of transient momentum fluxes (ITNE_M, Fig. 5e), a clear anti-cyclonic center, which resembles the observed one (Fig. 2b), becomes more evident than due to its direct nonlinear effect (INE_M, Fig. 5b). Considering the total nonlinear effects of orography, the extended anti-cyclonic system over the North Atlantic from the direct nonlinear effects of the orography is split into two centers and the southern center locates around the UK (Fig. 5f). The response due to the total effects of the diabatic heating exhibits a more complex structure with lesser magnitude, enhancing the total response over northwestern and southeastern Europe while mitigating the total response over the southwestern Iberian Peninsula (compare Fig. 5d with Fig. 5a).

When comparing these pattern correlations with their direct nonlinear effects, all the correlation values are increased by including their nonlinear interactions with other forcings. Especially for diabatic heating, the pattern correlation value increases from $-0.18$ to $-0.01$, from significant to insignificant. Due to the lowest correlation, the contribution from the total effects due to diabatic heating is more modest and, as shown in the streamfunction response, more complex than the other two forcings for maintaining the observed anomalous stationary wave. The significant positive pattern correlation of transient momentum fluxes slightly increases from 0.53 to 0.58, and as for the orography, the pattern correlation increases from 0.08 to 0.13 but remains insignificant.

Diagnosis of these pattern correlations combined with the spatial patterns of the streamfunction shown in Fig. 5, suggests that the total nonlinear effects of each forcing are more important than any of the individual direct forcings, especially for the investigation of diabatic heating. We then decompose and examine the role of the various nonlinear interactions between different forcings.

3.3.4 Nonlinear interactions

As shown in Eqs. 6 and 7, there are four nonlinear interaction terms that contribute to the total effects of nonlinear interactions in simulating and maintaining anomalous stationary waves, namely, TENI_QM, TENI_QO, TENI_MO, and TENI_QMO. These nonlinear interactions represent nonlinear interference between flows forced by different forcings. Note that if the dynamic process is completely linear, then the nonlinear interaction terms will be identically zero.

The contributions of each of the interaction terms are shown in Fig. 6. Note that, for display purposes, the range of the contour is 4 times smaller than for the previous figures. The nonlinear interferences between diabatic heating and orography (Fig. 6b) and transient momentum fluxes and orography (Fig. 6c) exhibit clearer spatial patterns with a larger magnitude over the northern hemisphere. These spatial patterns are broadly similar in North America as the downstream extensions emanated from the Tibetan Plateau. However, the responses over Europe have lesser amplitude. This indicates that the flow forced by orography plays an important role in these nonlinear interaction effects in maintaining the anomalous stationary waves. This is consistent with earlier findings (Ting et al. 2001) that emphasize the effect of orography in maintaining climatological and anomalous stationary waves during the NH winter. Our results suggest that the effect of orography is also important in maintaining anomalous stationary waves during heat waves via nonlinear interactions with other forcings.

It should be noted that the realistic orography is spectrally truncated in the NSWM. This truncation leads to a low spatial resolution of the orography in the model, for example, small orographic features such as those over the Arabian region and South Africa may not be well resolved (Ting 1994). The low representation of the orography therefore would potentially influence the related nonlinear individual and nonlinear effects with other forcings. As shown in Fig. 6, the existence of orography is indeed the most important source for the nonlinear interactions. Thus, the orography could induce significant changes in thermal forcing (Ting 1994) and affect the distribution of transient eddy flux (Held et al. 2002). However, we speculate that the lower spatial resolution of orography does not have a great impact on our results, at least for the generation of anomalous stationary waves over Europe. On one hand, the NSWM responses to the nonlinear interaction effects between orography and other forcings tend to be largest in the Pacific-North American sector, which is beyond the scope of this paper, in agreement with Held et al. (2002). On the other hand, Iles et al. (2019) suggest that increasing the resolution of the orography and coastlines in climate simulations have limited benefits for temperature extremes over Europe, except in reducing hot biases over mountainous regions.

Examination of the pattern correlation, the responses due to interactions involving diabatic heating and transient momentum fluxes (0.60), and orography and transient momentum fluxes (0.45) make a significant positive contribution to the simulation of the observed anomalous stationary waves. This significant pattern correlation of the nonlinear interference between diabatic heating and transient momentum fluxes proves our assumption in Sect. 3.3.2. Moreover, the contribution of the three-way interaction term is also non-trivial and surprisingly makes a significant damping contribution ($-0.63$). The streamfunction responses due to the nonlinear effect between diabatic heating and orography makes a modest contribution to the total response. However, the results of the pattern correlation should be interpreted cautiously as the magnitude of the responses is relatively low. By inspecting Fig. 5a–c and d–f, we can confirm that the nonlinear interference between diabatic heating and transient momentum fluxes and orography and transient momentum fluxes dominate more than the others (Fig. 5d–f). While the magnitude of these interaction terms are rather small, the larger magnitude of the response to the total forcing indicates that a nonlinear resonance is amplifying the anomalous stationary wave response.

Fig. 6

Same as Fig. 4, but for the NSWM response to nonlinear interactions between the three different forcing terms

3.4 Interrelationship between various nonlinear effects

To explore further the robustness of the above results, we examine the interrelationship between all the nonlinear combinations, direct, individual total, and interaction effects involved in the NSWM for all the simulations for all 38 heat wave events (Table 1). Table 2 shows the pairwise correlations for pattern correlations based on Spearman’s rank correlation. The pattern correlation is calculated for each simulation and decomposition between the response of the NSWM to every nonlinear effect and the observed anomalous stationary wave for the 38 events.

The dominant role of the transient momentum fluxes is corroborated. Obviously, the total response of the NSWM (driven by all forcings, TNE_QMO) has a significant positive correlation with most nonlinear effects of the transient momentum fluxes, for example, the direct (INE_M, 0.81) and individual total nonlinear effects (ITNE_M, 0.88), and its combinations with diabatic heating (TNE_QM, 0.92) and orography (TNE_MO, 0.87). The individual total nonlinear effect of orography (ITNE_O) has a significant positive correlation with the nonlinear interaction effect involving transient momentum fluxes and orography (TENI_MO, 0.59). The significant positive correlations also exist in the direct nonlinear effect of diabatic heating (INE_Q) with its total nonlinear effect (ITNE_Q, 0.62) and its combination with the orography (TNE_QO, 0.74). However, consistent with our previous results, other nonlinear effects including diabatic heating, which can cause negative contributions to the total response of NSWM, albeit these negative responses do not pass the significant test as shown in the second column in Table 2.

Table 2 Correlations between the various decomposed forcing components of the NSWM for simulated anomalous stationary waves of the 38 simulations

3.5 The role of transient momentum fluxes

The transient momentum fluxes consist of two components: vorticity and divergence fluxes as expressed in Eq. 4 and Eq. 5. We have shown the dominant role of transient momentum fluxes in maintaining the anomalous stationary waves. Thus, it is of importance to identify which of the components are more vital in this dynamic process. For this purpose, we conduct simulations similar to section 3.2.2. We employ the same global diabatic heating (Q) and orography (O) but replace the transient momentum fluxes (M) by the transient divergence (Fig. 7a) and transient vorticity fluxes (Fig. 7b), respectively. As shown in Fig. 7, the spatial pattern of the responses to both components is quite similar and shows an extended anti-cyclonic system over Europe. However, the response to the transient vorticity fluxes can capture the similar location of the center of the anti-cyclonic system compared with the simulated anomalous stationary waves (Fig. 4a) and the observed anomalous stationary waves (Fig. 2b). The result demonstrates the contributions from transient vorticity fluxes to the total responses in the NSWM are more important than those from divergence fluxes.

The pattern correlations support the importance of transient vorticity fluxes in maintaining the anomalous stationary wave over Europe. The pattern correlation between the response of transient vorticity fluxes with the observed anomalous stationary waves is 0.48 and is larger than the value of 0.07 of the transient divergence fluxes. This suggests that convective processes, driving divergent flows, play a rather minor role during European heat wave events.

We further decompose transient vorticity fluxes into high- (periods of less than 10 days) and low-frequency (periods between 10 and 30 days) eddies. Thus, Eq. 4 is decomposed into four nonlinear interaction components among different frequency scales:

$$\begin{aligned} TF_V= & {} -\nabla \cdot \overline{({\mathbf {v}}_L^\prime \zeta _L^\prime )} - \nabla \cdot \overline{({\mathbf {v}}_L^\prime \zeta _H^\prime )} - \nabla \cdot \overline{({\mathbf {v}}_H^\prime \zeta _L^\prime )} \nonumber \\&- \nabla \cdot \overline{({\mathbf {v}}_H^\prime \zeta _H^\prime )} \end{aligned}$$

(9)

where the subscript L denotes low- and H high-pass filtered eddies.

Broadly, the responses of the NSWM to these four components present a similar spatial pattern in terms of the magnitude and locations of anti–cyclonic and cyclonic systems (Fig. 8). The responses to the terms in Fig. 8b and Fig. 8d generate an anti-cyclonic system with weak positive NAO patterns (Li et al. 2020), but the other two terms just play a modest role. Moreover, the responses to high-frequency transient vorticity fluxes (Fig. 8d) exhibit an identical anti-cyclonic system compared with the observed anomalous stationary wave in Fig. 2b. By comparison of Fig. 7b, contributions from high-frequency transient vorticity fluxes dominate the response of the NSWM to the total transient vorticity fluxes regarding the location of the center of the anti-cyclonic system over Europe.

The above analysis is also confirmed by the calculation of the pattern correlation. The significant pattern correlation between the responses to high-frequency transient vorticity fluxes and the observed anomalous stationary wave is 0.78 over Europe, which is larger than other components and also than the total responses to transient vorticity fluxes. The other three components, representing the nonlinear interactions between low- and high-frequency transient waves, either play a damping role or make only modest contributions as shown in the pattern correlations (Fig. 8a–c). The pattern correlation of Figs. 8b and Fig. 2b has a negative value ($-0.29$) over Europe while the other two show no clear pattern correlations with the observed anomalous stationary wave.

Our results reveal that high-frequency transient vorticity fluxes dominate in the maintenance of observed anomalous stationary waves, which is consistent with Schubert et al. (2014). They found that the submonthly transient vorticity flux forcing dominates in the leading modes of boreal summer stationary Rossby waves which account for more than 60% of the surface temperature variability. Teng et al. (2013) also stated transient vorticity fluxes are instrumental during the life cycle of U.S. heat waves by contributing to the maintenance of streamfunction anomalies. Further research is needed to clarify the dynamic mechanisms of the transient vorticity fluxes in maintaining anomalous stationary waves, since this may be due to reduced low-frequency variability and persistent extreme events (Coumou et al. 2015) or by maintaining the waveguide teleconnections in summer (Teng and Branstator 2019).

Fig. 7

Same as Fig. 4a, but the transient momentum fluxes is replaced by separate transient divergence fluxes a and vorticity fluxes b

Fig. 8

Same as Fig. 7b, but now the transient vorticity flux is decomposed into different nonlinear interaction parts stemming from different frequency ranges: a $-\nabla \cdot \overline{({\mathbf {v}}_L^\prime \zeta _L^\prime )}$, b $-\nabla \cdot \overline{({\mathbf {v}}_L^\prime \zeta _H^\prime )}$, c $-\nabla \cdot \overline{({\mathbf {v}}_H^\prime \zeta _L^\prime )}$, and d $-\nabla \cdot \overline{({\mathbf {v}}_H^\prime \zeta _H^\prime )}$

3.6 The role of diabatic heating

3.6.1 Distribution of diabatic heating

Simulations with the NSWM show that global diabatic heating anomalies contribute to the maintenance of the anomalous stationary waves and heat waves mainly due to nonlinear interactions with the other two forcings. Both local and remote diabatic processes in the troposphere could contribute to the responses of the NSWM. Thus, we first investigate the spatial distributions of diabatic heating during European heat waves. Figure 9a shows the vertically averaged (900–150 hPa) diabatic heating anomalies for the composite of the 38 heat wave events (Table 1). Diabatic cooling dominates over Europe, except over the UK which experienced diabatic heating. The main reason for the diabatic cooling is the radiative cooling in the free atmosphere (Binder et al. 2017). Note that diabatic cooling is pervasive over Europe but with positive temperature anomalies (Fig. 2a). This is explained by the blocking systems in the upper-level troposphere. Zschenderlein et al. (2019) demonstrated that the descending motions of the air parcels lead to adiabatic warming which can overcompensate the diabatic cooling. Hence, the surface temperatures are anomalously high over Europe. As for the UK, diabatic warming may result from surface sensible heat fluxes, which is crucial in determining high temperatures near the surface.

Another important feature in Fig. 9a is that the diabatic heating anomalies are concentrated over the mid-latitude North Atlantic region ($30^\circ N$–$60^\circ N$) and over tropical regions ($30^\circ S$–$30^\circ N$). Heating in different regions may play different roles in influencing the large-scale circulation anomalies in nearby or remote regions (Teng and Branstator 2019). Hence, elucidating the impact of various heating regions on European heat waves is fundamental for a better understanding of the physical mechanisms and can also provide insight into improving the predictable skill of European heat waves.

Fig. 9

a Vertically averaged (900 hPa to 150 hPa) diabatic heating anomalies. b–e The NSWM response to different regional sources of mid-latitude and tropical diabatic heating. The dashed boxes represent the regions used for regionally decomposing the diabatic heating fields: Mid-latitude North Atlantic (MATL, 0–$60^\circ W$; red), East Tropical Pacific (ETPC, $70^\circ W$–$180^\circ W$; blue), Tropical Atlantic (TATL, $70^\circ W$–$20^\circ E$; black), and the Tropical Indian Ocean and West Pacific (TIWP, $60^\circ E$–$180^\circ E$; green). The contour interval is $1 \times 10^6\,\,\mathrm{{m}}^2\,\mathrm{{s}}^{-1}$

3.6.2 Diabatic heating in the mid-latitude North Atlantic

There is a strong positive diabatic heating anomaly over the mid-latitude North Atlantic (MATL) (Fig. 9a). This motivates us to perform NSWM simulations to determine the contribution of this heating source to the maintenance of European heat waves. A control simulation is first conducted, by running the NSWM with the climatological mean (60-year) of the JJA diabatic heating and keeping the transient momentum fluxes and the orography identical, similar to the simulation in Sect. 3.2.2. We then extract the anomalous diabatic heating over the considered MATL region from the global anomalous diabatic heating field which we used in the simulations in Sect. 3.2.2. This extracted diabatic heating is added to the 60-year mean diabatic heating data in the control simulations. Based on these data, the new simulation is set up and all other forcings remain the same as in the control simulation. Then the responses of the control simulation are removed from the new simulation, by doing so, we can obtain the sole contributions from the MATL heating source in maintaining the observed anomalous stationary waves over Europe.

The responses of the NSWM to the MATL heating regions are shown in Fig. 9b. A clear anti-cyclonic system is generated that extends from the mid-latitude North Atlantic to western Europe over the heat wave region (Fig. 2a). This anti-cyclonic system could partly contribute to the observed anomalous stationary waves over Europe in Fig. 2b. Our results confirm the role of heating over the MATL region in maintaining European heat waves, in agreement with Ghosh et al. (2017, 2019). However, the positive response over the MATL region is in contrast to the observed cyclonic system over this region (compare Fig. 9b with Fig. 2b). This indicates that heating over other regions makes likely substantial contributions as well. Our result is consistent with Black et al. (2004) to some extend. By using a wind trajectory model, they demonstrated that air masses trapped within the anticyclonic systems over western Europe can travel very short distances, mainly from the MATL region for example during the summer 2003 European heat wave event.

3.6.3 The various sources of tropical diabatic heating

We then focus on the contributions from different regions of tropical diabatic heating. As previous studies indicated, tropical heating can drive large-scale circulation anomalies and trigger intensified subsidence in remote regions over mid- and high-latitudes (Cassou et al. 2005; Bader and Latif 2005; Park and Lee 2019). Building upon these previous studies and our analysis in Sect. 3.1 about large-scale climate patterns and European heat waves, three tropical heating sources are considered to ascertain their contributions, as shown in Fig. 9a: the east tropical Pacific (ETPC), the tropical Atlantic (TATL), and the tropical Indian Ocean and West Pacific (TIWP).

With similar simulations as used for obtaining the response displayed in Fig. 9b, the responses of the NSWM to the three tropical heating regions are shown in Fig. 9c–e, respectively. Generally, the responses to the three sources of tropical diabatic heating all lead to a wave-like response over the mid-latitudes at the upper-level troposphere, regardless of the magnitude. One of the most prominent responses is the one where the NSWM is forced by heating over the TIWP region (Fig. 9e). The response to this forcing generates an approximately circumglobal pattern with different amplitudes and locations of anomalous stationary waves. An anti-cyclonic system is fostered with its center located over the Mediterranean region. The maintenance of this system could be responsible for high-temperature extremes in western and southern Europe, especially around the Mediterranean region. Our diagnosis of heating in the TIWP region confirms the important role of IOD events in modulating temperatures in Europe, as discussed in section 3.1. This tropical-extratropical connection is also consistent with Behera et al. (2013). They demonstrated that European heat waves are associated with a Rossby wavetrain via a circumglobal teleconnection. This wavetrain is determined by the anomalous diabatic heating over the TIWP region. This anomalous diabatic heating could trigger perturbations in the upper tropospheric Asian jet and transmit this signal to western Europe from the west through the mid-latitude wave-guide during the positive IOD events, and then making major contributions to heat waves in western Europe.

Concerning the anomalous diabatic heating over TATL, we find that an anti-cyclonic system is produced off the coast of north-western Africa (Fig. 9d). This system can reach the Iberian Peninsula and could be accountable for the high-temperature anomalies over there. This result is partly in agreement with Cassou et al. (2005), where they pointed out that tropical Atlantic anomalous heating could increase the probability for hot days to occur over Europe associated with the intraseasonal variability of the heating. They also stated that the tropical Atlantic heating acts to amplify the residence frequency of positive NAO events and the Atlantic low events, however, we do not capture these mechanisms by using the NSWM. Perhaps the main reason is that the NSWM seems to underestimate the role of transients in the positive feedbacks between the eddy mixing associated with wave breaking and both the deformation and deflection of jet stream (Held et al. 2002), which are crucial during the summer NAO life cycle (Luo et al. 2007; Franzke et al. 2004). Besides, our results demonstrate that the responses to tropical heating in the ETPC region are confined to local areas and the strong responses are mainly over the East Pacific (Fig. 9c). Although the positive responses appear over Europe, the magnitude is rather small.

To sum up, while previous studies indicated that the North Atlantic mid-latitude sea surface temperature forcing is important for European heat waves (Cassou et al. 2005; Ghosh et al. 2017, 2019), here we demonstrate that IOD and tropical heating in west Pacific are also responsible for long-lived European heat waves (38 cases at least with a threshold of 5 days). Thus, our study provides another new understanding of what leads to long-lived European heat waves.

4 Conclusions

In this study we examined the dynamical mechanisms involved in the maintenance of European heat waves from a nonlinear anomalous stationary wave perspective. For this purpose, we first identified 80 heat wave events, which lasted at least 3 days, for the period 1958–2017. Most events (76.3%) occurred after 2000, indicating the effect of global warming over Europe. We use 38 heat wave events, which lasted at least 5 days, to explore their relationships with large-scale teleconnection patterns. We find that European heat waves tend to occur during positive IOD events. This relationship is confirmed by examination of the two-dimensional probability distribution function of surface temperature anomalies and the IOD index.

Composite analysis of the 38 heat wave events shows a clear observed anomalous stationary wave located over Europe. We use the NSWM to simulate the observed anomalous stationary waves to better understand the maintenance mechanisms involved. While there are some discrepancies, the main features of the observed anomalous stationary waves over Europe are reproduced well when forced by global anomalous diabatic heating, transient momentum fluxes, and orography.

We then use the factor separation method to decompose the total responses from the NSWM. This method also enables us to further examine the nonlinear interactions in maintaining anomalous stationary waves and heat waves. Based on a visual inspection and pattern correlation, it is found that transient momentum fluxes play a dominant role in maintaining the observed anomalous stationary waves. But the NSWM responses to transient momentum fluxes are alleviated by the nonlinear interaction between the three forcings. Further decomposing the transient momentum fluxes reveals that transient vorticity fluxes are more important than transient divergence fluxes. We also find that transient momentum fluxes are dominated by the high-frequency rotational component, with the divergent component being much smaller.

Different from the maintenance of summer climatological stationary waves (Ting 1994), our results indicate that diabatic heating tends to make negative contributions to maintaining the observed anomalous stationary waves during heat waves, especially over the UK and eastern Europe. Furthermore, through nonlinear interaction with other forcings, orographic effects make moderate contributions to maintaining the observed anomalous stationary wave, especially via the nonlinear interaction with diabatic heating.

Our results show that the nonlinear interactions between the three forcing terms make substantial contributions. In particular, while diabatic heating as a individual forcing leads to rather small anomalous stationary wave responses, together with the transient fluxes it makes a substantial contribution.This indicates that all three forcings are necessary for maintaining anomalous stationary waves and European heat waves.

We further investigate the distribution of diabatic heating and find that strong diabatic cooling dominates over Europe during heat wave periods. This diabatic cooling indicates strong descending motions due to a persistent anti-cyclonic system. This system can lead to a large magnitude of adiabatic warming overcompensating the diabatic cooling and increasing the surface temperature. Further NSWM simulations show that tropical heating located in the tropical Indian Ocean and Western Pacific (TIWP) and tropical Atlantic (TATL) forces the observed anomalous stationary waves remotely, and the former plays a stronger role than the latter. Moreover, heating over mid-latitude North Atlantic also makes contributions to the anomalous stationary waves mainly over southwestern Europe.

Overall, our results in this study offer a new dynamical perspective on understanding the maintenance of anomalous stationary waves associate with heat waves. The successful application of the NSWM enabled us to examine the relative contributions due to various nonlinear effects in maintaining heat waves. From a forecasting perspective our research findings show the potential of improving seasonal forecasts or which predictors to select for statistical forecasting of European heat waves.