## Abstract

Idealized model simulations have long been established as valuable tools to gain insight into atmospheric phenomena by providing a simplified, easier to comprehend version of the complex atmospheric system. A specific subgroup of idealized simulations, such as baroclinic channel models, requires the initialization of the model with balanced atmospheric fields to investigate the evolution of an introduced perturbation. The quality of these simulations depends on the degree of balance of the initial state, as imbalances result in geostrophic and hydrostatic adjustment processes that potentially skew the results. In this paper, a general method to create geostrophically and hydrostatically balanced initial conditions is introduced. The major benefit of this method is the possibility to directly define a basic state wind field with the pertinent atmospheric fields being derived given appropriate boundary conditions. Application of the method is exemplified by constructing initial conditions for a baroclinic test case with WRF and analyzing a perturbed and unperturbed numerical simulation. The unperturbed simulation exhibits weak inertia–gravity wave activity and minimal adjustment of the initial state during a 5-day simulation, which confirms the high degree of initial balance provided by the initialization technique. In the perturbed simulation, baroclinic instability is initiated, resulting in a cyclogenesis event similar to previous idealized baroclinic channel simulations. The proposed method is compared with initial conditions formulated in a Boussinesq framework, illustrating the difference in imbalances and their effect on perturbation growth.

## 1. Introduction

Although state-of-the-art numerical weather prediction models have become extremely valuable tools for weather forecasting and research, these models often fall short of providing an optimal platform to disentangle the complicated interactions between the different physical and dynamical processes. This dilemma is often circumvented by employing a model hierarchy (Held 2006), where one aims to reduce the complexity of the model system and creates an appropriate setup for the respective problem. A specific, widely applied subgroup of idealized numerical simulations is the so-called baroclinic channel setup consisting of a baroclinic zonal mean state and a perturbation to initiate development. Extensive usage of this type of idealization includes, for example, the study of cyclogenesis and frontogenesis (e.g., Hoskins and West 1979; Thorncroft et al. 1993; Wernli et al. 1999; Schemm et al. 2013), the formation of jets (Baker et al. 2013), and chemical transport (Polvani and Esler 2007). For these types of studies, it is essential that the initial conditions ideally are in perfect balance to allow for an unambiguous analysis of the development of interest. In this study, we present an initialization technique allowing for an improved balance of the initial fields in a zonally symmetric setup.

Existing approaches to define a basic state can be divided into three main categories. The first category is using potential vorticity (PV) inversion, which has the advantage to implicitly demand hydrostatic and geostrophic balance (e.g., Rotunno et al. 1994; Zhang 2004; Tan et al. 2004; Schultz and Zhang 2007; Plougonven and Snyder 2007; Waite and Snyder 2009). However, this benefit comes with the disadvantage of reduced control over the wind and temperature field, rendering it difficult to predefine the vertical stability and horizontal or vertical wind shear. Furthermore, PV inversion often requires computationally more complex calculations and the imposed balance often requires iterative methods. The second category is integrating thermal wind balance, which allows for a greater degree of control over the temperature and wind field while keeping the computation rather simple and straightforward (e.g., Balasubramanian and Yau 1996; Wernli et al. 1999; Yamazaki and Peltier 2001; Conzemius et al. 2007; Beare 2007; Polvani and Esler 2007; Riviére 2008; Boutle et al. 2010; Wang and Polvani 2011). The third category is obtaining the three-dimensional pressure distribution by integrating hydrostatically upward and downward, given a horizontal pressure distribution at a midtropospheric level as well as the three-dimensional temperature field (Nuss and Anthes 1987; Olson and Colle 2007; Schemm et al. 2013). A variation of this approach is used by Chuang and Sousounis (2000), who define the horizontal pressure distribution at the surface instead of the midtroposphere.

To initiate development, baroclinic channel models often utilize a small-amplitude, normal-mode disturbance superimposed on a zonally uniform basic state (e.g., Simmons and Hoskins 1978; Thorncroft et al. 1993; Rotunno et al. 1994; Balasubramanian and Yau 1996). This type of setup captures synoptic-scale cyclogenesis as zonally periodic growing perturbations. Another type of setup features finite-amplitude perturbations to initiate cyclogenesis. The latter perturbation method extends the scope of investigation by allowing for upstream and downstream development and the influence of atmospheric precursors to cyclogenesis. Finite-amplitude perturbations can include more realistic and balanced perturbations, such as upper-level PV anomalies (e.g., Wernli et al. 1999; Tan et al. 2004; Schemm et al. 2013); surface temperature perturbations (Wernli et al. 1999); unbalanced perturbations, such as numerical noise (Yamazaki and Peltier 2001); pressure perturbations (Kunz et al. 2009); or barotropic temperature perturbations (e.g., Schultz and Zhang 2007; Polvani and Esler 2007; Boutle et al. 2010; Baker et al. 2013).

