Choosing between the two LES definitions of sub-grid fluxes

[dhruvbalwada]

This to catch up @rpa, @paigem and @TomNicholas, on the hole that @jbusecke and I are just emerging out of (I think successfully). I apologize that this is very long issue, and some of the details might not be fully present here and can be discussed in a zoom call.

Context: This issue arose out of the discussions based on #50, where we were looking for a clean way to decompose the sub-grid flux contributions ( $Q^{*} = \overline{Q_H} - Q_L$, where $\overline{Q_H}$ is the smoothed flux computed from high resolution fields and $Q_L$ is the flux computed from the low-resolution fields.) that come from the atmosphere vs the ocean.

What we tried and the problem with it: We tried to establish the contributions from just one component's small scale. E.g. for establishing the contributions from oceanic small scale, we computed the flux using all the full-resolution fields and remove the flux that was computed using the full-resolution atmosphere and smoothed ocean ($Q^C_{Ocean}$ = $\overline{QH} - \overline{Q{LOcean}}$ , where $\overline{Q{LOcean}}$ has been smoothed since it would have explicit small scale structure from atmosphere). The opposite was done to establish contributions due to atmosphere $Q^C{Atmosphere}$. The problem that we came up against was that the $Q^ = Q^C{Ocean} + Q^C{Atmosphere} + Q^_{Residual}$, where the residual term was very large and comparable to the other terms (also it was not completely obvious what this contribution was. [*I apologize for the funky notation switching between and C (should have been ), but that is because github has terrible math rendering.]

Solution: In the process of trying to figure out what was going on, we discovered that there are two complementary and equivalent ways to define the sub-grid contributions. The first is the definition we were using - $Q^{*} = \overline{Q_H} - Q_L$, and the second is $Q^{*} = \overline{Q_H} - \overline{Q_L}$ (Sagaut 2006: LES book*). The difference between the two is whether $Q_L$ is smoothed or not. They also have slightly different-looking evolution equations (even though both evolve smoothed $T$ in exactly the same way):

$$
\begin{aligned}
\partial_t \overline{T} & = Q (\overline{u},\overline{T}) + [\overline{Q(u, T)} - Q (\overline{u},\overline{T})] \\
&= Q_L + Q^*
\end{aligned}
$$

or 
$$
\begin{aligned}
\partial_t \overline{T}  & = \overline{Q (\overline{u},\overline{T})} + [\overline{Q(u, T)} - \overline{Q (\overline{u},\overline{T})}] \\
&= \overline{Q_L} + Q^{**}
\end{aligned}
$$.

There are a number of advantages to the second formulation:

• $Q_L$ can contain structure at sub-filter scales (due to the non-linearity) but $\overline{Q_L}$ will not have this, and so in many context $\overline{Q_L}$ is the preferred "large-scale" forcing. (This emergence of smaller scale structure than what is in the input fields is a known problem, and the reason why de-aliasing approaches are used in pseudo-spectral codes).
• $Q^{**}$ leads to a much cleaner decomposition, which can be physically interpreted more easily. (This point is described below.)

Here are two snapshots of the different formulations. Of course there are visual differences between the two small scale contributions (3rd column), but they are equivalent since they both lead to the exact same temperature evolution. This is because what we call large scale has also changed between the two (notice that column 2 in bottom row smoother).

Decompositon:

Here we formulate the discussion in terms of an advective non-linearity (AB), but the exact same reasoning follows for any non-linear function. As we checked for some arbitrary functions, and seems to be holding also for much more complex bulk-formulae as shown in figures below.

For an AB type non-linearity we have: $$ Q^* = \overline{A.B} - \overline{A} .\overline{B} $$

