# Note on air-sea flux project

## What we have

The small scale contribution from all mechanisms.
$$Q^* = \overline{Q_H} - Q_L$$

where $\overline{(.)}$ is a filtering operation, $Q_H = f(a, b)$ is the heat flux computed from high resolution fields, and $Q_L = f (\overline{a}, \overline{b})$ is the heat flux computed from filtered/low resolution fields.

The goal was to obtain some $Q^*$ that might be missing from low resolution coupled simulations, and understanding how large this might be.

## How can we understand more? (what we did first)

We are trying to arrive at a meaningful decomposition of $Q^*$ into contributions from atmospheric and oceanic components of the system (for now ocean atmosphere, but you could also look at only velocities or tracers from both atomsphere and ocean).

A natural thing to look at would be the following:

Compute the air-sea flux $Q_{L,Ocean}$ from an only partially smoothed field (e.g. only smoothing ocean fields) and construct a resulting small scale term similar to $Q^*$

$Q^*_{Ocean} = \overline{Q_H} - Q_{L,Ocean}$

We can do the same for the atmosphere ($Q^*_{Atmosphere}$), and diagnose a residual term ($Q^* = Q^*_{Ocean} + Q^*_{Atmosphere} + Q^*_{residual}$).

> Single time snapshots of $Q^*$ decomposition terms

> 5 year averages of $Q^*$ decomposition terms

### Issues
- The $Q^*_{Ocean}$ and $Q^*_{Atmosphere}$ (and subsequently the residual) terms include small spatial scales. This is problematic because we are trying to decompose a smooth field.
- This decomposition has no simple and clear mathematical interpretation (*this is based on working through how large and small scale parts of the fields travel through operations. Shown in handwritten notes below for a simple AB type non-linearity.* ), and the residual includes a bunch of messy terms which are not negligible (see picture above). (*Same holds true even in long term averages*)

## What we did instead (so far)
Since $Q^*_{Ocean}$ and $Q^*_{Atmosphere}$ have small scale structures, as an ad hoc solution we started to 'simply resmooth' the mixed large scale terms and compute
$Q^{**}_{Ocean} = \overline{Q_H} - \overline{Q_{L,Ocean}}$ and $Q^{**}_{Atmosphere} = \overline{Q_H} - \overline{Q_{L,Atmosphere}}$ terms and compute the resulting small scale contribution. But this does not give a good decomposition for $Q^*$ a mathematical sense (*see handwritten notes, where no terms drop out*), and this is revealed by the fact that the residual term ($Q^* = Q^{**}_{Ocean} + Q^{**}_{Atmosphere} + Q^*_{new res}$) is quite large and we do not have a clear phyiscal interpretation of it.

> Single time snapshots of the 'wrong' $Q^*$ decomposition terms

> 5 year averages of the 'wrong' $Q^*$ decomposition terms

One way to make progress, is to categorically look at a different definition of 'small scale contribution'

$Q^{**} = \overline{Q_H} - \overline{Q_L}$

If you define your small scale contribution like this, you can actually decompose $Q^{**}$ nicely into:
$Q^{**} = Q^{**}_{Ocean} + Q^{**}_{Atmosphere} + Q^{**}_{Res}$

Firstly, at an empirical level $Q^{**}_{Res}$ is small. Also, we believe (based on simple analysis shown in handwritten notes) that we have a clear physical interpretation of $Q^{**}_{Res}$. It is the heat fluxes contributed purely due to small scale coupled effects between ocean and atmosphere, while $Q^{**}_{Ocean}$ represents the effects of large scale atmosphere interacting with small scale ocean and vice versa for $Q^{**}_{Atmosphere}$. The coupled effects encompassed in $Q^{**}_{Res}$ also seem to be most directly connected to the air-sea interactions literature, where atmospheric wind anomalies are observed to be correlated with ocean temperature/SSH anomalies.

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

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

> Note the positive values in $Q^{**}$ which do not show up in $Q^*$

### Issues
- $Q^{**}$ decomposes well, and seems to be not primarily an atmospheric signal. However, it is not what one might want to inject in a low resolution simulations (given that $Q_L$ has been smoothed in its defintion, and this diverges from the LES style definition of small scale fluxes.)
- This problem is emerging because our filters are not Reynold's operators, and filtering functions composed of only filtered fields return a large anomalous component.
- We were very surprised to find that filtered $Q_L$ looked very different from $Q_L$. **Can we make a plot of this to see what is going on here (compare $Q_L$, $\overline{Q_L}$ and the anomaly term)? Some how the non-linearity is producing smaller scale structure that is absent in the original fields.**