The suitability of each approach depends on the problem of interest. Hence, it is usually the experimental design that determines the preferred choice of initialization method, both for the basic state and initiation of cyclogenesis. One of the main features of idealized experiments is the ability to perform controlled sensitivity experiments. Therefore, determining the appropriate approach often depends on the desired level of control of parameters used to define the initial state. In this paper we focus on an initialization method for the basic state in which the atmospheric fields are derived after defining the zonal velocity field, using height as the vertical coordinate. A similar framework has been used in several studies (e.g., Beare 2007; Wernli et al. 1999) and has the advantage that the initial wind field is easily defined analytically, allowing a direct control of the horizontal and vertical wind shear.

One way to derive the pertinent atmospheric fields, given a geostrophic zonal wind field , is to use thermal wind balance,

with *f* representing the Coriolis parameter, *g* being Earth’s gravitational acceleration, and *T* and *θ* representing the temperature and potential temperature, respectively. Scale analysis indicates that the first term on the right-hand side (rhs) tends to be small, although it becomes significant in regions of jet maxima or sharp temperature transitions in the vertical, such as the tropopause. In a Boussinesq framework the first term on the rhs can be neglected, and Eq. (1) reduces to

where refers to the background potential temperature. This equation is more straightforward to solve. However, initialization techniques adopting this approach while utilizing a non-Boussinesq numerical model introduce imbalances in the initial conditions. These imbalances subsequently initiate high-frequency variance in the model simulation. The resulting spurious inertia–gravity waves (IGWs) are arguably unphysical and potentially modify the disturbance development and are thus undesired.

In the following, we present our computationally inexpensive and mathematically straightforward solution to construct geostrophically and hydrostatically balanced initial conditions for a baroclinic channel setup. The equations are on a Cartesian grid with height as the vertical coordinate. The wind field can be freely defined by the user and thus allows for direct control of horizontal and vertical wind shear. Our method utilizes the full thermal wind Eq. (1), leading to an improved balance and hence a significant reduction of IGW activity. Furthermore, the choice of an arbitrary wind field with appropriate boundary conditions allows us to study a wide range of meteorological phenomena.

## 2. Method

First, we describe the general solution for a dry atmosphere, followed by the solution including moisture.

Given geostrophic,

and hydrostatic balance,

the pressure field can be expressed as

with , *p*, and *ρ* representing the geostrophic wind speed, pressure, and density, respectively, where both the Coriolis parameter *f* and the gravity parameter *g* are assumed to be constant. Given a wind field , Eq. (5) can be solved for pressure, provided appropriate boundary conditions in *y* and *z*. We choose one lateral boundary condition as a user-defined vertical pressure profile, which determines the stratification of the solution at this boundary. A possible realization for this boundary conditions is described in section 3b.

The simplest choice for the other boundaries is to define a vanishing wind field at these boundaries (i.e., ), implying at the upper and lower boundaries. Another choice is to allow for nonzero at the vertical boundary. For this case, we define the density at the boundary by linear extrapolation: for example, for the lower boundary , where indicates the lower boundary itself and the other indices indicate the numerical layers above the boundary. The lower boundary condition, using finite differences for the geostrophic and hydrostatic balance and the aforementioned extrapolation, can then be expressed in terms of pressure,

where the index *i* indicates the grid position in the *y* direction. Note that a similar approach can be used at the upper boundary to allow for a nonzero wind field at the top of the domain.

Given these boundary conditions, we use Gaussian elimination to solve Eq. (5), discretized with centered differences. The resulting pressure field is used to calculate density using hydrostatic balance Eq. (4), where the previously described extrapolation at the lower/upper boundary has to be applied for the case of nonzero wind at the lower/upper boundary. Given density and pressure, the temperature is calculated from the equation of state,

where *T* represents the temperature and is the gas constant of dry air.

After this procedure, the outlined method provides balanced atmospheric fields on a Cartesian grid, given a single vertical pressure profile at one of the lateral boundaries and an arbitrary wind field that vanishes at the lateral boundaries.

