feat: Implement SAL when activating TIDE GENERATING FORCE

Summary

SAL or Self Attraction and Loading is a really important feature for Global Models. This physical process accounts for the retro-active effects of the Earth Tides on the Oceans.

I would need info on where we should implement SAL-induced modified gravity field.

This issue follows a post on the forum where I show the comparison between options 0 and 2 to understand the effects of reduced gravity on global tide model.

Why is this feature useful?

1. We need to be able to run surge (so meteo forcing) contemporaneously with tides.
2. Reducing gravity in the CAS file is not a good practice and most probably influences other calculations
3. We would also get closer to the real physics happening at planetary scale, providing better estimations of water levels worldwide.
4. Finally, applying SAL would be a major step-up and would potentially make TELEMAC compete with the most renowned global models (like SCHISM, ADCIRC or even FES..). Leading potentially to important publications, benchmarks and further improved physics of the global model (tidal+storm surge 3D baroclinic being the end game).

How to implement it?

This is where I would need some help.

CAS file

There could be - like done in SCHISM - 4 different options for the activation of SAL:

• 0 no SAL
• 1 interpolated from a tidal database (like FES)
• 2 constant coefficient attenuation for the gravity
• 3 parametric law based on Stepanov & Hughes - 2004

I posted on the forum the comparison between options 0 and 2 to understand the effects of reduced gravity on global tide model.

Hence why I'd be tempted to give it only 2 options:

• 0 no SAL
• 1 parametric law based on Stepanov & Hughes - 2004

SAL

The implementation of SAL is straightforward: the gravity force is modified by a factor (1 - beta) with beta being a polynomial regression. Here the illustration of the law:

dp1 = np.linspace(0, 6500,100)
tmp1 = np.sqrt(dp1)
beta_bar = -9.8169e-3 + 1.8289e-3*tmp1 + 4.3787e-4*dp1 \
           - 2.9042e-5*dp1*tmp1 + 6.6038e-7*dp1**2 \
           - 4.7393e-9*dp1**2*tmp1 - 1.9354e-11*dp1**3 \
           + 2.6969e-13*dp1**3*tmp1

beta_bar=np.maximum(0,beta_bar)

df = pd.DataFrame({"depth": dp1,"beta_bar": beta_bar})
df.hvplot(x="depth",y="beta_bar", label="beta_bar") * hv.Curve(((0,6500),(0.12,0.12)),label="loadtide_coef")

Its implementation could even be in MARAST

GRAV

Here's where things get more complex, the modified gravity field (let's call it GRAV2), would need to become an output of MARAST function (additionally to the already existing GRAV input).

GRAV2 would then be used in other functions. In SCHISM, GRAV2 is implemented in the bottom drag coefficient (l. 2061), in the density at boundaries (l. 5449) and the wave equation (l. 5802)

I need someone to point me towards where these changes should happen, and maybe physical processes that I - most probably - have forgotten.

Inviting @alexander.breugem, @christophe.coulet, @sebastien.bourban @noemie.durand and @michael.turnbull to share thoughts on this. Feel free to ping other people interested in large/global models as well.

Additional resources

• Stepanov & Hughes - 2004
• Forum post
• I think it's time to also have a proper test case for global models, to evaluate the implementation of new features. I will try to do this beforehand. More info on FESOM's pimesh, considered for the test cases, here.