Atmospheric Circulation Regimes in a Nonlinear Quasi-Geostrophic Model

Atmospheric low-frequency variability and circulation regime behavior are investigated in the context of a quasi-geostrophic (QG) three-level T63 model of the wintertime atmospheric circulation over the Northern Hemisphere (NH).Themodel generates strong interannual and decadal variability, with the domination of the annular mode of variability. It successfully reproduces a satisfactory model climatology and the most important atmospheric circulation regimes.The positive phase of the Arctic Oscillation is a robust feature of the quasi-geostrophic T63model.Themodel results based on QG dynamics underlie atmospheric regime behavior in the extratropicalNHand suggest that nonlinear internal processes deliver significant contribution to the atmospheric climate variability on interannual and decadal timescales.


Introduction
The atmospheric circulation possesses preferred quasistationary planetary-scale flow configurations, the so-called circulation regimes.The existence of these regimes was presumed several decades ago [1][2][3][4][5] and statistically proven later on (e.g., [6][7][8][9]).According to a modern paradigm [9][10][11], the atmospheric circulation regime behavior, owing to irregular transitions between the regimes, substantially contributes to atmospheric natural variability on interannual and decadal timescales.However, understanding of dynamical mechanisms underlying these regimes remains a challenge to dynamical meteorologists and climate theorists, despite the significant efforts that have been made.This paper continues investigations by Sempf et al. [12,13], where the atmospheric regime behavior was investigated by means of a quasi-geostrophic (QG) three-level T21 model of the wintertime atmospheric circulation over the Northern Hemisphere (NH), and extends the research with a higher-resolution hemispheric QG three-level T63 model.In Kurzke et al. [14] the Southern Hemisphere (SH) version of the three-level quasi-geostrophic T21 model was coupled to a global ocean circulation model with 2 ∘ × 2 ∘ resolution and simplified physics, and the consequent SH wintertime extratropical climate variability has been studied.Here a step toward increasing the horizontal spectral resolution of the atmospheric model has been made by enhancing the resolution to T63.This resolution is often used by comprehensive climate models, for example, by ECHAM5 [15].Because the planetary-scale, quasi-barotropic flow configuration is primarily studied in the context of atmospheric regimes, a coarse vertical resolution with three levels may be conjectured as adequate for capturing essentials of the regime behavior.An increase in horizontal resolution alone can be considered as a meaningful sensitivity study to check the robustness of model simulation with respect to atmospheric circulation regimes.The T63 model used thus serves as an idealized tool that avoids substantial complexities (including computational costs) associated with multilevel model structure and allows us to focus primarily on the effects of an increased horizontal resolution comparable to medium-range models [15].A similar approach was used by Zidikheri and Frederiksen [16] who in their study of stochastic subgrid parameterizations applied a global QG two-layer model at three different resolutions: T126, T63, and T31.In our paper we follow similar lines of thinking and investigate the extent to which a three-layer QG model with substantially increased horizontal spectral resolution reproduces the essential features of low-frequency variability and atmospheric circulation regime behavior.