Moisture is included by considering the virtual temperature , where the equation of state becomes

with *q* representing the water vapor mixing ratio. We replace Eq. (7) with Eq. (8) to obtain the balanced virtual temperature field and account for the added moisture by modifying *T*. This ensures that the unsaturated moist Brunt–Väisälä frequency,

remains unchanged for added moisture. Thus, changes in moisture content do not change the static stability of the initial conditions. A possible procedure for defining the mixing ratio is described in section 3b.

## 3. Illustrative implementation of a baroclinic test case

In this section, we apply the outlined method, describe a possible realization of initial conditions, and introduce the performed simulations. The example consists of a relative narrow baroclinic zone located in a typical high-latitude environment.

### a. Model setup

The Weather Research and Forecasting (WRF) Model (Skamarock et al. 2008) is utilized to illustrate the usage of our method. The idealized model configuration consists of a zonally periodic channel, with a computational domain size of 5000 km × 2000 km. The model equations are solved on a staggered Cartesian Arakawa C grid with a grid spacing of 20 km. We employ 61 vertically stretched *η* levels with a model top at 25 km. The model is integrated forward in time using the third-order Runge–Kutta scheme with a time step of 120 s. We use Rayleigh damping at the upper 4 km of the domain and a free-slip lower boundary. Subgrid-scale parameterizations are limited to the Purdue–Lin scheme (Chen and Sun 2002) for microphysics in simulations including moisture.

Conversion from height coordinates of the initial conditions in *η* coordinates is performed using a modified version of the original WRF idealized initialization routine. Output is collected every 1 h for the unperturbed simulations and every 3 h for the perturbed simulations.

### b. Basic state

The configuration of the atmospheric fields is determined by defining the geostrophic wind field , the water vapor mixing ratio , and a vertical pressure profile at (i.e., the southern boundary), . Below is an outline of a baroclinic setup to define these parameters with the wind field:

This definition results in a zonally symmetric jet, where in the center of the domain. The width of the jet is denoted by and is the maximum wind speed of the jet located at . The vertical extent of the jet is described by and , representing the lower and upper bound of the jet, respectively. Assigning negative values to results in nonzero wind speeds at the surface. Similarly, defining above the model domain results in nonzero wind speeds at the top of the domain. The parameters *t* and *s* control the shape of the vertical wind profile below and above the maximum wind speed, respectively.

To solve Eq. (5) for , we define a vertical pressure profile at the southern boundary, require at the northern boundary and at the top of the domain, and use the option outlined in section 2 with density extrapolation at the lower boundary.

As determines the stratification at , we actually choose an approach where we derive from a prescribed static stability profile at this boundary. Given the values for tropospheric and stratospheric Brunt–Väisälä frequencies, and , the potential temperature profile is estimated with

where *j* is the index for the *z* levels, if , and if , where corresponds to the tropopause height. The lower boundary condition has to be given (see Table 1). Combining hydrostatic balance Eq. (4), the equation of state Eq. (7), and the definition of potential temperature, we can estimate the corresponding pressure,

with and where hPa and represent a reference pressure and the heat capacity of dry air at constant pressure, respectively. By selecting a surface pressure (see Table 1) we obtain a vertical pressure profile .

To determine the water vapor mixing ratio, we define a horizontally uniform relative humidity profile,

where represents the surface relative humidity, is the scale height of the relative humidity profile, and the shape of vertical decay of moisture is regulated by *n* (see Table 1). The water vapor mixing ratio corresponding to the relative humidity profile is calculated using

where is the gas constant for water vapor and is the saturation vapor pressure. Above and below the freezing level K, is estimated by

### c. Perturbation

We add a circular surface-based warm-cored temperature perturbation for the perturbed simulations,

where , K, km, and km are the amplitude, radius, and depth of the perturbation, respectively. The perturbation is placed at the surface below the center of the jet, at km and km.

### d. Simulations

Previous initialization techniques adopting a similar framework (i.e., defining the wind field followed by derivation of the pertinent atmospheric fields using height as the vertical coordinate) often utilized a Boussinesq (BQ) approach (e.g., Adakudlu 2011; Beare 2007; Wernli et al. 1999). To illustrate the reduced IGW activity and the effect on cyclogenesis for initialization by the method outlined in section 2 [referred to as the thermal wind approach (TW)], we perform three sets of simulations which contrast the BQ and TW initialization. The first set of simulations addresses the intensity and frequency of IGW activity and therefore only includes the basic state. Hence, no cyclogenesis takes place in these simulations. The other simulations illustrate the effects of IGW activity on cyclogenesis. These simulations include both the basic state and a thermal perturbation to initiate cyclogenesis. For both initialization approaches, we perform a simulation with and without moisture. A summary of the performed simulations and their respective names is given in Table 2.

