Influence of the spatial distribution of gravity wave activity on the middle atmospheric circulation and transport

Analyzing GPS radio occultation density profiles, we have recently pointed out a localized area of enhanced gravity wave (GW) activity and breaking in the lower stratosphere of the Eastern Asia/North-western Pacific (EA/NP) region. With a mechanistic model of the middle and upper atmosphere experiments are performed to study a possible effect of such a localized GW breaking 15 region on the large-scale circulation and transport and, more generally, a possible influence of the spatial distribution of gravity wave activity on the middle atmospheric circulation and transport. The results indicate an important role of the spatial distribution of GW activity for the polar vortex stability, formation of planetary waves (PW) and for the strength and structure of the zonal mean residual circulation. Also, a possible effect of a zonally asymmetric GW breaking in the longitudinal 20 variability of the Brewer-Dobson circulation is analyzed. Finally, consequences of our results for a variety of research topics (sudden stratospheric warmings, atmospheric blocking, teleconnections and a compensation mechanism between resolved and unresolved drag) are discussed.

derived a first dynamically consistent two-dimensional picture of the mean-transport streamlines for the middle atmosphere that is often used as a basic BDC concept.However, Demirhan Bari et al. (2013) found a 3D structure of the circulation in the middle atmosphere to be in good correspondence with tracer fields, especially in relation to the zonal wave-one pattern observed in the stratosphere and mesosphere, although their study did not give a comprehensive dynamical explanation of the 5 discovered circulation structures (enhanced downward branch of BDC over north-eastern Asia, waveone pattern).
PWs are usually thought to be created in the troposphere and then vertically propagating into the middle atmosphere.The theoretical possibility of PW creation by zonally asymmetric IGW breaking was first numerically analyzed by Holton (1984) and later on, e.g., by Smith (2003) and Oberheide et 10 al. (2006), and experimentally verified by Lieberman et al. (2013).There is a building agreement in the literature on the role of wave activity in preconditioning sudden stratospheric warming (SSW; e.g.Ayarzaguena et al., 2011) events.
SSWs belong to the most pronounced atmospheric phenomena, as they cause abrupt changes of middle atmospheric circulation and tracer distribution, and they also affect tropospheric weather patterns (e.g.further statistically confirmed by Cohen et al. (2013), who interpreted it as a response of the resolved waves to maintain a ''sensible'' stable circulation.Such a response is expected, since all processes in the atmosphere are driven by the tendency to reach an energetically more favorable, stable state.In addition to the stability constraint, Cohen et al. (2014) proposed two additional mechanisms using a potential vorticity (PV) concept, PV mixing and refractive index interaction.

5
In this study, we focus on the physical mechanism and structure of the atmospheric response to the zonally asymmetric forcing represented by an artificially injected GWD in the stratosphere.We are following Šácha et al. (2015), who described a localized area of enhanced GW activity and breaking in the lower stratosphere over the Eastern Asia/North-western Pacific (EA/NP) region and discussed possible implications of this GW hotspot for large scale dynamics and transport.By artificially 10 enhancing the GWD in a 3D mechanistic circulation model of the middle atmosphere we examine the hypothesis that such a robust breaking region plays a role in forcing longitudinal variability of the BDC and can generate PWs.Further, we investigate possible implications for the polar vortex stability and the role of the GWD distribution and of the artificial forcing components (direction of the force).
The structure of the paper is as follows: in section 2 we describe the model and sensitivity simulation 15 set up together with the observational motivation and justification for an artificial GWD enhancement.
The section closes with a brief description of tracer data used in this study.Section 3 starts with an illustration of the geopotential response to different GWD injections with particular focus on effects in the polar region.We also present the dynamical impact, structure and modes of the PWs generated by the artificial GWD.Finally, we show the differences of the BDC due to the geometry of the GWD 20 modulation and analyze the 3D residual circulation spatial patterns in relation to the GWD distribution.
In Section 4, we give a summary of our results, discuss potential implications of our findings and outline future directions of our work.

25
We use the Middle and Upper Atmosphere Model (MUAM), which is a nonlinear 3D mechanistic global circulation model.It has a horizontal resolution of 5°x5.625° and extends in 56 vertical layers up to an altitude of about 160 km in log-pressure height (Pogoreltsev et al., 2007).At 1000 hPa, the lower boundary of the model, we prescribe stationary PWs of wave numbers 1, 2 and 3 obtained from decadal monthly mean ERA Interim (ERAI) temperature and geopotential reanalysis data (ECMWF, 2016).Up 30 to an altitude of 30 km the model zonal mean temperature is nudged to ERAI zonal mean temperature.This is necessary because MUAM does not account for, e.g., orography or some radiative processes including 3D water vapor or surface albedo.However, the troposphere is necessary for stationary PW forcing and the generation and propagation of traveling PWs and tides and therefore it cannot be neglected.The assimilation of stationary PWs and zonal mean temperatures is not only active during 35 the spin-up of 330 model days but also during the 30-day analysis period.
The effect of nudging during the analysis period is dependent on the strength of the artificial forcing.In the reference simulation the nudging effect is lower than 1 K/day everywhere and for the simulation with strongest forcing it reaches locally to magnitudes around 2 K/day in a zonal mean (as shown in Fig. 1c).Because in MUAM simulations only the zonal mean temperatures are nudged to the zonal mean, nudging has no direct effect on the wave structure of the response to the forcing, but is likely to reduce the magnitude of the zonal mean response.
The time step of the model is 225 s following a Matsuno (1966) integration scheme.For simulations, the model starts with a globally uniform temperature profile and no wind.During a spin-up period, the mean circulation is built, and PWs and tides are generated.After that, a time interval of 30 model days with a temporal resolution of 2h is analyzed.Since the lower boundary conditions are taken as a decadal mean January mean, this interval refers to an average January climatological state.Monthly zonal means of wind, temperature and GWD are given in Fig. 1.Owing to the constant forcing with time in the lower atmosphere, the standard deviation of temperature within these 30 days is smaller than 3K near the stratopause and mesopause and smaller than 1K elsewhere.The standard deviation of the zonal wind is largest within the jets reaching 4 m/s in the summer easterlies.These values do not have a meteorological meaning and are provided here to demonstrate that MUAM has a rather small variability within the analysis interval.

15
GWs are parameterized after a linear Lindzen (1981) type scheme updated as described in Fröhlich et al. (2003) and Jacobi et al. (2006), and they are initialized at an altitude of 10 km with six different phase speeds ranging from 5 to 30 m/s, each propagating in eight different azimuth angles, and with GW vertical velocity amplitudes with an average value of 0.01 ms −1 .As input for the GW parameterization scheme, we modified the GW source function to reflect a distribution based on the 20 mean January field of the potential energy of disturbances computed from FORMOSAT3/COSMIC radio occultation density profiles between the tropopause and 35 km altitude taken from Šácha et al. (2015).The GW weights are calculated from these data by dividing the potential energy at each grid point by its global mean.This setup has a positive impact on some climatological features in MUAM.
Nevertheless, the effect on the horizontal distribution of GWD in the stratosphere is negligible.We will 25 refer to this setup as the reference simulation.Zonal (gcu) and meridional (gcv) flow acceleration as well as the heating due to breaking or dissipation of GWs (gt) is calculated by the parameterization scheme.
To examine and to demonstrate the effect of spatial distribution of the GW activity we performed a set of sensitivity simulations (Table 1) with artificially changed GWD imposed on the model by 30 modulating the GW parameterization output.Note that this change of GWD is only added after the spin-up so that only the 30 model days incorporate GWD changes.Thus, the simulation period also includes the temporally delayed response for the adaption from reference conditions to enhanced GWD (gcu/gcv/gt) values.The naming convention (Table 1) is given by "Gcu+distribution+gcv", where the basic value of gcu of 0.5 m/s/day is not stated .

35
The enhancement is performed for a certain 3-dimensional box in the lower stratosphere (about 18-30 km) above the EA/NP region (37.5-62.5°N/112.5-168.8°E),according to the area of enhanced GW activity described by Šácha et al. (2015).This refers to the "box" distribution in Table 1 (an example is shown in Fig. 2, left panel).There are no exceptional GWD values in the reference simulation in this region.In a second version we additionally averaged the respective GWD parameters zonally within the same latitude range like the box.This way, we obtain a zonally uniform distribution, i.e. a ring of enhanced GWD parameters instead of a box but with a smaller local magnitude.We refer to this configuration as ring or "Zon" simulations (see Table 1).For all simulations, the GWD parameters outside the box or the ring, respectively, remain unchanged and are not influenced by the enhancement.

5
We are not smoothing the boundaries of the artificial enhancement area and the step between artificial and background GWD values is dependent on the horizontal location, the time step and, most importantly, the altitude level.To illustrate the sudden and localized effect of GW breaking, we have chosen to enhance the GWD in our simulations stepwise and rather abruptly.As suggested by Cohen et al. (2013), such a sharp change (as at the boundaries of our enhancement) leading to dynamic 10 instabilities is likely to induce compensation processes.
Although it is impossible to directly compute the GW drag force from current satellite measurements alone (Alexander and Sato, 2015), Ern et al. (2011) gave a methodology to estimate absolute values of a "potential acceleration" caused by GWs (maximum zonal mean values of 3 m/s/day below 40km).
Using ray tracing simulations Kalisch et al. (2014) estimated a zonal averaged GWD to be around 20 15 m/s/day in the lower stratosphere.In our model simulations we are injecting three values of additional artificial zonal component of GWD, -0.5 m/s/day as a conservative enhancement and -10 m/s/day to demonstrate a big impact of the injection.In addition, an extreme case with -70 m/s/day is added to force substantial circulation changes.
Depending on the GW type and on the direction of background winds the GWD has also a meridional 20 component, which is usually poorly constrained by observations.We performed simulations with three different values of meridional GW induced acceleration (-0.5 m/s/day, -0.1 m/s/day, 0.1 m/s/day).
Directions of the zonal and meridional GW induced acceleration were chosen based on the assumption that the majority of GWs in the EA/NP region in January will be of orographic origin and taking into account the prevailing directions of horizontal winds in the EA/NP region (see Šácha et al., 2015).On

25
this basis we argue that the 5:1 ratio between the zonal and meridional GW induced acceleration is the most realistic and therefore we choose the Box0.1 (and Zon0.1) simulation as a representative conservative enhancement for most of the analyses in this paper.A comprehensive discussion of our sensitivity simulation set-ups is given in the Discussion section.

Residual circulation 30
To highlight the importance of the stratospheric research in the EA/NP region and to show the robustness of our claim of an enhanced branch of the BDC in this region we present in the Supplement the 1978 to 2008 average total ozone January mean distribution from the ozone Multi Sensor Reanalysis version 1 (MSR1; van der A et al., 2010van der A et al., ) data (TEMIS, 2016)).Additionally, in the Supplement, the comparison is shown between the vertical structure and longitudinal variability of the 35 residual circulation and zonal cross sections of Michelson Interferometer for Passive Atmospheric Sounding (MIPAS) CH 4 volume mixing ratio profiles (KIT, 2016;see von Clarmann et al., 2009;Plieninger et al., 2015).However, the interpretation of the differences of the distributions must be done with care, since the tracer distributions result from several different processes in the atmosphere, namely advective transport, mixing, and chemical reactions (Garny et al., 2014).Also, the residual velocities are closely related to Lagrangian-mean velocities up to O(α 2 ) only for small amplitude (α ) steady waves (Bühler, 2014).
In the section 3.3, we study consequences of the GW hotspot for the longitudinal variability of the residual circulation (and BDC consequently) by means of the time mean 3D residual circulation according to Kinoshita and Sato (2013).The time averaging inserts additional uncertainty in the 3D residual circulation concept.Unlike Demirhan Bari et al. (2013), who based their analysis on monthly means and daily eddies, we are employing a 5-day running average on the 6 hourly MUAM output fields.This configuration gives the strongest zonally averaged Stokes drift from several choices of the running mean.But, it is still smaller than the value of the Stokes drift resulting from transformed 10 Eulerian mean equations, which is computed in this study according to Hardiman et al. (2010) for logpressure height vertical coordinate models.

