This post takes a look at the default behavior of ICAR with regards to the calculation of the Brunt-Väisälä frequency N and discusses an alternative method that is arguably better suited from a theoretical point of view. Note that the data and results here are not peer reviewed. However we have a manuscript in preparation that discusses aspects of this. I’ll add the corresponding links once they are available.

The underlying theory

Linear mountain-wave theory assumes that the state of the atmosphere may be described by a background state and perturbations to it. In case of mountain waves the perturbations are induced by the topography and their properties may be studied by analyzing the Taylor-Goldstein equation. It is derived from the linearized and Boussinesq approximated Navier-Stokes equations.

ICAR is based on linear mountain-wave theory, more specifically on the set of equations derived by Barstad and Grønås (2006). The equation relevant for this post is the buoyancy equation (A5):

    \[ \frac{\partial b}{\partial t} + \vec{U}\cdot \nabla b + N^2w = 0 \]

Note that here N is treated essentially as a constant, independent of the horizontal and vertical coordinates. While not explicitly stated, this already indicates that N is a property of the atmospheric background state. In this case, neglecting the influence of moisture, it is given by

    \[ N^2 = -\frac{g}{\theta_0}\frac{d\theta_0}{dz} \]

where \theta_0 is the horizontally homogeneous potential temperature profile of the atmospheric background state.

If that’s not enough to convince you, this is stated more explicitly in other works that go into the details of the derivation of the Taylor-Goldstein equation. See, for instance, Nappo (2012) on page 32. Here N is derived from the density profile of the background state (the referenced equation 1.67):

    \[ N^2 = -\frac{g}{\rho_0}\frac{\partial \rho_0}{\partial z} \]

This is where potential problems arise. ICAR uses the potential temperature profile and moisture distribution of the perturbed state of the atmosphere to calculate N. Put in an equation it would look like this (neglecting moisture once again):

    \[ N^2 = -\frac{g}{\theta_0+\theta_1}\frac{d(\theta_0+\theta_1)}{dz} \]

Here \theta_0 is the potential temperature profile of the background state and \theta_1 the perturbation to it. Note that \theta_0 = \theta_0(z) and \theta_1 = \theta_1(x,y,z). As illustrated there’s no basis in the underlying theory to use the above equation.

We’re essentially leaving the safe harbor of linear mountain wave theory. Even if you were to run an ICAR simulation with a perfectly linear background state you lose an important part of the justification of why you’d think that your results are correct. As shown above, the standard procedure of ICAR to calculate N violates an assumption of linear theory.

Another argument in defense of calculating N from the state of the atmosphere in ICAR might be that it must be a better approximation of what is actually going. The ICAR domain has a finer grid spacing than the forcing data set and, therefore, more detailed information about the moisture and temperature distribution. However, that this information is correct has yet to be shown conclusively. In fact I think it’s even harder to argue since all these distributions result from a wind field that was calculated from an essentially invalid underlying assumption. Why would that result in more realism?

The difference between both approaches

Lets have a look at how large the difference between the two methods actually is! We’re going to take a look at a Northwest-Southeast cross section through the South Island of New Zealand during a comparatively large precipitation event. During the two day period we’re investigating, north-westerlies prevail and ICAR predicts large amounts of precipitation to fall out along the upwind slopes of the topography. The simulations are based on Horak (2019), however, I’ve rerun them with a model top at approximately 15\,\text{km} above topography.

All following panels show the northwest-southeast cross section (oriented parallel to the background flow). Northwest is on the left side, Southeast on the right. The top row shows the vertical velocity w and Brunt-Väisälä frequency when N is calculated from the current state of the ICAR domain, while in the bottom row N is calculated from the forcing data set (the method compatible with linear theory). All plots were done for the 4th of May 2015, 06:00 UTC.

Even from a purely qualitative point of view the differences are obvious. While the basic structure of the w-field is similar, the magnitude and shape of the up- and downdraft regions clearly has changed. This is not suprising since the Brunt-Väisälä frequencies used to calculate the wind fields differ as well (see panel b and d).

My explanation of the differences is that if N is calculated from the state of the ICAR domain, there’s something I’d loosely refer to as feedback mechanism happening. First, N is determined by the potential temperature field and the moisture distribution in the ICAR domain. In a next step N determines the structure of the wind fiel. Then, all scalar atmospheric quantities, such as temperature and moisture, are advected in this wind field. And from these modified distributions the next N field is calculated which in turn modies the wind field and so on.

If, however, N is calculated from the forcing data, the advection of scalar quantities in the wind field does not feed back into the calculation of N. While the distribution of these quantities in the ICAR domain changes, all input data relevant for N is taken from the forcing, which is what ICAR assumes to be the background state.

Conclusion

While we’ve seen that linear mountain-wave theory states that N is a property of the background state, we’ve not shown that ICAR results are indeed “better” (whatever that means) if the calculation method of N is adapted to satisfy this requirement imposed by theory. We’ve only shown that they are different. I do indeed claim that the representation of the wind field is improved in a manuscript that’s currently in preparation, but that’s not out yet (I’ll update this post once it is though).

ICAR is interesting for downscaling because it is based in laws of nature. Adhering to these laws and their assumptions is important – otherwise any result obtained with the model is automatically less reliable due to a crack in the foundations. It doesn’t really matter whether the model so far was able to correctly model measurements. It is possible to draw correct conclusions from the wrong assumptions. A comparison to measurements alone is not enough to conclusively validate a model (e.g. Schlünzen, 1997). But once we find something in a model that contradicts the theory we should “fix” it – which is what we did.

Again, just to make it clear – this has not (yet) been peer reviewed. I’ll update this post accordingly in the future.

Modifications to ICAR

At the time of writing the calculation of N from the forcing data set is available as an option in the development branch of ICAR. Additionally, more recently, I’ve added another option that allows for the specification of a variable in the forcing data set from which the Brunt-Väisälä frequency should be read directly.

Update 23/01/2020: All these modifications have found their way into the master branch of ICAR.

If you want ICAR to calculate N from the forcing data set, set the following flag in the linear theory options (lt_parameters namelist):

N_from_forcing = true

If you want to specify N yourself, you need to additionally add the field where you do that to the variable names (var_list namelist):

nvar = "N2ICAR"

In the example above the forcing data set contains a variable called “N2ICAR” that contains the square of the Brunt-Väisälä frequency in every grid cell at every time step. Note that it needs to be the square of N. Currently ICAR reads in those values and checks whether they are within the boundaries specified by min_stability and max_stability. If not it adjusts the read-in data accordingly

It goes without saying that all of these are fairly recent adaptions of ICAR that need more testing. While both seemed to work for test cases I ran, there’s definitely more that is, so far, untested.

References

Barstad, I. & Grønås, S., Dynamical structures for southwesterly airflow over southern Norway: the role of dissipation, Tellus A, Munksgaard International Publishers, 2006, 58, 2-18

Nappo, C. J., Nappo, C. J. (Ed.), The Linear Theory, An Introduction to Atmospheric Gravity Waves, Academic Press, 2012, 102, 29 – 56

Horak, J., Hofer, M., Maussion, F., Gutmann, E., Gohm, A., and Rotach, M. W.: Assessing the added value of the Intermediate Complexity Atmospheric Research (ICAR) model for precipitation in complex topography, Hydrol. Earth Syst. Sci., 23, 2715–2734, 2019.

Schlünzen, K.: On the validation of high-resolution atmospheric mesoscale models, Journal of Wind Engineering and Industrial Aerodynamics, 67-68, 479 – 492, Computational Wind Engineering, 1997.