This can be further decomposed using $A = \overline{A} + A'$ as $$ Q^* = [\overline{ \overline{A}.\overline{B}} - \overline{A} .\overline{B}] + [\overline{A'.\overline{B}} + \overline{\overline{A}.B'}] + \overline{A'.B'} $$

where the bracketed terms are called the Leonard stress, cross-stress, and Reynold's stress respectively.

Note that for Reynold's operators, the Leonard stress and cross-stress are zero for AB type non-linearities, and for Reynold's operators only Leonard stress are zero for more complex non-linearities (can check easily be working through $A^2B$ for example). As the next immediate step, @jbusecke is planning to test the above idea by using coarsened fields to compte Q* and Q**. They should be similar (exactly same?) if the above remark is true for any non-linearity.

But also note that the above expression can be written as

$$
Q^* = [\overline{ \overline{A}.\overline{B}}  - \overline{A} .\overline{B}] +  Q^{**}
$$

as

$$
Q^* = \overline{Q_H} - Q_L = \overline{Q_H} - \overline{Q_L} + \overline{Q_L} - Q_L = Q^{**} + (\overline{Q_L} - Q_L).
$$ 

The difference between $Q^*$ and $Q$ is the Leonard term, which is purely a result of large-scale fields - it is the small-scale structure that results when two large-scale fields interact in a non-linear way. On the other hand there is a nice decomposition for
$$ Q^{} = [\overline{A'.\overline{B}} + \overline{\overline{A}.B'}] + \overline{A'.B'} $$ which has contributions from (1) how the small scale A interacts with large scale B, (2) how large scale A interacts with small scale B and (3) how the small scale A interacts with small scale B. We can interpret A and B as the ocean and atmopshere parts, so contribution (3) seems to the focus of most of the air-sea interaction literature, which is focussed on analyzing how wind anomalies form over mesoscales SST anomalies.

This decomposition is shown in figures below. We find that the cross-terms have large contributions, while the Reynold's term is small (Does this mean that the part the air-sea interaction literature has been focussed on is a second order effect in the air-sea flux?).

Single timestep terms of the $Q^{**} Decomposition$

5 year average terms of the $Q^{**} Decomposition$

An extended set of notes on this topic (as our thinking evolved) can be found here

[jbusecke]

I modified the coarsen interval for our comparison to 50 (matching the filterscale) and got this very encouraging result:

Global mean small scale contribution for coarsening (Note that for a reynolds operator $Q^{*} = Q^$) and smoothing.

To me this means we can confidently say that our results are fairly insensitive to the method chosen 🚀. This also gives an interesting preview of the scale dependence of the results: With a coarsening scale of 30 we got ~half of the small scale contribution. To be explored in a follow up study.

[jbusecke]

I will go ahead and compute the $Q^{**}$ terms for the full time now.

[dhruvbalwada]

Looking good. For completeness also add smoothed (Q*) to above plot.

[jbusecke]

Before I do that I actually took a quick look at the spatial comparison between the methods:

Top-Left: Results calculated purely from coarsening; Top-Right: Results calculated purely from filtering; Bottom-Left: Same as Top-Right but coarsened. Bottom-Right: Difference between Top-Left and Bottom-Left.

We can see the same negative features in both methods ✅ but the filtering has some prominent positive anomalies that do not appear in the coarsened output. Some of these I interpret as small scale noise (e.g. the smaller patches in the SO), but the effect of the coasts and the equator seems to persist.

I do not think this is a blocker, I am just wondering how to interpret these. Are these effects of very strong local gradients which influence the filtering differently than the coarsening?

The differences show this large scale alternating patterns, which seems to mostly average out in the global average, but is not negligible locally!

Thoughts on how to deal with this?

Keep in mind this is only computed on a very short timeframe (~1year)

[jbusecke]

This also made me notice that my coarsening does not treat the coordinates properly. It replaces the values of e.g. xt_ocean (-270-90) with integer indices.

EDIT: Here is a partially fixed version for better visual comparison. #56 will be applied to all data in the next full processing round.

[jbusecke]

Here is a longer version of the timeseries above:

This clearly shows that interannual variability is minor (and I agree with @rabernat now that we should not show this full timeseries). But I would argue that the seasonal variability is significant?

How do folks feel about changing the 20 yr timeseries to a seasonal climatology?

[rabernat]

How do folks feel about changing the 20 yr timeseries to a seasonal climatology?

:+1:

[jbusecke]

Here we go:

The line is the ecmwf algo and the shading represents the range of 5 algos (~20% variation for Latent and ~35% for sensible).