Results
To establish the timescales of the response, in Fig. 3 we show Hovmöller diagrams of the zonal mean zonal wind and its variance.The time evolution is presented at the 6.25 hPa level (around 35.5km logpressure height, 13th model level).This level was chosen for our analyses because it is above the location of the artificially modified area and above the nudging extent, so it contains the atmospheric response only.In Fig. 3a, a Hovmöller diagram is given for the zonal mean zonal wind at 6.25 hPa level in the Ref simulation documenting that the model is essentially steady.Fig. 3b shows the time evolution of variance of the zonal mean wind anomaly (Box0.1-Ref)and Fig. 3c shows the time 20 evolution of zonal mean zonal wind for the 10box simulation.We can see the buildup of the response until approx.day 7 after the GWD injection and after that the structure of the response remains quasisteady, with small variations of the magnitude only.
In contrast to this, the zonal mean zonal wind time evolution from the so-called SSW simulations (Fig. 4a and 4b) do not reach a steady state in the course of the 30 days simulation and therefore the results

25
based on those simulations are presented at particular time steps or as animations in the Supplement.
Results of other simulations (Table 1) are averaged across the quasi-steady state (7th-30th day of the simulation) and are supplemented with the estimate of statistical certainty or standard deviation of the mean.
Except for the SSW simulations, our study is focused mainly on the mean response to a monthly mean 30 GWD distribution, because from observational analyses we usually have information on the IGW activity distribution on a monthly or seasonal basis (Šácha et al., 2015).The short-term response, which would be arguably more relevant to the real atmosphere taking into account the intermittency of large amplitude GWs (e.g.Hertzog et al., 2012;Wright et al., 2013), is not well captured by the mechanism of constant GWD injection, which will be discussed in the final section.Still, there are some interesting  3D)) shows the same as the 2 nd row, but for the artificial ring GWD configuration.Note the different scaling of the color bars, which is chosen according to the maximal and minimal value of geopotential (anomaly), so that the labels of the color bar give direct information 10 on the magnitude of the differences in geopotential response.
The anomalies and differences are analyzed with special focus on the polar vortex response, since it will be shown below that the dynamical response on GWD changes is strongest in the polar region.This comparison demonstrates not only the importance of the role of the longitudinal distribution of the zonal mean drag force but it also highlights an important and different effect of each of the individual The respective geopotential responses in the Box0.1gcv and the Box0.1gtsimulations have almost exactly opposite features, as the positive gt enhances geopotential in the upwind and northward 25 direction from the GWD region, while artificial northward deceleration has an opposite effect.
Although we used a nonlinear model, the additivity of effects of different GWD components (Figs.

5(2B))
, where the forcing constitutes of exactly these components.Also, the differences between simulations with different meridional drag (compare Figs.Another two important results are visible from the comparison of the plots in the second and third row of Figure 5. Firstly, there are much bigger anomalies for the box enhancements (second row) than for 35 the corresponding ring enhancements (third row).This is true locally as well as in the zonal mean (compare the sum of geopotential response given in the legend for Figs.5(2A), 5(2B), 5(3A) and 5(3B)).In the box simulations (Figs.5(2A) and 5(2B)) the response is typically dominated by a rather meridionaly oriented dipole pattern with a localized positive geopotential anomaly at the center of the polar vortex and negative geopotential anomaly at the location of Aleutian High.In the correspondent ring simulations (Figs.5(3A) and 5(3B)) the geopotential response is more zonally uniform.
Secondly, there are large and significant differences (50% or 25% of the magnitude of the anomaly) between box simulations with slightly different setups of the meridional drag (Fig. 5(2B) vs. Figs.
As noted above, the magnitude of the geopotential response is larger for the box enhancements than for the ring enhancements.For the Box0.1 simulation, the geopotential anomaly at the 6.25 hPa level 10 reaches about 20 gpm in a monthly mean.The horizontal wind anomaly for the Box0.1 simulation (Fig.

5(2B)
) reaches maximal values slightly below 1m/s.Anomalies for the 10box simulation (Fig. 5(2A), 20 times bigger eastward deceleration than for Box0.1) are almost exactly 20 times stronger and show a very similar dipole pattern.Although locally the difference between these two simulations may seem to be linear, this comparison is misleading, since both simulations (10box and Box0.1) have different ratio between the strength of GWD components.This means, for example, that the drag force has different orientation between these two simulations.
Unexpectedly, the box simulations lead to anomalies that would contribute to weakening rather than amplification of the Aleutian High.Based on the results and discussion of Šácha et al. (2015), who argued that the EA/NP hotspot (high GW activity already in October/November) may play a role in the 20 onset of the winter circulation in the stratosphere in this region, we expected a positive contribution of the GWD response to the background climatology (e.g.contribution to the unusually hot temperatures in the stratosphere in the EA/NP region by induced subsidence).
In Figs.4a and 4b we presented a time evolution of the zonal mean zonal wind at 6.25 hPa for the SSWbox and SSWzon simulations with signs of a wind reversal at polar latitudes at particular time 25 steps suggesting an occurrence of a minor SSW.Now we show additional results from the SSWbox and SSWzon simulation in the form of animations of geopotential and horizontal wind field response at 6.25 hPa (Animation 1a and 1b in the Supplement) and a snapshot at 280 hours after the injection to show the response at a developed stage of the SSW.In response to a strong GWD increase in a box we observed a vortex displacement (Fig. 6a and Animation 1a in the Supplement) and in response to a strong GWD increase in a ring we obtain a vortex split like event (Fig. 6b and Animation 1b in the Supplement).
In the SSWbox simulation (Animation 1a), immediately after the spin-up period when the GWD starts to be artificially modified (injection of GWD), a geopotential ridge begins to form above the Northern Pacific (northward from the GWD area).This anomaly strengthens and shifts a little westward above Siberia, where, from approximately 5 days of the GWD injection, we observe an evolution of a pressure high.All the time the vortex is shifting towards the northern boundary of Northern America where it stays till the end of the simulation.
In the SSWzon simulation we observe a slow (compared to the SSWbox simulation) creation of a pressure high above the Northern Pacific together with a high-pressure ridge above the Northern Atlantic.This pressure high is almost stationary (in contrast to the SSWbox) leading to the vortex split approximately ten days after the injection.This is a potentially very interesting result suggesting that a symmetric forcing favors vortex split and localized forcing favors displacement events, but the 5 robustness of this claim needs to be tested in future work for various initial vortex states.
For illustration, in Figure 6 we show the geopotential field and horizontal wind speed 280 hours after the GWD injection, when the vortex split develops (Fig. 6b) and the vortex displacement is in its mature state (Fig. 6a).The vortex displacement event develops more quickly, as seen from comparison of Fig. 4a and 4b or from the animations 1a and 1b in the Supplement.However, both events have 10 limited vertical extent, and do not disturb the entire vortex (only up to 60 km of log-pressure height; not shown).

Creation of planetary waves and dynamical impact
In this section we compare PW activity and amplitude structure of the leading PW modes between reference, box and ring simulations.We show results of the E-P flux diagnostics and Fourier transform 15 (FT) analysis of geopotential anomalies.Fig. 7 shows the mean (7th-30th day) E-P flux and its divergence for the Ref simulation (Fig. 7a), Box0.1 simulation anomalies (Box0.1-Ref;Fig. 7b), Zon0.1 simulation anomalies (Fig. 7c), the difference between the Box0.1 and Zon0.1 simulations (Fig. 7d) and mean mean E-P flux and its divergence for the 10box simulation (Fig. 7e) and respective anomalies (Fig. 7f).Note that we show 20 the E-P flux divergence as a force per unit area (units [ !!  !! ]), not as an induced acceleration (units [ !! ]), as in Hardiman et al. (2010), because otherwise upper stratospheric and mesospheric effects would dominate the plots due to the density decrease with height.The statistical significance of the mean E-P flux divergence differences has been computed by a t-test and regions with p values < 0.05 are stippled.

25
In Fig. 7b, for the Box0.1 and Ref simulation differences, we find an anomalously weak E-P flux convergence (positive difference to the Ref simulation) centered at the equatorward flank of the GWD enhancement area and an anomalous convergence in a broad area around 60°N.This pattern is similar for the Zon0.1 simulation anomalies (Fig. 7c), but much weaker and with the anomalous convergence starting more poleward.It is also similar in the 10box simulation anomalies (Fig. 7f), but much 30 stronger in magnitude (approx.20 times).In all of those simulations, this anomalous pattern is limited in altitude and only slightly exceeds the vertical boundaries of the GWD area (especially in the polar region).
Taking into account the reference E-P flux field (Fig. 7a), the anomalies can be caused by two different mechanisms.The first one is an indirect mechanism, when the artificial GWD drag modifies the winds 35 causing changes (with respect to the Ref simulation) in propagation conditions for PWs propagating from below (for more details on the refractive index interaction see Cohen et al. (2014).According to this mechanism, the E-P flux and its divergence anomalies and differences (Figs.7b, c, d, f) would be associated with a stronger poleward and weaker/stronger upward propagation of PWs in the 10box/Box0.1 simulation along the northern edge the polar night jet in the northern part and northern boundary of the GWD area.The E-P flux divergence anomaly at the southern flank of the GWD would be associated with a suppression of upward and equatorward PW propagation elsewhere in the GWD region.But this mechanism fails to explain some features in Fig. 7, e.g. the E-P flux divergence emerging in the E-P flux field in the 10box simulation (Fig. 7e).Therefore, although the changes in the refractive index will definitely be present in the artificial GWD simulations, we have an indication that another mechanism is dominant.
This second mechanism is associated with the evidence given by Holton (1984) that a zonally asymmetric GW breaking possibly generates PWs in the mesosphere.In connection with that the 10 artificial GWD enhancement in a box would cause displacements of fluid particles (in the initially balanced predominantly zonal flow) and thus generate a broad spectrum of waves of different types depending on background conditions and on the geometry of the drag region.We can find support for this mechanism from the E-P flux difference between Box0.1 and Zon0.1 simulation (Fig. 7d).
In the previous section, we have shown that the box enhancement induces a stronger zonal mean 15 geopotential response than the corresponding ring enhancement.So we can assume that the first mechanism has bigger effect in the box simulations, which is true for the E-P flux divergence difference (Fig. 7d).However, regarding E-P flux vectors, Fig. 7d reveals that there are not only differences in magnitude between Box0.1 and Zon0.1 E-P flux anomalies, but also that the Zon0.1 simulation lacks the horizontal component of the anomalous E-P flux, with biggest differences in the 20 latitudinal band encompassing the artificial GWD area.This latitudinal band is not significant in Figs.7b,c and d, because the plotted t-test results are based on the difference of the E-P flux divergence (not on the magnitude of the E-P flux vector difference).From Figs. 7 b, d, f, we see that the anomalous PWs are generated at the southern flank of the GWD area and propagate predominantly northward (with a small downward component), where they cause anomalous convergence between 60°-80°N.

25
For the conservative Box0.1 simulation, the anomalies in the E-P flux divergence are about 5% of the reference values.Zon0.1 E-P divergence anomalies (Fig. 7c) reach only 1-2% of the reference values, locally.Anomalies of the 10box simulation (Fig. 7f) exhibit the same pattern as Box0.1 anomalies, but the magnitude is much stronger -more than 50% of the reference E-P flux divergence values.
Therefore, we observe an influence of the 10box GWD enhancement also in the mean field in Fig. 7e, where the artificial GWD box demonstrates itself as an E-P flux divergence area on the southern flank of the GWD enhancement region.This is another supporting argument that the box enhancement generates PWs, with further evidence given below.
In Fig. 8, E-P flux diagnostics is presented at particular moments (1 and 5 days) after the GWD injection for the SSWbox and SSWzon simulations.The anomalous E-P fluxes in those highly 35 nonlinear simulations absolutely overcome the reference fields, so that we can directly observe the generation and propagation of PWs generated by the artificial GWD.However, for these simulations the structure of the E-P flux divergence area changes with time and also the propagation directions of PWs created in this region are time dependent.So, we have chosen to present snapshots from the 1 st and 5 th day to demonstrate particular features of the box GWD enhancement.For interested readers, the full time evolution is given as Animation 2 in the Supplement.
In Figs.8a, b, one can clearly see the generation of PWs by the box enhancement.Five days after the GWD enhancement (Fig. 8b), the E-P flux divergence region extends almost over the whole GWD area.Anomalous PWs propagate equatorward, poleward and upward with two major E-P flux 5 convergence regions around 30° N and between 60° and 80° N. One day after the GWD injection (Fig. 8a), the E-P flux divergence area is located at the southern flank of the GWD and generates horizontally, southward propagating PWs only.In Fig. 8a, in the majority of the GWD region, we can see also the first mechanism (refractive index interaction) being active, as the GWD region influences propagation of PWs propagating from below.This is the most dominant effect of the ring enhancement 10 (Fig. 8c, d), where in the SSWzon simulation we can hardly observe any anomalous PW generation and the dominant effect of this ring enhancement is altering the propagation conditions for the upward propagating PWs from the troposphere.There is a weaker propagation through the GWD region, with deflection of PWs northward and southward at the southern GWD flank.
Further indication of the creation of PWs by the GWD region is provided by the FT analysis of  f.To quantify the dispersion of the monthly mean differences, the dotted lines show the standard deviations.
The amplitude anomaly is positive starting at the northern flank of the artificial GWD (37.5-62.5°N)and further poleward.The maximum is gained between 70-75°N.Another smaller, but still significant, 25 region of positive wave-1 amplitude anomaly is located around 30°N south of the GWD.Smaller negative wave-1 amplitude anomaly lies inside the GWD area.In the Box0.1 simulation, wave-2 (Fig. 9c) has a pronounced negative amplitude anomaly inside the latitudinal belt encompassing the enhancement region.For wave-3 (Fig. 9d), we find positive anomalous amplitudes starting from central latitudes of the GWD region and ending around 80°N, although inside the GWD region the positive 30 anomaly is locally not significant.There is a negative wave-3 amplitude anomaly starting at the southern flank of the GWD region with the end around 10°N.The effect on wave-4 amplitudes is almost negligible (Fig. 9d).The ring enhancement in the Zon0.1 simulation has a negligible effect on amplitudes of harmonics, as is visible from the similarity of the Box0.1 anomalies (Figs.9c, d) and differences with Zon0.1 simulation (Figs.9e,f ).These results suggest that the box GWD enhancement 35 generates preferentially wave-1 and -3 modes in comparison to the reference and also the ring GWD configuration.
Another indication that the PWs are indeed generated by the GWD box enhancement is given in Fig. 10, where the time evolutions of the anomalous wave-1 and wave-3 amplitudes are presented.
Especially in the first approximately seven days from the GWD injection, we can observe a slow propagation of anomalous wave-1 (Fig. 10a) and wave-3 (Fig. 10b) amplitudes from the GWD region to the north.For wave-3 this propagation is visible later than for wave-1 (from approx.day 3).The oscillating patterns in Fig. 10 most likely originate from a non-linear interaction between anomalously generated inertia GWs and solar tides (see e.g.Waltersheid, 1981).Those inertia GWs are responsible for propagation of the anomalous wave activity through the Rossby wave critical layer in the tropics, across the equator, and into the Southern Hemisphere (Fig. 11).

Residual circulation response
The first row of Fig. 11 shows the mean (over the whole 30days) residual circulation mass fluxes for the reference simulation and the snapshot at 5days from the GWD injection for SSWbox simulation on 10 the right.Mean (7 th -30 th day) anomalies and differences with the respective ring configuration are given in the second and third row for the Box0.1simulation on the left and 10box simulation on the right.
There are some remarkable results visible.Firstly, even for a conservative drag enhancement (Box0.1 simulation) there are significant (dashed) differences in the magnitude of the residual mass flux between box and ring GWD distribution of up to 3% in the lower stratosphere (Fig. 11e).For the 10box simulation the differences reach about 40% and create a similar pattern as for the conservative enhancement (Fig. 11f).The largest differences between the two artificial GWD configurations are found poleward from the GWD enhancement region in the altitude range between 20 and 30 km corresponding approximately to the vertical extent of the area and are associated with a stronger subsidence north of the enhancement region in the box simulations.

20
There is a smaller region of significant differences at the southern flank of the enhancement region associated with lesser downwelling in the box simulations.These two regions of significant differences constitute together a butterfly like pattern in the box-ring differences centered at approximately 45°N (the center of the enhancement region) and influencing a shallow BDC branch.Taking into account the reference field (see Fig. 11a) we can explain this pattern by a faster northward advection starting at approx.45°N and stronger subsidence northward of 60°N.On the other hand there is less upwelling in the equatorial region (not significant for the Box0.1 simulation) and slower advection from the tropics.
The continuity is satisfied through smaller downwelling south of 60°N.
We observe a similar but stronger pattern in the anomalies (Figs.11c, d), with the mean residual circulation mass flux anomaly reaching up to 5% for the Box0.1 simulation and more than 60% for the 10Box simulation.The position of the anomalous residual circulation patterns corresponds with the E--P flux divergence anomalies (Fig. 7b, f), where, for the box simulations, we observed anomalous E-P flux divergence at the southern flank and convergence north of the GWD region.In Figs.11c and 11d the butterfly like pattern is centered more southward (35°N) than in box-ring differences and the anomalous pattern on the south of the GWD region is not as well pronounced and 35 appears to be shifted above the GWD region for the 10box simulation (Fig. 11d).
In the upper stratosphere there are anomalies up to 2% only for the Box0.1 (Fig. 11c) and locally around 25% for the 10box simulation (Fig. 11d).The box simulations (not significant for Box0.1)show weaker subsidence towards the polar vortex center than the Ref simulation in the upper stratosphere and there is also anomalously low mass flux poleward and downward between 30 and 40 km of height above the GWD enhancement region.For both box enhancements, there is a large area of statistical significant anomalies giving a weak hint of lesser upwelling in the SH stratosphere (Figs.11c, d).The differences between the two sets of box and ring GWD configuration are not significant in the SH
The fact that the mean response of the upper BDC branch is rather weak and for the most part not significant can be explained by the effect of the artificial GWD region acting like an obstacle for northward flowing wind.The GWD enhancement region (Fig. 11b, snapshot for a SSWbox run) is constantly flown around inducing a significant mean anomaly (Fig. 11c, d) with anomalous upwelling 10 in its southern part and downwelling on the northern flank.But, the GWD region (obstacle) creates also a lee wave like pattern with oscillating anomalies in the upper stratosphere and in the SH.Considering a time mean, these anomalies are small and not significant, but, at particular time steps, the magnitude of the anomalies is comparable regardless of the BDC branch.Supporting information is given in Animation 3 in the Supplement, which presents the time evolution of the zonal mean residual 15 circulation associated mass flux for the 10box simulation (on the left) together with its anomaly (on the right).One can see here the global nature of the response and gain insight into how quickly the residual circulation gets affected by the NH anomalous forcing.After few time-steps, the response is constituted by a constant anomaly corresponding roughly to an accelerated shallow BDC branch sloping down from approx.30km at the North Pole to the lowest analyzed levels at the equator.Except for this 20 region, the entire domain is dominated by anomalies seemingly descending downward from the mesosphere associated with the obstacle analogy.
The zonal structure of the induced flow, and possible consequences of the GW hotspot for the longitudinal variability of the BDC were studied by means of 3D residual circulation analysis according to Kinoshita and Sato (2013).5-day running averaging was performed.Šácha et al. (2015) 25 pointed out unusually high temperatures in the EA/NP region at 30 hPa in winter and concluded that there could be an enhanced downwelling above the EA/NP region penetrating to lower levels than elsewhere.This is in agreement with Fig. 3 in Demirhan Bari et al. (2013).Supporting results highlighting the importance of future research in this region are given in the supplement.In Fig. 1S in the supplement we present a thirty-year average January MSR total ozone column field with a total 30 ozone column maximum located in the EA/NP region.In Fig. 2S in the supplement, longitudinal crosssections of MIPAS CH 4 volume mixing ratios show a peak of subsidence around 15 km in the EA/NP region (at 140°E) and the interesting massive upwelling branch east of it.
To evaluate the possible role of the GW activity in the longitudinal variability of the BDC, we present longitudinal cross-sections of the reference 3D vertical residual velocity and Box0.1 anomalies going

35
from the northern to southern part of the artificial GWD (Fig. 12).From longitudinal cross-sections of the reference vertical residual velocity (left side of Fig. 12), we see that MUAM vertical residual velocity field is dominated by a wave-2 pattern, with the maximum subsidence branch penetrating to the lower stratosphere in the EA/NP region and with an abrupt switch to upwelling on the east.Ridges and troughs of the wave show a characteristic westward tilt with height.Šácha et al. (2015) hypothesized that the collocation of the EA/NP GW hotspot and the enhanced BDC branch can be partly a consequence of the circulation induced by the GW breaking.But the results are rather contradictory.In agreement with the zonal mean residual circulation analysis, we can see that in the southern part of the area (Fig. 12f), the GWD induces predominantly anomalous upward flow.

5
Anomalous subsidence strengthens when going further northward (Fig. 12b, d).In line with the obstacle analogy, we observe subsidence in the eastern part of the GWD region only, while anomalous upward flow dominates the western part of the GWD region, and then again eastward and slightly above the anomalous subsidence area.Similar structure of an Eulerian mean vertical velocity field has been found by Shaw and Boos (2012) as a response to an artificial torque placed in the troposphere 10 around 30°N.These results show that GWs can contribute to longitudinal variations in the BDC and not only the downwelling, but also upwelling patterns may be related with GWs.
The magnitude of the vertical residual velocity anomalies maximizes around 2% of the reference value for the Box0.1 simulation (Fig. 12b, d, f).For the 10box simulation (Fig. 3S in the Supplement) the distribution of upwelling and subsidence is identical and the magnitude reaches 30% locally.

15
Physically, such an anomalous pattern can be explained by considering the dominant background horizontal north-eastward wind together with the previously mentioned small obstacle analogy, with induced upward flow upwind and downward flow downwind from the GWD box.However, for the SSWbox simulation we can observe a completely different distribution variable with time, with subsidence dominating directly above the GWD area in the later stages of the simulation (Animation 4 20 in the Supplement).When the artificial GWD is strong enough to induce significant dynamical changes (SSW simulations) the anomalies cannot be directly explained as being GW induced because also the dynamical state of the atmosphere changes (e.g. the anticyclonic evolution in Animation 1a).
Therefore, the explanation of residual vertical wind cross-section patterns for both SSW simulations is much more complicated and requires future research allowing at least the GWD enhancement to reflect 25 the changing background conditions.

Discussion and conclusions
We will begin this section by giving a brief summary of results.Then, we will discuss limits of our results stemming from the construction of the sensitivity simulations and afterwards, we will give some conclusions for different research topics in the middle atmosphere.

30
In this paper, we presented results of a set of sensitivity simulation to find out the possible role of a localized GW hotspot and also, generally, to demonstrate the influence of spatial distribution of GWD on the middle atmospheric dynamics.The focus was on a mean response to a steady GWD perturbation injected into climatological January condition.Except for the strongest GWD enhancement (SSW simulations; Fig. 3), all simulations (Table 1) have reached a quasi steady state approximately 7 days after the GWD enhancement (Fig. 3).The average across this state was considered as a mean response later in the text.Section 3.1 was concerned with a mean geopotential response at the 6.25 hPa level (Fig. 5).Mean anomalies (differences with reference) were found to be largest in the polar region and larger for the box GWD enhancements (both globally and locally) than for the corresponding ring enhancements.The important role of a purely constraint (from observations) meridional GWD component, especially for the polar vortex response, was highlighted.Most importantly, for simulations with the strongest GWD enhancement (SSWbox and SSWzon; Table 1), we observed different types of polar vortex events, namely a vortex split in response to the ring GWD enhancement 5 and a vortex displacement for a localized forcing (Fig. 6).
In section 3.2 we studied the influence of the artificial GWD and of its distribution on the PW activity.
We have found (Fig. 7) mean E-P flux convergence anomaly centered at the equatorward flank of the GWD enhancement area and an anomalous convergence in a broad area around 60°N in response to the artificial GWD.The anomalies are bigger for the box enhancements and in the box simulations we 10 have also identified anomalous, predominantly horizontal PW propagation indicative of in-situ PW generation.This is further supported by the results of FT analysis of the geopotential anomalies (Fig. 9), where, for the box simulations we have found especially the wave-1 and also wave-3 mean amplitude anomalously enhanced.Also, the short-term response (Fig. 10) showed the origin of the enhanced amplitudes to be in the GWD area.

15
Section 3.3 has been concerned with a residual circulation response.It was shown that there are significant differences in a zonal mean residual circulation between different distributions of the same zonal mean GWD (Fig. 11).A butterfly like pattern in the box-ring differences was identified centered at approximately 45°N (the center of the GWD region), with a stronger/weaker subsidence north/south of the enhancement region in the box simulations between 20 and 30 km log-pressure height.Evidence 20 was given that the artificial GWD in our model acts like a small obstacle for the flow, which was further supported by the 3D residual circulation analysis (Fig. 12).We have found downwelling to the northeast (downwind) and upwelling to the southwest (upwind) of the GWD box showing that GWs can contribute to longitudinal variations in the BDC.
The biggest limit of our analysis is naturally the artificiality of our GWD enhancement.The GWD 25 enhancement introduces an additional artificial constant momentum sink in the model.The concept of the artificial GWD enhancement leaves us also no chance to reflect any feedback between GWs and background conditions (changes in background winds, evolving PW field, etc.).Therefore, for example, our simulation of a vortex displacement differs from reality by not reflecting the background changes, as the GWs are known to be significantly filtered during SSWs (e.g.Holton, 1983; 30 Limpasuvan et al., 2012).Considering the intermittent nature of GWs (e.g.Hertzog et al., 2012;Wright et al., 2013), another inaccuracy of our sensitivity simulation set-ups arises from the constancy of the artificial GWD.In particular in the EA/NP region, where we expect mountain wave forcing to be prominent in January, variations of more than an order of magnitude from day to day are to be expected (Schroeder et al., 2009).A multiple (during a month) pulse like injection of the artificial GWD would 35 be arguably more realistic, but on the expense of absense of any steady response during the whole simulation.It is also a question, what is a more realistic illustration of the GW effect on the atmosphere, a sudden GWD injection or smooth increase and decrease with e.g. a 10-day e-folding time to minimize the initial adjustment noise as proposed by Holton (1983)?Also the spatial distribution of our artificial GWD is highly idealized (in both the horizontal and the vertical).We must note that we compare two "extreme" GWD longitudinal distributions only.It is also very likely that the sharp boundaries of the GWD enhancement in the 10box/zon and SSWbox/zon simulation are influencing some minor patterns of the response (e.g. the lee wave pattern in Fig. 10b).
In future work it is therefore necessary to take into account more realistic GWD distributions to address 5 e.g. the efficiency of PW creation.For example, it is possible that a configuration of GWD taking into account the EA/NP and e.g. the Greenland GW activity hotspot would favour enhanced wave-2 instead of wave-1 activity, and for comparison a chessboard-like or random distribution of GWD would possibly be more appropriate for comparison.Generally, the fact that the PW activity depends on the longitudinal GWD distribution (Fig. 7) suggests that the rate of compensation between resolved and 10 unresolved drag (Cohen et al. 2013(Cohen et al. , 2014) ) can be variable in dependence on the GWD distribution influencing the efficiency of PW creation.
Another motivation for future research is to concentrate on the position of the GW hotspots relative to the climatological stationary wave location in the stratosphere and to analyze the interaction between the GWD effects and the climatological waves.For example, the EA/NP hotspot lies in the region of the phase transition between a trough and ridge of the climatological wave-1 and our results show (Fig. 8) an anomalous amplification of wave-1 amplitude for a box GWD enhancement in this region.The importance of standing waves for polar vortex strength is well recognized (Watt-Meyer and Kushner, 2015;Yamashita et al., 2015).
In the atmosphere, the most natural, immediate and fastest way for communication of information in 20 the vertical are the GWs (apart from acoustic and acoustic-gravity waves with effects much higher in the atmosphere).We can argue that any change in the troposphere resulting in changes of sourcing, propagation or breaking conditions for GWs will almost immediately influence the distribution of GWD in the stratosphere, with possible effects demonstrated in our paper (in-situ generation of PWs in the lower stratosphere, anomalous vertical movements, etc.).For example, on the interannual scale, the occurrence and strength of the EA/NP GW hotspot can be dependent on the Pacific Decadal Oscillation (PDO) phase and can play a role in the relationship between PDO and SSW occurence frequency (Kren et al., 2015;Woo et al., 2015;Kidston et al., 2015).
There are more conclusions relevant for the SSW research in our results.It is common methodology (see e.g.Albers and Birner (2014) for a review of SSW preconditioning concepts) to estimate e.g. the relative impact of GWs and PWs on polar vortex preconditioning from zonal mean values of zonal forces only.But our results show that the dynamical effect of forcing depends also on its distribution.
The impact connected with a localized area connected with a higher value of drag can be much stronger than one would expect from the zonal mean value only.Importantly, we have found that for a sufficiently strong artificial zonal mean zonal force there is a vortex split response to the ring artificial

35
GWD configuration and vortex displacement for a localized forcing.We aim to investigate this in more detail and also for more realistic forcing distributions, but it seems to be clear at this stage that the SSW type may be determined also by the geometry of the forcing, not only by the vortex geometry.On the other hand, vortex geometry can to a large extent influence the distribution of the forcing, e.g.spontaneous emission processes connected with the jet (Plougonven and Zhang, 2014).
Blocking connection with SSWs is a well-known correlation (e.g.Andrews et al., 1987;Martius et al., 2009;Nakamura et al., 2014;Albers and Birner, 2014) but the mechanisms standing behind are still rather elusive.The geographical location and evolution of the stationary positive geopotential anomaly 5 with anomalous anticyclonic horizontal winds upstream of the GWD area is a remarkable feature of the atmospheric response to a localized GWD (Fig. 5) suggesting that GWs can be one of the missing mechanisms behind this relationship.This is connected with the important role of the meridional GWD component, especially for the polar vortex response.Interestingly, this feature becomes apparent for the localized enhancement only and has an almost negligible effect in simulations with ring 10 enhancements.To our knowledge, the effect of the meridional component of GWD on the middle atmospheric circulation has not been studied yet.Also, horizontal GW propagation is neglected in most climate model parameterizations (Kalisch et al., 2014).Thus, it is not surprising that there are only few modelling constraints regarding the horizontal propagation directions, although some information is available from ray tracing simulations (Preusse et al., 2009).In most studies based on satellite data,

15
GW propagation directions have not been analysed, because the information needed for such computation (e.g.hodograph analysis) is not available for most of the global observational instruments and their combinations (Wang and Alexander, 2010).
Finally, regarding polar vortex effects, the anomalous PW generation and breaking may be the physical justification for disturbing the vortex in its central levels which was a mechanism hypothesized by 20 Scott and Dritschel (2005).Traditionally, PWs are thought to be generated in the troposphere and propagate up on the polar vortex edge.But, as Scott and Dritschel (2005) pointed out, when wave amplitudes become large and nonlinear effects become important, the notion of upward propagation ceases to be appropriate.Therefore, they considered an option of some in situ disturbance at a given level, with a possible explanation being what we propose -localized GW breaking inducing anomalous 25 PW activity.
Regarding residual circulation, a general conclusion of this paper is that for the same magnitude of an artificial zonal mean zonal force (zonal mean meridional force as well) there are significant differences (depending on the magnitude of the GWD enhancement) in the zonal mean residual circulation between different distributions of this force (localized vs. zonally uniform).Also our results indicate 30 that the distribution of GWD may play a role in zonal asymmetries of the BDC.This is a clear signal that e.g. in the research of future BDC changes from climate models we need to be concerned not merely by the magnitude or latitude-height profile of the zonal mean GWD but also by its zonal distribution.In particular, the models should be able to mimic the main GW activity hotspots.This suggests the need for improvement especially of the nonorographic GW parameterization (though 35 nonorographic GW are usually assumed to have significant effect at higher altitudes than in the vertical range analyzed in this paper), since many global climate models use e.g. a globally uniform gravity wave source function (Geller et al., 2013).

3. 1
Fig. 5(1A) shows the mean (7th-30th day) horizontal wind and geopotential field at the 6.25 hPa level (13th model level) for the Ref simulation and the remaining plots in the first row show anomalies (i.e.differences from the results of this run) caused by different components of GWD with artificial values corresponding to the Box0.1 simulation.The second row (Fig. 5(2A) -Fig 5(2D)) shows horizontal Figs. 5(1B), 5(1C), and 5(1D) we see that among the GWD components modified in the Box0.1 simulation the response to the gcu component is the strongest.It induces a dipole structured anomaly with negative geopotential anomaly downwind from the region of GWD enhancement and positive anomaly north of this region (Fig. 5(1B)).The gt component alone induces a 20 positive anomaly of smaller magnitude northward and upstream of the area (Fig. 5(1C)).In contrast to that, meridional drag induces a negative geopotential anomaly northward and downwind of the area, which has the smallest magnitude of all three components, but is still significant (Fig. 5(1D)).
5(2C) and 5(2D)) show the same pattern as 30 induced by the meridional drag only (Fig. 5(1D)).The distribution of the response to the meridional component suggests that a box gcv enhancement in this geographical position can influence the geopotential response in the area of location of the Aleutian High.
15 geopotential anomalies at the 6.25 hPa level.FT provides information about the representation of different harmonics in the anomalous wave activity revealed by the E-P flux diagnostic, and about the spatiotemporal distribution of their amplitudes.The mean (7 th -30 th day) latitudinal structure of reference amplitudes of leading PW modes is given in Figs.9a, b.Anomalous amplitudes (Box0.1-Refsimulation) are presented in Fig. 9c, d and differences from the Zon0.1 simulation are shown in Fig 9e,

5 Fig. 1 :
Fig. 1: Mean January zonal means of temperature (a), zonal wind (b), GW induced heating (d) GW induced zonal wind (e) and meridional wind (f) acceleration for the reference simulation.Additionally, January mean zonal mean nudging strength for the strongest GWD injection (SSWbox simulation in Table 1.) is shown (c).