Model Description
The QG model, described in detail by Weisheimer et al. [17] and Sempf et al. [12,18], simulates the spatiotemporal evolution of the NH large-scale atmospheric circulation in terms of atmospheric stream functions in three homogeneous vertical layers of equal mass under perpetual winter conditions.In the following, we refer to the lowest layer as the 833-hPa level, to the middle layer as the 500-hPa level, and to the upper layer as the 167-hPa level.This three-level model has the minimum vertical resolution necessary to represent interactions between the troposphere and the lower stratosphere.It is a hemispheric model with a T63 spectral resolution.The spectral interaction coefficients are found by transformation from the spectral space onto a grid with 192 points in the zonal direction, which corresponds to a grid spacing of 1.875 ∘ in longitude and with 96 points in the meridional direction (the Gaussian latitudes).The latter points are symmetric with respect to the equator; that is, there are 48 points in the NH, but the distance between them varies according to the relative weight of different Gaussian latitudes.The Northern Hemisphere's topography adapted to the spectral T63 resolution, displayed in Figure 1, provides the orographic forcing.Diabatic heating is established by thermal relaxation towards predefined radiative equilibrium temperature fields at the model's 333 and 667-hPa pressure levels.The relaxation timescale is 22.7 days.An additional surface friction forcing mechanism damps the 833-hPa stream function towards a predefined zonally symmetric surface forcing function [19] with a timescale of 1.08 days.This zonally symmetric forcing mimics surface baroclinicity due to the meridional gradient of surface temperature and helps to strengthen low-level westerlies that would be too weak otherwise.A horizontal scale-selective ∇ 6 hyperdiffusion attenuates the shorter waves.Those with maximal total wavenumber  = 63 are damped with an e-folding time of 48.5 h (≈2 days).The vertical temperature lapse rate has been fixed to 3.0 K km −1 at 333 hPa and 6.5 K km −1 at 667 hPa.
The nonzonal components of the radiative equilibrium temperature fields have been adjusted in a way that, in the wintertime mean, leads to realistic patterns of non-zonal extra-tropical diabatic heating.The zonal components of the radiative equilibrium temperature fields and the surface forcing have been tuned to generate zonal wind climatology as realistic as possible in order to ensure adequate westerly flow against orography and appropriate flow instability conditions.After tuning, the diabatic heating field in this simplified model is kept constant over time.
The adjustments of forcing are made by an automatic iterative procedure, described by Sempf et al. [12,13,18].The nonzonal parts of wintertime (December, January, and February) diabatic heating rates over the Northern Hemisphere for isobaric surfaces at 300 hPa and 700 hPa are derived from observations by Nigam et al. [20] and have been attenuated near the equator (Figure 2).Compared with the T21 model of Sempf et al. [12,13,18], the magnitude and structure of the nonzonal thermal forcing fields remain very similar.Diabatic warming occurs over the oceans and cooling over the land regions.The zonal thermal forcing fields have the same structure too, but the amplitude of the radiative equilibrium temperature has increased following the enhanced resolution, while the amplitude of the surface forcing has been decreased.This change is a result of the increase in energy in the model system, since during the course of model development all other parameters of the T21 model were adjusted to the higher T63 resolution to ensure realistic model simulations.The adjusted parameters include the scale-selective horizontal diffusion and the surface friction coefficient.The intensity of surface friction has a large impact on the model behavior as shown for the T21 model version by Sempf et al. [12].This sensitivity is observed for the T63 model as well and is studied here in terms of changes of atmospheric regime behavior in Section 3.3.To summarize, for the implementation of the T63 QG model all parameters and forcing fields have been selected based on their ability to represent the zonal wind climatology as realistic as possible in a similar way as it was done for the T21 model [18].Whether this common way of model tuning provides consistent subgrid scale parameterizations with respect to a correct representation of the nonlinear energy and enstrophy cascades as it is derived in Zidikheri and Frederiksen [16] has not been studied yet.

Model Results
Using the forcings described in Section 2, a 550-year model integration was performed.This section describes the model  climatology, followed by the evaluation of the low-frequency variability based on the results of empirical orthogonal function (EOF) analyses.Finally, the dominant atmospheric circulation regimes are analyzed.To compare the model results with the observations, we use the geopotential height fields computed from NCAR-NCEP reanalysis data [21], averaged with mass weighting over three vertical layers of approximately equal mass, and their statistical analyses by Sempf et al. [12,13].

Model Climatology.
The time-mean zonal wind profiles of the model simulation are displayed in Figure 3, together with zonal winds taken from the NCEP-NCAR reanalysis data, which were brought to the three-layer structure with the above-described averaging procedure and subsequently averaged over time for the winter (DJF) months.The agreement between the modeled and observed profiles is very good from the equator to the mid-latitudes.Modeled winds from the polar regions are slightly weaker than those observed.This can be attributed to the formative effect of decadal scale variations during the long-term integration, since after 1000 days of model tuning (see Section 2) the modeled and observed zonal wind profiles were more closely matched.The wind maxima at the middle and upper levels represent the subtropical jet, and the strongest wind maximum occurs in the upper troposphere at 30 ∘ N.
For consistency with the reanalysis data, the daily stream function output data for each model level were converted into time series of geopotential height fields by means of the nonlinear balance equation (see also, e.g., [22][23][24]).The diagnostic balance equation is derived from the divergence equation by applying scale analysis.It is common to convert the stream function data from QG models into geopotential height data by means of the balance equation (e.g., [25]).The averaged geopotential height fields over the simulated  The model correctly reproduces key characteristics in the observed geopotential height field, as can be found at the lower model layer over each side of the Rocky Mountains and the Himalayas.There is a trough on the leeward side of both mountain systems, which is shifted westward over North America compared with the observations.Another pronounced trough is present over Europe, which like the other two extends to the upper layer but attenuates with height.As in the observations, Icelandic and Aleutian minima are recognizable at the lower level, but the latter minimum is rather weak.This model bias has been shown already by the spectral T21 model version [12,13], but for the T63 model the Aleutian low is even weaker.At the middle and upper tropospheric levels the observations show minima near the North Pole.In the model, these minima moved to the south and extended across the northern Siberian coast of the Kara Sea up to the Laptev Sea.