Initial conditions for the TW simulations are generated using the method outlined in section 3b with the parameters given in Table 1. The method utilized to generate initial conditions for the BQ simulations is included in the appendix. The zonal mean state consists of a baroclinic jet in thermal wind balance on an *f* plane located at 70°N. Typical high-latitude conditions are mimicked by choosing a relatively low tropopause ( km) and a low surface temperature ( K). Figure 1 shows the vertical cross section of the resulting initial zonal wind, potential temperature, and water vapor mixing ratio for the mTWp simulation.

## 4. Results

The main aim of this section is to compare the different initialization methods in terms of IGW activity and the effects on perturbation growth. Thus, we do not provide a detailed analysis of the dynamical evolution in the perturbed simulations but rather highlight features related to IGW activity.

One of the advantages of our initialization method is the geostrophically and hydrostatically balanced initial fields, resulting in minimal initiation of IGW by the initialization method. Characteristics of IGW within the domain are projected in fluctuations in the sea level pressure. Figure 2 displays the sea level pressure deviation from the time mean for both initialization methods. Although both simulations feature fluctuations, the fluctuations for BQ are several orders of magnitude larger. The total wind speed at the height of the jet maximum (Fig. 3) provides an indication of the magnitude of the wind speed fluctuations initiated by spurious IGWs. The maximum wind speed of the jet is 20 m s^{−1} and for BQ the wind speed fluctuations are approximately 2 m s^{−1} with a period of about 8 h. In contrast, fluctuations in the wind speed are nearly absent for the simulation initiated with our proposed method. Note that the fluctuations are not filtered out by the model; they persist throughout the entire simulation. Hence, the basic state generated with our proposed method is less prone to produce significant imbalance and remains balanced during the time of integration. Therefore, we can conclude that any development within the domain will only result from a perturbation added to the basic state.

To illustrate that the described setup is indeed baroclinically unstable for small perturbations, we briefly describe the results for the moist perturbed simulation, where a weak initial perturbation initiates baroclinic growth resulting in cyclogenesis. Similar to previous idealized baroclinic channel simulations perturbed with a finite-amplitude perturbation (e.g., Wernli et al. 1999; Schemm et al. 2013), the main cyclogenesis event is accompanied by weaker cyclonic and anticyclonic disturbances both upstream and downstream of the main cyclone (not shown). Figure 4 shows a snapshot of the development of the main cyclone after 3, 4, and 5 days, displaying the sea level pressure deviation from the zonal mean and equivalent potential temperature at 950 hPa. After 3 days, the development of a warm front in the northeastern quadrant and a cold front in the southwestern quadrant are visible. Both features strengthen during the following 2 days, resulting in a mature cyclone after approximately 5 days. Structurally, the cyclone exhibits similarity with the previously mentioned simulations in terms of the relative location of the frontal structures and the development of a warm air seclusion. The major differences compared to previous simulations are the time and length scales of cyclogenesis, which are both considerable smaller in our simulations. These differences can be attributed to the nature of the setup, as the configuration mimics typical high-latitude settings, featuring a larger Coriolis parameter and a lower tropopause height, leading to smaller wavelengths and larger growth rates (e.g., Vallis 2006). The example of cyclogenesis in the perturbed simulation indicates that the initialization method can be successfully used to simulate rapid cyclogenesis with minimal spinup requirements because of a significant reduction of spurious IGW activity. For further demonstration of the setup see Terpstra et al. (2015).

By performing perturbed simulations with both initialization methods, it is possible to quantify the effects of these spurious IGWs on developing disturbances. The total vertically integrated eddy kinetic energy (Ke) is used as a proxy for the intensity of the disturbances within the domain. For simulations without moisture, Ke features small differences between both initialization methods, whereas for the simulations including moisture these differences are more prominent (Fig. 5). In both cases simulations initiated with the TW method exhibit larger values of Ke. Hence, spurious IGWs within the domain weaken the cyclone intensity and this effect becomes larger when moisture is included. As moist cyclogenesis features more nonlinear processes than the dry version, it is not surprising that the effects of IGWs are stronger for simulations including moisture. Hence, these experiments illustrate that IGW activity due to imbalances in the initial state can result in a modification of the developing cyclone. As these IGWs are mainly introduced because of a noncomplete version of the thermal wind equation, we argue that they are undesirable and should therefore be minimized by the initialization method.