## Where do we go from here?
- Use an alternate (or modified interpretation of) decomposition where we say:
$$
Q^* = \overline{Q_H} - Q_{L}
$$
$$
Q^* = \overline{Q_H} - ( \overline{Q_{L}} + Q_{L}')
$$
$$
Q^* = Q^{**} - Q_{L}'
$$
where $Q^{**} = \overline{Q_H} - \overline{Q_L}$.
$$
Q^* = Q^{**}_{Ocean} + Q^{**}_{Atmosphere} + Q^{**}_{Res} - Q'_{L}
$$
- Here $Q^{**}_{Res}$ is the impact of small scale coupled components.
- Here $Q^{**}_{Ocean}$ and $Q^{**}_{Atmosphere}$ are components that explain how large scale in one fluid interacts with small scale in other fluid.
- Here $Q^{*}_{Res}$ is the hardest component to understand, and seems to be a result arising due to small scale structure that emerges in $Q$ in the presense of smooth fields, which is not present in either of the input fields.
- Also $Q^{*}_{Res}$ is strongly anti-correlated with $Q^{**}$.
- How can we understand this further:
- **Maybe looking into $Q_L'$ might help.**
- The size of $Q_L'$ should shrink as the filter gets sharper - multiple applications of the filter produces the same field.
- Is this some property of a coupled system? (Can looking at simple Hasselman model help?)

- If we would use a reynolds operator instead of a filter, we would actually be able to use the $Q^{**}$ decomposition because there $Q^*$ should be the same as $Q^{**}$.
- There are some technical challenges but Dhruv and I are convinced that its technically possible to use a box average (some form of `.coarsen` or '.groupby().mean()') to achieve this quickly. It is probably worth checking.
- If changing the filter will change the qualitative behavior of Q* that would be really weird and interesting.

> Commentary:
- JULIUS:It (using coarsening) would be kind of a bummer though because we would basically toss all of our work, but I guess that is science for ya.
- DHRUV:
- I think we learnt a lot. Specially the fact that we need to be careful about how we filter, as it can make qualitative changes.
- Also eventually we might want to combine different filters with coarse graining, to get closer to some operator that is more comparable to a low-resolution model. The goal is to do something that helps low resolution simulations, where the grid is much coarser than the grid our filtered fields live on.
- Would be curious to see if coarsened fields resemble more of Q* or Q**.

## Appendix

### A1: Hand written notes
Working through the math for a simple AB type nonlinearity. While the exact expressions would not hold true for heat fluxes, the general principles about which term is a combination of large scales or large scale+ small scale or small scale+small scale would still apply.

### A2: LES decompositions and terminology

As oceanographers we are used to thinking of eddies in terms of Reynold's operators ($\overline{c} = \overline{\overline{c}}$ and $\overline{c'} = 0$), but most filters don't follow these properties. A more general definition of the eddy or sub-grid fluxes for an advective non-linearity is
$$
\tau = \overline{A.B} - \overline{A} .\overline{B}
$$

This can be further decomposed using $A = \overline{A} + A'$ as
$$
\tau = [\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 1: For a Reynold's operator and advective non-linearity the Leonard stress and cross-stress are 0.
> Note 2: For a Reynold's operator and a general non-linearity (can check easily with something like $A^2B$) the Leonard stress is 0, but not the cross-stress or Reynold's stress. (this is the relevant case for the air-sea flux formula)
>Note 3: A sharp spectral filter approaches a Reynold's operator, as subsequent filtering does not impact any of the un-filtered scales in the first go.

What we called $Q'_L$ above is simply the negative of the Leonard stress, can be seen for the advective non-linear form as

$$
\tau = \overline{A.B} - \overline{\overline{A} .\overline{B}} + [\overline{\overline{A} .\overline{B}} - \overline{A} .\overline{B}]
$$

we have called the brackted term here as $-Q'_L$, but notice that it is simply the Leonard stress. The first two terms make up what we call $Q^{**}$.

The Leonard stress is often parameterized, in this $AB$ context, using simple models (explicit filtering, Bardina) as it is purely a function of resolved/large scale flows and can in principle be directly computed from large scale flows. Pantano and Sarkar 2001 showed how to do this some particular non-linear models (also discussed slightly more pedagogically in Saugat book).

> Some LES papers do define the sub-grid stress purely as what we call $Q^{**}$. Probably with the idea that the resolved terms should not have smaller scales than filter scale.

### A3: Two forms of LES equations

Following the discussion Sagaut 2006 or Winckelmans et al 1996, the equations for the filtered (low-resolution) temperature field can be formulated in 2 ways:

$$
\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}
$$.
Here we are considering the $Q(u,T)$ to be the surface flux, but it could just as well have been the advection term.