Advances in Meteorology
Figure 5 displays the standard deviations of geopotential height at the three model levels (left) and for the layeraveraged NCEP reanalysis data (right).While the maxima of standard deviation of the observed geopotential height on the lower and middle levels are concentrated in the vicinity of the Aleutian and Icelandic lows, the latter one does not exist in the model.The model Aleutian standard deviation maximum is strongly spatially pronounced and extends over the entire North Pacific, being larger than in the observations.
In Figure 6 we display the simulated meridional eddy heat flux at 333 hPa (a) and 667 hPa (b) of the 550-year integration.The model yields the maximum of the heat flux at the right position of 40 ∘ N, although its magnitude, which is equivalent to a total meridional heat flux of about 10 PW, is two times stronger than in reality (e.g., Figure 1 in [26]).We attribute this difference to the specific form of the applied zonal thermal forcing, which was used to reproduce a realistic mean zonal wind structure over the NH hemisphere.Also, as the model has only three levels, a more realistic simulation of the heat flux amplitudes requires an improved vertical resolution of the atmospheric circulation as well.

Low-Frequency Variability.
The spatiotemporal structure of the modeled low-frequency variability is studied by means of an empirical orthogonal function (EOF) analysis of the geopotential height data for the whole hemispheric domain (e.g., [27,28]).Equal-area weighting is ensured by multiplying all fields by the square root of the cosine of latitude [29].Here, we show renormalised EOF patterns carrying the physical units and corresponding dimensionless time series of leading principal components (PC) (cf.[30]).To determine the low-frequency variability patterns of geopotential height fields, 30-day averages of the modeled fields and monthly averages of the reanalysis data fields were analyzed.For the reanalysis, the data preprocessing included calculation of anomaly fields by removing the mean seasonal cycle.To quantify the similarity between EOF patterns spatial correlation coefficients have been calculated.
The first EOF patterns obtained for the modeled levels are very similar to the well-known patterns from NCEP reanalysis data (Figure 7, right).Small departures from patterns shown, for example, in Kimoto and Ghil [6], Wu and Straus [31], and Itoh and Harada [32], are due to the applied layer averaging described above.Throughout the whole troposphere, the dominant patterns display an annular structure.However, as mentioned by, for example, Molteni et al. [33], the first mid-tropospheric EOF (Figure 7, right, middle) has more pronounced zonally asymmetric features than the first EOF patterns below and above.The zonally asymmetric features indicate wave train structures like the PNA pattern over the eastern Pacific Ocean.The stronger deviations from zonal symmetry in the mid-troposphere are also indicated by the spatial correlation coefficients between the levels (cf.Table 1), with higher correlation between the upper and lower levels (0.88) than between the upper and middle levels (0.79).The modelled patterns are slightly southward shifted compared to reanalysis data, especially over the Atlantic.The first EOF of the model (Figure 7, left) is described by a pronounced annular structure at all three levels.The EOF pattern at the upper model level is similar to the observed one with high pattern correlation of 0.93 (see Table 3) but is slightly more zonally symmetric and explains  a larger percentage of variance than its observed counterpart (cf.Table 2).In particular, the central vortex is more elongated in SW-NE direction in the reanalysis data.As in the observations, the first mid-tropospheric EOF (Figure 7, left, middle) has more pronounced zonally asymmetric features than the first EOFs below and above, which is quantified again by the higher correlation between the upper and lower levels (0.76) than between the upper and middle levels (0.72, see Table 1).In particular, the observed patterns at the middle level display stronger deviations from the zonal structure.The location and strength of the mid-latitude centers of action show larger differences between reanalysis and model data, in particular in the middle level and over the North Pacific corresponding to a decreased pattern correlation between the observed and modelled patterns of 0.63 (see Table 3).
In the lower and middle levels, the percentage of explained variance by the annular pattern is reduced compared to the observations (see Table 2).The annular pattern visible both in the observations and in the model corresponds to the Arctic Oscillation (AO) pattern inferred from EOF analysis of the wintertime sea-level pressure field at the sea level [34,35].
The second EOFs for the reanalysis data are shown in Figure 8 (right).The pattern for the upper troposphere displays a wavenumber-1 pattern, which is typical for stratospheric second EOFs (e.g., [32]).In accordance with Molteni et al. [33] and Wu and Straus [31], the low-and midtropospheric patterns are similar to the Cold-Ocean/Warm-Land pattern (COWL [36,37]).The COWL pattern is obtained as a regression pattern upon NH mean SAT anomalies over land and results from contrast in thermal inertia between continents and oceans [37].Over the Atlantic, the COWL pattern includes a north-south dipole similar to the North Atlantic Oscillation (NAO) pattern.The change from the COWL pattern in the troposphere to the wavenumber-1 pattern in the lower stratosphere is indicated by decreased pattern correlations between the levels (0.81 between lower and middle levels versus 0.35 between the middle and the upper levels).
The second EOF of the model at the upper level (Figure 8, left) shows a wavenumber-1 pattern similar to the observations but slightly shifted to the east, resulting in a pattern correlation with the observed pattern of 0.51 (see Table 3) and similar value of explained variance (see Table 2).At the other two levels the model differs compared to the observed patterns and shows only low pattern correlations with the observed patterns (Table 3) and reduced explained variance (Table 2).The modeled pattern does not resemble the COWL pattern but is characterized by two main centers of action, located over East Siberia and middle Asia.
To summarize, the observed annular mode and its vertical changes are partly reproduced by the quasi-geostrophic model with increased T63 resolution.Quantified by spatial correlation patterns between modeled and observed EOF patterns, there is an improvement for the T63 compared to the T21 resolution, in particular in the upper and middle levels (Table 3).Moreover, for the T63 resolution, there is high correlation between the spatial structures of EOFs in the middle and the lower model levels.This high correlation between the two lower levels is not restricted to the leading EOF but occurs also for higher EOFs (Table 1).Compared to the reanalysis data and to the T21 model, the quasibarotropic behavior in the lower and middle troposphere is more pronounced (cf.Table 1).This is most likely due to the increased T63 spectral resolution of the model and can be explained by the additional energy created at the lower model level and afterwards upward energy propagating towards the middle level.
The spatial structure of model variability is determined by the annular mode.To assess the extent to which the temporal behavior is associated with the annular mode, or with another spatial pattern of variability, the time series of leading principal components (PCs) are analyzed by performing wavelet analyses.Figure 9 (left) shows the continuous local and global wavelet spectra of the first PC (PC 1) at each of the three model levels.For the local wavelet spectrum, the contour lines give information on the relative power at a certain period and a certain time, whereas the global spectrum gives information on the time-integrated power at a certain period.At the lower and middle level, the annular mode presents bands of enhanced energy at 7-10 yrs and 20-30 yrs.At the middle level, the enhanced energy on the decadal scale results in a statistically significant spectral peak in the global spectrum.In contrast, at the upper level, several energy bands are present, which sum up to significant peaks in the global  spectrum in the range of 2-3 yrs and 7-14 yrs and at about 35 yrs and 70 yrs.
Compared to the wavelet analyses of the observed NAMindex given in Rossi et al. [38], the model reproduces the significant fluctuations on a decadal scale (7-10 yrs), whereas the fluctuations on interannual scale (2-3 yrs, 4-6 yrs) are underestimated.This underestimation can be related to the time-constant diabatic heating fields which force our model simulation.According to Frederiksen and Branstator [39], a substantial part of interannual variability in the atmosphere is induced by interannual variability of the SST.The variability on multidecadal scales (30-50 yrs, 70-90 yrs) is strongest at the upper model level.
In the wavelet spectra of the second PC (PC 2) maxima at 7-10, 14, 35, and 70 years occur (Figure 9, right).These significant levels of variability extend over almost the entire model integration.The period of 8-9 years is mostly pronounced during the second half of the model integration.The wavelet spectra at the middle and lower levels are again very similar, the signal being stronger at the middle level.Whereas only the period of 8-9 years is statistically significant in the global spectrum at both heights, the period of 14 years is additionally significant at the middle level.
Wavelet spectra of PC 1 show a strong vertical coupling of the lower and middle levels, whereas for PC 2 this similarity is less pronounced.For both PCs the lower stratospheric level is more decoupled and the multidecadal variability dominates there.At all heights, the temporal variability of the leading patterns is characterized by strong fluctuations on decadal scales (7-10 yrs), which indicates the importance of internally generated variations by the nonlinear atmospheric dynamics at these timescales.
The results show that the model atmosphere, in the absence of any external time-dependent forcing but solely due to internal nonlinear dynamics, generates pronounced variability on ultralong timescales.This confirms the robustness of results of previous and more idealized model studies [12,13,17,[40][41][42][43][44] for a horizontal resolution of T63, also applied in state-of-the-art GCMs.
The real atmosphere is, of course, not a closed dynamic system, such as the model atmosphere, but is influenced by many external factors, such as the solar radiation, anthropogenic influences, or interactions with the cryosphere and the oceans.Thus, the decadal variability determined by an analysis of observed time series of the Arctic Oscillation [38,45] arises due to internally generated and externally forced processes and their interactions.Using our model as an analogue of the real atmosphere, where the decadal variability associated with the annular mode is found at all model levels, it can be conjectured that the nonlinear internal dynamic processes make significant contributions to the atmospheric climate variability on decadal timescales.

Regime Behavior.
To analyze the modelled atmospheric circulation regimes, the regime detection in the twodimensional state space spanned by the first two PCs at the lower model level has been performed as in Kimoto and Ghil [6] (see also [12,46]).Areas of the state space with an unexpected high recurrence probability are related to regimes.Therefore, the two-dimensional probability density function (PDF) has been calculated.Areas with significantly higher probability density than what to be expected are estimated by means of Monte Carlo simulations of bivariate autoregressive processes of order 1 [AR(1) processes] with the same means, variances, and coefficient matrix of the AR(1) process as the bivariate original PC data.For the calculation of the coefficient matrices of the AR models and the simulation of the Monte Carlo PC surrogates the MATLAB package ARfit has been used [47].Finally, spatial anomaly patterns of the recurrent climate regimes have been reconstructed from the PC 1-PC 2 coordinates of the centre of the estimated areas with an unexpected high recurrence probability.The probability density calculation was performed with the 30-day averaged data for the entire model integration.The comparison with the real atmosphere is based on the layer-averaged and then monthly averaged NCEP-NCAR reanalysis data for the winter months (DJF) described above (55 years).Figure 10 (left) shows the results of the probability density estimates for the model run.The thick lines enclose areas with unexpected high recurrence probability (at the 95% level) and thus mark the atmospheric circulation regimes.For further analysis, four areas (labeled as A, B, C, and D in Figure 10) are considered.The corresponding regime anomaly patterns are displayed in Figure 11.
Regime A is related to a pattern with positive geopotential height anomalies over East Siberia and over the Chukchi, East Siberian, and Laptev Sea, accompanied by negative anomalies over the Pacific Ocean and the Asian mid-latitudes and subtropics.Related to the mean geopotential height field at the lower model level (Figure 4), this implies a weakening of the meridional pressure gradient over the North Pacific and Asia, leading to stronger meridionalization of the atmospheric flow over the North Pacific and Asian midlatitudes.The counterpart of regime A is detected as regime D. Due to the decreased amplitude of the related anomaly pattern seen in Figure 11, the resulting zonalization of the flow over the North Pacific and Asia is weak.
Regime B is characterized by four anomaly centers, over the northern North Pacific (Bering Sea), West Russia, Greenland/Iceland, and the subtropical North Atlantic.The latter two centers resemble the NAO in its positive phase.Due to its well-expressed projection onto EOF1, regime C resembles the positive phase of the AO (AO+) with hemispheric-wide increase of the zonal flow components over the mid-latitudes.Analysis of the reanalysis data reveals four regimes as well, labeled in Figure 10 (right).Due to the rather small sample size of 165 (53 winters) a confidence level of 85% has been applied to detect the four regimes.The corresponding geopotential height anomaly fields are shown in Figure 12.All regimes have well-expressed projections onto either EOF1 or EOF2; therefore, the regimes show up as the AO in its two phases (AO+ as regime A and AO− as regime C) and the COWL pattern in two phases (COWL+ as regime D and COWL− as regime B).
To summarize, the regime analyses for the reanalysis data as well as for the control run data have given evidence for the existence of recurrent atmospheric circulation regimes.For each data set we detected four regimes.Regime C displays well-expressed projections onto the annular mode pattern in its positive phase for both data sets.
We performed two additional 550-yr model runs: the first experiment with increased surface friction, now with a timescale of 0.91 days, and the second experiment with decreased surface friction characterized by a timescale of Figure 13: (Left) estimated probability density in the space spanned by the first two PCs of the 833 hPa geopotential height for experiment I with a surface friction timescale of 0.66 days.Contour interval is 0.02; zero contour omitted.The thick lines enclose areas with significantly higher probability than to what be expected in the case of adapted red noise process (at the 90% significance level) and thus mark the circulation regimes.The regimes are denoted as indicated in the text.(Right) as left figure, but for experiment II with a surface friction timescale of 1.73 days and 99% level of significance.
1.73 days.The resulting PDFs for experiments I and II are shown in Figure 13, left and right, respectively.The areas with unexpected high recurrence probability have been determined with significance levels of 90% for experiment I and 99% for experiment II.Dramatic enhancement of the regime behavior accompanying the surface friction reduction in experiment II is observed and in agreement with Sempf et al. [12].
As in the control run, the regime analysis for experiment I reveals four regimes as well, labeled in Figure 13 (left).The corresponding regime anomaly patterns are shown in Figure 14.The anomaly pattern of regime A of experiment I bears similarity with the pattern of regime A of the control run, but with increased zonal symmetry.Thus, the pattern implies a hemispheric-wide weakening of the meridional pressure gradient between polar and subtropical latitudes leading to stronger meridionalization of the atmospheric flow over the mid-latitudes.The counterpart of regime A is detected as regime C. Due to the decreased amplitude of the related anomaly pattern seen in Figure 14, the resulting zonalization of the mid-latitude flow is weak.The remaining two regimes B and D are counterparts of each other as well.Regime B is characterized by three anomaly centers, which form a dipole over the Western Hemisphere leading to decreased meridional pressure gradients and decreased zonal flow.The third anomaly center indicates positive anomalies over Siberia related to a strengthening of the Siberian High.Conversely, the anomaly pattern of regime D is related to increased zonal flow over the Western Hemisphere and a weakening of the Siberian High.
For experiment II, three pronounced regimes have been detected.The corresponding anomaly patterns are displayed in Figure 15.They are characterized by strong zonally symmetric structures related to the strong zonally symmetric character of the two leading EOFs spanning the state space.The annular structures of regimes A and C correspond mainly to the AO− and AO+ regimes, respectively.Regime B has strong projections onto EOF2 and displays an annular structure as well.In contrast to regimes A and C, the annular structure of regime B is restricted to the middle and high latitudes.This pattern corresponds to the Polar Annular Mode (PAM) [48].Black and McDaniel [48] detected the PAM in reanalysis data as second EOF of the zonally averaged zonal wind fields.The fact that not only the first dominant pattern but also the second dominant pattern reveals an annular structure underlines the enhanced zonality as a result of reduced friction for experiment II.
The differences between the regime anomaly patterns of the three model experiments are due to different positions of the regimes in the respective state spaces and due to pronounced changes of the vector (EOF1, EOF2) spanning the state space for the regime detection.In particular, the spatial structure of EOF2 changes considerably by reducing the surface friction.In contrast, the spatial structure of EOF1 remains nearly unchanged, but the amplitude of this annular pattern increases (see [49]).
The dynamical origin and the dramatic enhancement of the regime behavior by reducing the surface friction have been studied by Sempf et al. [12,13] using the same model with T21 resolution.These studies revealed that the circulation regimes emerge from the unification of multiple attractors.The results presented here support the hypotheses that changes in the large-scale geometry of the model attractors due to changes of the friction parameter can explain the detected atmospheric regime behavior.

Summary and Discussion
A QG, hemispheric, three-level model using the T63 orography of the NH has been forced by an adjusted thermal forcing using an automated iterative procedure, described in detail in Sempf et al. [12,13,18].The zonal part of the thermal forcing has been adjusted to produce a realistic zonal wind structure, whereas the nonzonal part of the thermal forcing has been tuned in such a way that the time-averaged nonzonal extratropical diabatic heating acting in the model coincides with wintertime observations.In approximately 550 yrs of perpetual winter simulations, the model has been shown to reproduce the extratropical largescale wintertime circulation with adequate accuracy, though variability forced by interannually varying diabatic heating sources is not included in this kind of model simulations (cf.[39]).The model tuning appeared substantially more difficult than in the T21 model version and though the result shown in Figure 3 is satisfactory it is not as perfect as that in the T21 model (Figure 3 in [12]).Whether the applied model tuning provides consistent subgrid scale parameterizations between the T21 and T63 model versions and ensures a correct representation of the nonlinear energy and enstrophy cascades as it is suggested in Zidikheri and Frederiksen [16] has not been studied yet and is a topic for future research.
Using a 550-yr perpetual wintertime integration, the model shows multiple circulation regimes in the space spanned by the first two PCs of the 833-hPa geopotential height, which are discriminated as areas with unexpected high probability, that is, significantly higher probability than what to be expected in the adapted red noise process (at the 95% level).The four most pronounced regimes resemble the AO+, AO−, NAO+, and PNA regimes, as observed in the real atmosphere.In the reanalysis data, by using basically the same methodology, four regimes, AO+, AO−, COWL+, and COWL−, have been identified.The T21 model by Sempf et al. [12] shows two circulation regimes (AO+ and NAO−) in the same (PC 1, PC 2) space.It can be concluded that our model results represent a satisfactory model climatology and indicate that circulation regimes, especially the AO+ pattern, are robust features of the QG model and survive a very substantial increase in the horizontal spectral resolution to T63.
We have shown that our T63 three-level QG model not only successfully reproduces the most important NH circulation regimes, but also depicts a significant level of variability on interannual and decadal timescales due to inherent nonlinear atmospheric dynamics.A prominent feature of the model variability is the domination of the annular mode and the high degree of linkage between the dynamical processes at the two lower model levels, the upper level behaving more independently.On the one hand, such an exaggerated quasi-barotropic behavior at the two lower levels might be a model artefact.On the other hand, an increased horizontal spectral resolution in the T63 model permits the reproduction of nonlinear dynamical processes like vorticity advection with higher accuracy than in models with coarser spectral resolution, because significantly more degrees of freedom are explicitly involved in model dynamics.
The most obvious explanation for the deterioration effect mentioned at the beginning of this section consists in necessity of enhancement of vertical resolution together with increased horizontal resolution.According to Lindzen and Fox-Rabinovitz [50], the consistent combination of horizontal (Δ) and vertical (Δ) scales for quasi-geostrophic flows is derived from the relationship Δ = (/ 0 )Δ, where  0 is the reference value of Coriolis parameter and  the Brunt-Väisälä frequency.Lindzen and Fox-Rabinovitz [50] emphasized that an insufficiently fine vertical resolution leads to incorrect solutions, because even if the vertical resolution may be adequate for the horizontal scales one is physically concerned with, smaller scales are inevitably generated in the course of integration.This argument can be reformulated as the requirement of approximate equality of the magnitude of the relative vorticity and the baroclinic stretching term in the mathematical expression for quasigeostrophic potential vorticity (QGPV), the latter being a linear elliptic operator applied to the geostrophic stream function (e.g., [51][52][53]).Restriction to three levels, as in our case, is equivalent to limiting to the barotropic and two first baroclinic vertical modes.In fact, the latter introduces biases in reproducing the phase velocities of smaller horizontalscale baroclinic Rossby waves but does not make the problem ill posed.By substantially increasing the spectral resolution in the model to T63, we approach the intrinsic limit of applicability of the QG approximation, because the Rossby number Ro = /( 0 Δ) ( is the characteristic wind speed weakly dependent on the spatial scale of motion) on the smallest resolved scale becomes comparable to unity.From this perspective, keeping a coarser vertical resolution than would follow from the consistency relationship permits the decoupling of different vertical levels for the motions of the smallest horizontal scales, when the thermal wind equation, which is physically responsible for this coupling, ceases to be valid.In effect, the smaller-scale motions obey barotropic dynamics at each model level, whereas it is known (e.g., [51,52]) that a barotropic model is in fact applicable under the quasi-solenoidal approximation when the Mach number is much less than unity, which is significantly weaker requirement than the smallness of the Rossby number, and therefore is nearly always filled.As such, we have a seamless transition of quasi-geostrophic dynamics operating at larger baroclinic scales, which are satisfactorily described within the framework of our three-level model, to barotropic dynamics acting on the smaller scales.The latter scales nonlinearly interact through the vorticity advection with larger-scale motions.Due to the aforementioned model performance limitations, enhancing vertical resolution is desirable.This is planned in future work but, in view of the above considerations, it will be accomplished to a lesser extent than what follows from the ideal Δ = (/ 0 )Δ relationship.

Figure 2 :
Figure 2: Wintertime (DJF) extratropical nonzonal diabatic heating (K day −1 ) at (left) 300 hPa and (right) 700 hPa derived from NCEP-NCAR reanalysis data and used in the QG three-level model.See text for details.

Figure 3 :
Figure 3: Modeled and observed time-mean zonal wind profiles (m s −1 ) for the three model levels and corresponding vertical layers, respectively.

Figure 9 :
Figure 9: (First column) global and (second column) local wavelet power spectra for the monthly averaged PC 1 of the geopotential height of the (top) 167 hPa, (middle) 500 hPa, and (bottom) 833 hPa model level.The contour lines give information on the relative power at a certain scale and a certain time.At both ends, dash-dotted lines separate regions where edge effects become important.The thick black contour envelopes areas of greater than 95% confidence for a corresponding red noise process.(Third column) global and (fourth column) local wavelet spectra for the monthly averaged PC 2 of the geopotential height of the (top) 167 hPa, (middle) 500 hPa, and (bottom) 833 hPa model level.

Figure 10 :
Figure10:(Left) estimated probability density in the space spanned by the first two PCs of the 833 hPa geopotential height field.Contour interval is 0.02; zero contour omitted.The thick lines enclose areas with significantly higher probability than what to be expected in the case of adapted red noise process (at the 95% significance level) and thus mark the circulation regimes.The regimes are denoted as indicated in the text.(Right) as left column, but for the layer-averaged NCEP-NCAR reanalysis data with significance at the 85% level.

Figure 11 :
Figure 11: Composite geopotential height (gpm) anomaly patterns of the modeled regimes of the control run (areas A-D in Figure 10 left) at the 833 hPa level.

Figure 12 :
Figure 12: Composite geopotential height (gpm) anomaly patterns of the observed regimes (areas A-D in Figure 10 right) at the 833 hPa level.

Figure 14 :
Figure 14: Composite geopotential height (gpm) anomaly patterns of the modeled regimes of experiment I (areas A-D in Figure 13 left) at the 833 hPa level.

Figure 15 :
Figure 15: Composite geopotential height (gpm) anomaly patterns of the modeled regimes of experiment II (areas A-C in Figure 13 right) at the 833 hPa level.

Table 1 :
Spatial correlation coefficients of dominant variability patterns (EOF1 and EOF2) between the atmospheric levels for reanalysis data and the control runs of the quasi-geostrophic model with T63 and T21 resolution.

Table 2 :
Explained percentage of total variance by the dominant variability patterns (EOF1 and EOF2) for reanalysis data and the control run of the quasi-geostrophic model with T63 resolution.

Table 3 :
Spatial correlation coefficients of dominant variability patterns (EOF1 and EOF2) between the reanalysis data and the control runs of the quasi-geostrophic model with T63 and T21 for the atmospheric levels.