## 5. Summary and conclusions

In this paper, we outlined a method to obtain hydrostatically and geostrophically balanced atmospheric fields suitable for initialization of idealized numerical simulations. A feature of this method is the possibility to freely define a wind field and moisture distribution in height coordinates, whereafter balanced atmospheric variables (pressure, density, and temperature) are derived. While this feature has been the preferred method of initialization for several studies, we introduce a new method to arrive at the initial fields, which shows a considerable reduction of spurious IGW activity.

Similar improvements in balance might be possible by utilizing other methods, such as a dynamical initialization method, where the model is integrated backward and forward until a satisfying atmospheric balance is achieved (e.g., Lynch and Huang 1994). Another option is applying filtering to remove spurious IGWs. However, a major drawback of these filters is that they are not selective and hence modify both the unwanted IGWs and the developing feature of interest. Although both methods are able to enhance the balance of the initial conditions, they are somehow detrimental to the key idea behind idealized simulations: that is, the possibility to perform experiments with direct control over the atmospheric environment. It is therefore preferable to utilize an initialization method that is able to reduce unwanted IGW activity.

We performed both perturbed and unperturbed numerical simulations to illustrate the applicability and balance of the basic state utilizing a baroclinic setup constructed with our outlined method. Furthermore, we contrasted our proposed method with an initialization utilizing a previously used Boussinesq framework. Unperturbed simulations initiated with our proposed method feature only minimal adjustments of the atmospheric fields over a 5-day period, indicating that imbalances in the initial conditions are nearly absent, while the degree of imbalance for the Boussinesq approach is several orders of magnitude larger. The imbalance yields a reduced cyclone intensity for perturbed experiments. The latter effect is enhanced when moisture is included.

Another advantage of our method is that the basic state wind field can be defined directly, thus allowing for studies of the influence of both horizontal and vertical wind shear on cyclogenesis. One drawback of our outlined method is the restricted ability to control the inner domain temperature and static stability. The main control over the static stability is at one of the lateral boundaries, whereas the interior static stability is a result of the definition of the wind profile together with the boundary conditions. However, this restriction in controlling the interior static stability can be overcome by defining appropriate boundary conditions to obtain the desired static stability over the entire domain.

The baroclinic example provided in section 3 is only one possibility to define specific initial conditions. The solution method presented in section 2 is more general and allows for a large degree of freedom in defining the wind field and boundary conditions. This flexibility aids the potential applicability of this initialization method to a wide range of research questions, such as Rossby wave breaking, shallow baroclinic zones, reverse shear conditions for cyclogenesis, barotropic instability, embedded low-level jets, etc. The method is especially suitable for studying phenomena that include multiple jets: for example, low-level, tropospheric, and stratospheric jets. As the initial conditions are balanced, the initialization technique is particularly adequate for studying mesoscale features with relatively short time scales. Furthermore, the perturbation method is independent of the basic state, allowing for any type of perturbation to be added, depending on the phenomenon of interest.

## Acknowledgments

This work is supported by the Norwegian Research Council as part of the project High Impact Weather in the Arctic (Project 207875).

### APPENDIX

#### Implementation Initial Conditions Applying Boussinesq Approximation

Following the implementation of the lateral boundary condition for in section 3b, the background potential temperature and pressure are defined by the Brunt–Väisälä frequency [see Eqs. (11) and (12)]. The background density is derived from the equation of state together with the definition of the potential temperature . Potential temperature and pressure deviations from the background state are derived by horizontally integrating the thermal wind Eq. (2),

and the geostrophic equation,

assuming zero deviation in and at . The density deviation is obtained by vertically integrating the hydrostatic equation,

The total fields (i.e., , etc.) are then passed on to WRF for initialization of the baroclinic channel.

## REFERENCES

**138,**1297–1307, doi:.

**140,**96–110, doi:.

**86,**1609–1614, doi:.

*Tellus,*

**46A,**583–597, doi:.

*J. Geophys. Res.,*

**112,**D23102, doi:.

**66,**1569–1592, doi:.

*Quart. J. Roy. Meteor. Soc.,*doi:, in press.

*Atmospheric and Oceanic Fluid Dynamics.*Cambridge University Press, 768 pp.

*J. Geophys. Res.,*

**116,**D05108, doi:.