While both formulations are for the temporal evolution of $\overline{T}$, and equivalent, the RHS look very different. In the first form, the system is driven by the large scale flux $Q_L$ and the corresponding small scale flux contribution $Q^*$, while in the second case these terms are replaced by $\overline{Q_L}$ and $Q^{**}$.

So, one might ask how to choose between one form of the other. The choice to some degree depends on the practicality of how the problem needs to be solved.

The difference in the two forms lies in what scales can have variance/ be active in the forcing terms. Notice that in the second forumlation, all the RHS terms are filtered. This means that none of them will have contributions coming from sub-filter scales. However, this is not true for the first form, where the non-linearity inherent in the functional form of $Q$, might produce variability at or amplify scales smaller than either of the input fields (*see Appendix A4 for a toy example*).

Since a low resolution model does not have the grid points available to resolve these amplified small scales, what happens to these non-linearly amplified small scales? This question can alternatively be framed as - what does a low resolution model really get forced by: $Q_L$ or $\overline{Q_L}$? There is likely no perfect answer to this question. If there is signficant variance near the grid scale in a low-resolution model, then computation of $Q_L$ will predict small-scale structure. Since this structure can't actually be resolved, it would get aliased onto large scales. So the best one can hope for is $\overline{Q_L}$ as atleast the small scale would be well defined, but even that requires some appropriate reconstruction scheme (e.g. dealiasing in spectral codes). However, this choice is to some degree dependent on the model developers.

A small practical advantage of the second formulation is that one does not need to consider describing or understanding the Leonard terms ($Q_L - \overline{Q_L}$) (since they are just baked into the definition of large scale).

> Note: it is not sufficient to simply look at long term structure of only $Q^*$ or $Q^{**}$, but rather it should always be considered in combination with the appropriate large scale flux. *Think about this when considering the equatorial signal.* **In the figures below we can look at framework of the top row or bottom row, the difference between the two comes from the fact that in the second row the second column is filtered. **


### A4: Impact of non-linearities on scales

From https://www.astro.auth.gr/~vlahos/GravitoplasmaWS1/pseudo-spectral_2.pdf in the section of aliasing and de-aliasing.

Consider a nonlinearity of the form $AB$.
Let $A = a\sin(k_ax)$ and $B = b\sin(k_bx)$ be two signals, considered here for illustration purpose. Then the non-linear function gives us

$$
\begin{aligned}
AB & = a\sin(k_ax)b\sin(k_bx) \\
& = \frac{1}{2}ab [cos ((k_a - k_b)x) + cos ((k_a + k_b)x)]
\end{aligned}
$$

Notice that the non-linear function has a low-frequency part ($k_a - k_b$) and a high-frequency part ($k_a + k_b$), where the higher frequency part has a frequency greater than either of the input signals.
> For input signals of equal frequencies, the output signal can have variability at twice the input frequency (**for this particular non-linearity**).

*Imagine if you had a grid that could just barely resolve $k$, then how would you ever resolve $2k$ on it? If anything, this high-frequency signal would get aliased to lower frequencies on a discrete grid.*

### A5: Impact of scales on flux
**Future thoughts.**
One way to assess the role of different scales on a non-linear flux is by considering the impact in wave-number space. For example for AB type non-linearity we can use the cross-spectra,

$$
\overline{AB} = \sum_{physical space}{AB} = \sum_{Fourier space} \hat{A} \hat{B},
$$
e.g. as done in Abernathey and Wortham 2013 or Kong and Jansen 2017 (MDPI-Fluids). Often the idea being that the dominant contributions peak at dominant eddy scales.

Can we do some similar decomposition for air-sea heat fluxes as well? Does the contribution similarly peak at the largest scales?