Part 2 - Find fronts

Prepared by Robin Noyelle.

The objectives of this part are the following:

compute the gradient of \(\theta_e\)
create a mask for frontal structures
separate between cold and warm fronts

Preparation

Create a new notebook for this part, select the correct kernel and give a name to the first cell using Markdown format
Add a new cell and import the required packages as in the previous part
Load data as previously, this time also load the data for meridional and zonal wind speeds

Compute the gradient of \(\theta_e\)

We saw in the previous part that the computation of \(\theta_e\) allows to visually separate the air masses: warm air masses have a high \(\theta_e\) while cold air masses have a low \(\theta_e\). As a consequence, the fronts, i.e. the regions at the interface between two such air masses are characterized by a high value of the gradient of \(\theta_e\): in other words, there is a large change in \(\theta_e\) where fronts occur. We will use this observation now to systematically and objectively identify the frontal structures.

The first step is to compute the gradient of \(\theta_e\). On a \((x,y)\) plane this gradient is defined as:

\[ \nabla \theta_e = \frac{\partial \theta_e}{\partial x}\vec{e_x} + \frac{\partial \theta_e}{\partial y}\vec{e_y} \]

where \(\vec{e_x}\) and \(\vec{e_y}\) are the unit vector in the \(x\) and \(y\) directions. Can you recall what is the intuitive meaning of the gradient operation? In our particular case, what would be the meaning of the norm and the direction of the vector \( \nabla \theta_e \)?

The L2 norm of the gradient is:

\[ || \nabla \theta_e ||_2 = \sqrt{ (\frac{\partial \theta_e}{\partial x})^2 + (\frac{\partial \theta_e}{\partial y})^2 } \]

The L2 norm of a vector is just its length. Thererefore, the higher the L2 norm of the gradient of \(\theta_e\), the more \(\theta_e\) changes rapidly on the horizontal plane: a high value of this norm therefore can be the mark of a front !

As a first step and for illustration, we will treat the lat-lon coordinates of our fields as if they were \((x,y)\) coordinates (which they are not).

a. Write a function which takes a lat-lon field as input and returns its gradient in the lat and lon directions. You may be interested by the differentiate method of xarray.

b. Write a function that compute the L2 norm of a vector on the plane.

c. Write a third function, making use of the two previous ones to compute the norm of the horizontal gradient of a lat-lon field. Test the function on the \(\theta_e\) field computed for the first day of our data. Plot the resulting norm of the gradient and describe what you see. Can you identify the frontal structures ?
As you may know, we do not live on a flat planet. This has important consequences when computing the horizontal gradient of \(\theta_e\): can you see what is the issue of assuming that the lat-lon coordinates correspond to a flat 2D plane? In the following we will use the following function to compute the horizontal gradient on a sphere:
```
def compute_gradient(field):
    R_earth = 6381 * 1e3

    _, LON = np.meshgrid(field.lon, field.lat)
    deg2rad = 180 / np.pi

    grad_lon = field.differentiate("lon") * deg2rad / (R_earth * np.cos(LON * deg2rad))
    grad_lat = field.differentiate("lat") * deg2rad / (R_earth)
    return grad_lon, grad_lat
```
What are the main differences with your function to compute the gradient? Compute the norm of \(\nabla \theta_e\) field for the first day of our data with this new function. Compare with the result you obtained at the previous question. What is the unit of \(\nabla \theta_e\)?

Identify frontal structures

Frontal structures can be identified as the interface between two different air masses. This would then correspond to a high value of the norm of \(\nabla \theta_e\).

We will now compute a mask to identify frontal structures. The idea is to have a lat-lon xarray which, at each grid point, has a value 1 if there is a frontal structure and 0 if there is no frontal structure. To do so, we will apply a filter depending on the size of \( || \nabla \theta_e ||_2 \): if \( || \nabla \theta_e ||_2 \) is bigger than a user-chosen value threshold_gradient then there is a frontal structure. We first take threshold_gradient=2.5 K/(100km). Compute this mask and make a plot of it. Are the frontal structures you find coherent with what you would expect from a visual inspection ?
Do the same using, threshold_gradient=2.0 K/(100km) and threshold_gradient=3.0 K/(100km). What differences do you observe ?

A lot of the structures that are found with this procedure are spurious fronts, especially in the tropics. Additional constraints apply to the definition of a front: it must be long enough for example and therefore isolated grid points cannot constitute fronts. There are ways to get rid of them but here we will not do that and we will assume that they are indeed fronts. In the following we will use threshold_gradient=2.5 K/(100km) by default.

Separate between cold and warm fronts

As explained in the introduction, fronts usually come in two forms: cold fronts and warm fronts. The distinction between the two is straightforward: a cold front occurs when cold air replaces warm air, whereas a warm front occurs when warm air replaces cold air. In this last section of the second part we will let you search a little bit for how to do the cold vs warm fronts separation. You may find interesting to draw a plot by hand for a warm and a cold front to understand what is happening. Here are some intermediate questions that may help you implement this in python:

At one grid point, how can you know in which direction the air is flowing?
At one grid point that you identified as a front, how do you know in which direction is the warm air and in which direction is the cold air?
How can you compute whether two vectors are pointing in the same direction or in two different directions?
In python, if you have two masks mask1 and mask2 representing arrays with True and False values, you can compute their intersection (i.e. determine where both masks are True) as mask1 & mask2

Once you have implemented a function to separate between cold and warm fronts, make a nice figure of warm and cold fronts for the first day of our data set. You can try to overlay the two with the \(\theta_e\) field for example by using the .plot.contour() method for your fronts data and putting red contours for warm fronts and blue contours for cold fronts. Save this figure as a png file. Do the same for one other date per season for the first year of the data set. Comment what you see in the figures.

Bonus question: when detecting fronts for the four previous dates you may have noticed that some regions seem to be always associated with fronts (which regions for example?). Those are called stationary fronts and they are usually not of major interest for dynamical studies of the climate model. You can filter them out by imposing that at your frontal structure the wind should be big enough, typically larger than 5 m/s. Try to implement this third category of fronts and redo your previous figures.