Assignment 4: Basin-scale waves and overturning circulation

For this assignment you will do quantiative analysis of basin-scale waves (see Assignment 3: Basin-scale flows and waves), and a qualitative analysis of overturning circulations.

Q1.1 Super-inertial Poincarre waves.

Q1.1.a Determine the wavelength and frequency of waves

From the poincarre12.py simulation, estimate the wavelength and frequency of the waves produced.

Estimate the phase speed and compare to theory. Estimate the group speed and compare to theory. Recall that \(f=10^{-4} \ \mathrm{rad\,s^{-1}}\), and that the high-frequency wave speed is given by \(\sqrt{\frac{g'h_1 h_2}{h_1 + h_2}}\), where \(h_1 = 60 \ \mathrm{m}\), \(h_2 = 140\ \mathrm{m}\), and \(g' = 0.2 g / 1000 = 1.96\times10^{-3}\ \mathrm{m\,s^{-2}}\)

import xarray as xr
import numpy as np
import matplotlib.pyplot as plt

def subplots_shared(n, m, sharex=True, sharey=True, layout='constrained', figsize=None):
    return plt.subplots(n, m, sharex=sharex, sharey=sharey, layout=layout, figsize=figsize)

def pcolormesh_sym(x, y, z, vbound, ax=None, cmap='RdBu_r'):
    if ax is None:
        ax = plt.gca()
    return ax.pcolormesh(x, y, z, vmin=-vbound, vmax=vbound, cmap=cmap, rasterized=True)

def tohours(Time):
    return Time.values.astype('float') / 1e9 /3600

# use the matplotlib widget
%matplotlib widget

fig, ax = plt.subplots()
with xr.open_dataset('Poincarre12.snapshot.nc') as ds:
    ds = ds.sel(yt=100, yu=100, method='nearest')
    time = tohours(ds.Time)
    pc = ax.pcolormesh(ds.xt, time, ds.v.isel(zt=-2) * (ds.xt-320) / 25,
                      vmin=-0.01, vmax=0.01, cmap='RdBu_r')
    ax.contour(ds.xt, time, ds.v.isel(zt=-2) * (ds.xt-320) / 25, levels=[-1e-3, 1e-3])
    fig.colorbar(pc, ax = ax, shrink=0.4, extend='both', label=r'$v \ \mathrm{[m\,s^{-1}]}$')
    ax.set_yticks(np.arange(0, 100, 12))
    ax.set_xlabel('x [km]')
    ax.set_ylabel('Time [h]')
    ax.grid('on')
    
    T = 12  # h
    lam = 242-165  # km
    #ax.plot(ds.xt, (300 - ds.xt) * 12 / lam + 32, ls='--', color='C0')
    #ax.plot(ds.xt, (300 - ds.xt) * 12 / lam + 32 + 12*3, ls='--', color='C0')
    ax.set_ylim([0, 96])
    
    # group speed:
    dx = 221 - 141
    dt = 28 - 54
    
    cg = dx/dt
    ax.plot(ds.xt , -(300-ds.xt) / cg , ls='--', color='C1')
    
    gp = 20 * 0.167 * 9.81 / 1024
    h1 = 60
    h2 = 140
    ch = np.sqrt(gp * h1 * h2 / (h1 + h2))
    print('Chigh', ch)
    om = 2 * np.pi / 12 / 3600
    f = 1e-4
    k = np.sqrt((om**2 - f**2) / ch**2)
    cp = om / k  # m/s
    print('phase speed', cp, 'm/s')
    cp = cp / 1000 * 3600  # km / h
    print('phase speed', cp, 'km/h')
    
    ax.plot(ds.xt, (300 - ds.xt) / cp + 32, ls=':', color='C0')
    ax.plot(ds.xt, (300 - ds.xt) / cp + 32 + 12, ls=':', color='C0')
    
    cg = ch**2 / om * k
    cg = cg * 3.6   # km/h
    ax.plot(ds.xt +30 , (300-ds.xt) / cg, ls=':', color='C1')

    
    

/var/folders/p1/6grcm4fx3tx_dyzph3cbwk4c0000gn/T/ipykernel_10318/3351537445.py:2: FutureWarning: In a future version of xarray decode_timedelta will default to False rather than None. To silence this warning, set decode_timedelta to True, False, or a 'CFTimedeltaCoder' instance.
  with xr.open_dataset('Poincarre12.snapshot.nc') as ds:

Chigh 1.1592641456436925
phase speed 1.596479878032792 m/s
phase speed 5.747327560918051 km/h

The phase speed is quite easy to pick off as the slope of the phase lines. I found it easiest to contour a very small amplitude to pick off the propagation of energy across the basin, or the group speed. Note that the group speed is substantially less than the phase speed.

Theoretically,

\[ \omega = \sqrt{f^2 + c_H^2 k^2} \]

\(c_H = 1.15 \ \mathrm{m\,s^{-1}}\). We know that \(\omega = 2 \pi / 12 h = 1.45\times10^{-4}\ \mathrm{rad\,s^{-1}}\) so we can invert for \(k\) and then get the phase speed as \(\omega / k = \lambda / T = 1.74 \ \mathrm{m\,s^{-1}}\) (note that plot above is km/h so convert). The eyball estimate and the dashed line are in very good agreement.

The group speed is given by \(\frac{\partial \omega}{\partial k} = \frac{c_H^2}{c_p} = 0.921 \ \mathrm{m\,s^{-1}}\). Note that this groups speed is somewhat faster than what I picked off by eye. The reason is that the amplitude of the spreading waves drops off with distance from the paddle because of cylindrical spreading. In the above plots I tried to correct for this by multiplying the velocity signal by \(r/r_0\).

We can also see by looking at the plots we made last time that the wave signal interferes with itself as Kelvin waves are created when the Poincarre waves hit the southern boundary, and carry energy back to the paddle region.

Q1.2a Kelvin wave

Do as above, but for the Kelvin wave: Kelvin40.py.

fig, ax = plt.subplots()
with xr.open_dataset('Kelvin40.snapshot.nc') as ds:
    ds = ds.isel(xt=-3, xu=-3)
    time = tohours(ds.Time)
    pc = ax.pcolormesh(ds.yt, time, ds.v.isel(zt=-2),
                      vmin=-0.01*2, vmax=0.01*2, cmap='RdBu_r')
    # ax.contour(ds.xt, time, ds.v.isel(zt=-2), levels=[-1e-4, 1e-4])
    fig.colorbar(pc, ax = ax, shrink=0.4, extend='both', label=r'$v \ \mathrm{[m\,s^{-1}]}$')
    ax.set_yticks(np.arange(0, 100, 12))
    ax.set_xlabel('y [km]')
    ax.set_ylabel('Time [h]')
    ax.grid('on')
    
    T = 12  # h
    ax.set_ylim([0, 96])
    
    ax.plot(ds.yt, (ds.yt) / ch / 3.6, ls=':', color='C0')
    ax.plot(ds.yt, (ds.yt) / ch / 3.6 + 40, ls=':', color='C0')
    

/var/folders/p1/6grcm4fx3tx_dyzph3cbwk4c0000gn/T/ipykernel_10318/3153273277.py:2: FutureWarning: In a future version of xarray decode_timedelta will default to False rather than None. To silence this warning, set decode_timedelta to True, False, or a 'CFTimedeltaCoder' instance.
  with xr.open_dataset('Kelvin40.snapshot.nc') as ds:

The Kelvin waves very consistenly move at the high-frequency phase speed, and that phase speed is also the group speed.

Q1.2b Kelvin wave assumptions:

Check the assumptions of Kelvin wave theory. How large is the along-sidewall velocity versus cross-sidewall?
What happens when the wave reaches a corner? What is the Rossby radius and how does it compare to the observed signals?

fig, ax = subplots_shared(1, 2)
with xr.open_dataset('Kelvin40.snapshot.nc') as ds:
    ds = ds.isel(Time=110, zt=-2, zw=-2)
    pcolormesh_sym(ds.xu, ds.yt, ds.u, vbound=0.01, ax=ax[0])
    pc = pcolormesh_sym(ds.xt, ds.yu, ds.v, vbound=0.01, ax=ax[1])
    for a in ax:
        a.contour(ds.xt, ds.yt, ds.psi, levels=np.arange(-0.103, 0.1, 0.005), vmin=-0.01, vmax=0.01, cmap='RdBu_r')
    
        a.set_aspect(1)
    fig.colorbar(pc, ax=ax, shrink=0.4, extend='both')
    ax[0].set_ylabel('Y [km]')
    ax[1].set_xlabel('X [km]')
    ax[0].set_title('u [m/s]')
    ax[1].set_title('v [m/s]')
    
    
    Ro = ch / f / 1000
    print(Ro)
    ax[1].axvline(320-Ro)
    ax[0].axhline(320-Ro)
    
    

11.592641456436924

/var/folders/p1/6grcm4fx3tx_dyzph3cbwk4c0000gn/T/ipykernel_10318/107751205.py:2: FutureWarning: In a future version of xarray decode_timedelta will default to False rather than None. To silence this warning, set decode_timedelta to True, False, or a 'CFTimedeltaCoder' instance.
  with xr.open_dataset('Kelvin40.snapshot.nc') as ds:

So the Kelvin wave signal is as expected - most of the velocity is parallel to the coasts, and there is almost no velocity perpendicular. The velcoity signal falls off in approximately one Rossby radius (\(c_H / f = 11.6 \ \mathrm{km}\)).

Q1.2c geostrophic balance?

Calculate the north-south geostrophic flow and compare to the flow in the Kelvin waves in the upper layer. Comment on any difference.

fig, ax = subplots_shared(1, 2, figsize=(3, 6))
with xr.open_dataset('Kelvin40.snapshot.nc') as ds:
    ds = ds.isel(Time=110, zt=-2, zw=-2)
    vgeo = ds.psi.differentiate(coord='xt') / f / 1000
    pc = pcolormesh_sym(ds.xt, ds.yu, vgeo, vbound=0.01, ax=ax[0])
    
    pc = pcolormesh_sym(ds.xt, ds.yu, ds.v, vbound=0.01, ax=ax[1])
    for a in ax:
    
        a.set_aspect(1)
    fig.colorbar(pc, ax=ax, shrink=0.4, extend='both')
    ax[0].set_ylabel('Y [km]')
    ax[1].set_xlabel('X [km]')
    ax[0].set_title('$v_{geo}$ [m/s]')
    ax[1].set_title('v [m/s]')
    ax[0].set_xlim([280, 320])
    
    ax[1].axvline(320-Ro)
    ax[0].axhline(320-Ro)
    
    

/var/folders/p1/6grcm4fx3tx_dyzph3cbwk4c0000gn/T/ipykernel_10318/939435257.py:2: FutureWarning: In a future version of xarray decode_timedelta will default to False rather than None. To silence this warning, set decode_timedelta to True, False, or a 'CFTimedeltaCoder' instance.
  with xr.open_dataset('Kelvin40.snapshot.nc') as ds:

So, the velocities are not quite in geostophic balance - with the simulated velocities somewhat weaker. This makes sense, as the flow has not reached equilibrium, and indeed it cannot because of the wall. The velocities cause the sea level to start rising between the highs and lows creating a new high, and this slows the flow down in the existing high.

Q1.2 Rossby waves

Q1.2.1 expected wavelengths

We ran three simulations for Rossby waves T=200, 300 and 400 h. Considering just phase propagation in x, along the center of the model domain, what are the expected x-direction wavelengths? Recall that the water depth is 2000 m, and that \(\beta = 10^{-11}\ \mathrm{m\,s^{-1}}\).

The dispersion relation witn \(m=0\) is

\[ \begin{align}\begin{aligned} \omega = -\beta_0R^2 \frac{l}{1 + R^2l^2} $$ which is a quadratic in $l$ with solution:\\$$l = -\frac{\beta_0}{2\omega} \pm \frac{1}{2}\sqrt{\frac{\beta_0^2}{\omega^2} - \frac{4}{R^2}} \end{aligned}\end{align} \]

beta = 1e-11
R = np.sqrt(9.81 * 2000) / f
for T in [200, 300, 400]:
    omega = np.pi * 2 / T / 3600

    lp = -(beta/2/omega) + 0.5*np.sqrt(beta**2 / omega**2 - 4 / R**2)
    lm = -(beta/2/omega) - 0.5*np.sqrt(beta**2 / omega**2 - 4 / R**2)
    print(f'T {T}: {2 * np.pi / lp / 1e3} m, {2 * np.pi / lm / 1e3} km')
    
    # group speeds for later:
    cgp = beta * (-R**2 + R**4*lp**2)/(1+R**2*lp**2)**2
    cgm = beta * (-R**2 + R**4*lm**2)/(1+R**2*lm**2)**2

T 200: nan m, nan km
T 300: -16493.370309850576 m, -4696.229690149424 km
T 400: -25176.219981350416 m, -3076.5800186495926 km

/var/folders/p1/6grcm4fx3tx_dyzph3cbwk4c0000gn/T/ipykernel_10318/3227838967.py:6: RuntimeWarning: invalid value encountered in sqrt
  lp = -(beta/2/omega) + 0.5*np.sqrt(beta**2 / omega**2 - 4 / R**2)
/var/folders/p1/6grcm4fx3tx_dyzph3cbwk4c0000gn/T/ipykernel_10318/3227838967.py:7: RuntimeWarning: invalid value encountered in sqrt
  lm = -(beta/2/omega) - 0.5*np.sqrt(beta**2 / omega**2 - 4 / R**2)

So T=200 h has no Rossby wave solutions. Note that T=300,400 h, the longer wavelengths are very long indeed. The shorter wavelengths are shortest for T=400 h.

Q1.2.2 Hovmoller diagram

Show using a Hovmoller diagram for all three Rossby wave simulations that the shorter wavelengths are close to expected. Comment on the energy propagating to the left.

fig, axs = subplots_shared(1, 3)
for nn, T, lambm, lambp in zip([0, 1, 2], 
                               [200, 300, 400], 
                               [np.nan, -4696.2, -3076.6], 
                               [np.nan, -16493.37, -25176.22]):
    ax = axs[nn]
    with xr.open_dataset(f'Rossby{T}.snapshot.nc') as ds:
        ds = ds.isel(yt=65, yu=65, zt=-2, zw=-2)
        time = tohours(ds.Time)
        pcolormesh_sym(ds.xt, time, ds.v, vbound=0.002, ax=ax)
        ax.plot(ds.xt, ds.xt / lambm * T + 5*T, color='C0')
        ax.plot(ds.xt, ds.xt / lambm * T + 6*T, color='C0')
        if nn > 0:
            ax.text(ds.xt.mean(), 4*T, 'Phase speed', color='C0')
        ax.plot(ds.xt, ds.xt / lambp * T + 5*T, color='C1')
        ax.plot(ds.xt, ds.xt / lambp * T + 6*T, color='C1')
        if nn > 0:
            ax.text(ds.xt.mean(), 5.2*T, 'Phase speed', color='C1')
        if nn == 2:
            ax.plot(ds.xt, (ds.xt - ds.xt.mean())/ cgp / 3.6 + T / 2, color='C2', ls='--')
            ax.plot(ds.xt, (ds.xt - ds.xt.mean())/ cgm / 3.6 + T / 2, color='C2', label='Group speed')
            ax.text(ds.xt.mean(), T, 'Group speed', color='C2')
        ax.set_ylim(bottom=0)
        ax.set_title(f'{T} h')

/var/folders/p1/6grcm4fx3tx_dyzph3cbwk4c0000gn/T/ipykernel_10318/2690420010.py:7: FutureWarning: In a future version of xarray decode_timedelta will default to False rather than None. To silence this warning, set decode_timedelta to True, False, or a 'CFTimedeltaCoder' instance.
  with xr.open_dataset(f'Rossby{T}.snapshot.nc') as ds:

/var/folders/p1/6grcm4fx3tx_dyzph3cbwk4c0000gn/T/ipykernel_10318/2690420010.py:7: FutureWarning: In a future version of xarray decode_timedelta will default to False rather than None. To silence this warning, set decode_timedelta to True, False, or a 'CFTimedeltaCoder' instance.
  with xr.open_dataset(f'Rossby{T}.snapshot.nc') as ds:

/var/folders/p1/6grcm4fx3tx_dyzph3cbwk4c0000gn/T/ipykernel_10318/2690420010.py:7: FutureWarning: In a future version of xarray decode_timedelta will default to False rather than None. To silence this warning, set decode_timedelta to True, False, or a 'CFTimedeltaCoder' instance.
  with xr.open_dataset(f'Rossby{T}.snapshot.nc') as ds:

The 400-h simulation is the most clear, with the best agreement between the forcing and the response. Note that in these simulations, the forcing region is relatively small compared to the wavelengths, so the response is not always very good.

To the east of the forcing, there is a strong response. Note that the phase velocity is still to the left, but the response grows to the right. There is a slower phase propagation to the left, but Rossby waves refelct off the wester boundary and interfere with the westward radiating signal.

Note that there is a response in the 200-h forcing, and that it looks a lot like the 400-h forcing. This is an example of frequency leakage, where the 200-h forcing is not a perfect sine wave, and hence leads to some lowere frequency forcing and generates weak Rossby waves at the 400-h period.

Q1.2.3 Group speed

Calculate the x-direction group speed for both sets of waves in the 400-h simulation and compare to the simulation’s Hovmoller diagram (I did this by just adding to the Hovmoller diagram above and replotting it.

From

\[ \omega = -\beta_0R^2 \frac{l}{1 + R^2l^2} $$ we can calculate $$c_{gx} = \frac{d\omega}{dl} = \frac{-\beta_0 R^2 + \beta_0 R^4 l^2}{(1 + R^2 l^2)^2}\]

cgp = beta * (-R**2 + R**4*lp**2)/(1+R**2*lp**2)**2
cgm = beta * (-R**2 + R**4*lm**2)/(1+R**2*lm**2)**2

print(f'Positive branch: wavelength = {2*np.pi/lp / 1000:1.2f} km, group speed = {cgp *3.6:1.2f} km/h ')
print(f'Negative branch: wavelength = {2*np.pi/lm / 1000:1.2f} km, group speed  = {cgm * 3.6:1.2f} km/h')
print(cgm)

Positive branch: wavelength = -25176.22 km, group speed = -49.23 km/h 
Negative branch: wavelength = -3076.58 km, group speed  = 6.02 km/h
1.6712038454938625

So we see that the shorter wavelenth has a group speed to the right at 6 km/h, whereas the longer wavelength has a group speed to the left at 49 km/h. If we add these speeds to the plots above, then we get reasonable estimates for the speed at which energy propagates across each basin..

Q1.2.4 Compare dq/dt to q

For Rossby400.py, calculate the time derivative of the relative vorticity, and compare to the relative vorticity and argue if this make sense for the direction of propagation of the Rossby wave.

with xr.open_dataset('Rossby400.snapshot.nc') as ds:
    ds = ds.isel(Time=slice(30,32)).isel(zt=-1)

    # calculate the vorticity: dv/dx-du/dy:
    vort = ds.v.diff(dim='xu').values + ds.u.diff(dim='yu').values
    vort = vort / 40000
    dvort_dt = np.diff(vort, axis=0)[0] / 24 /3600

    fig, ax = plt.subplots(2, 1, sharex=True, sharey=True, layout='constrained')
        
    pc = ax[0].pcolormesh(ds.xt, ds.yt, vort[0], vmin=-1e-7, vmax=1e-7, cmap='RdBu_r')
    ax[0].grid('on')
    ax[0].set(title='$\zeta\ \mathrm{[s^{-1}]}$', ylabel='y [km]')
    fig.colorbar(pc, ax=ax[0], shrink=0.7)
    
    pc=ax[1].pcolormesh(ds.xt, ds.yt, dvort_dt, vmin=-4e-13, vmax=4e-13, cmap='RdBu_r')
    ax[1].grid('on')
    ax[1].set(title='$d\zeta/dt\ \mathrm{[s^{-2}]}$', xlabel='x [km]')
    fig.colorbar(pc, ax=ax[1], shrink=0.7)
    

<>:13: SyntaxWarning: invalid escape sequence '\z'
<>:18: SyntaxWarning: invalid escape sequence '\z'
<>:13: SyntaxWarning: invalid escape sequence '\z'
<>:18: SyntaxWarning: invalid escape sequence '\z'
/var/folders/p1/6grcm4fx3tx_dyzph3cbwk4c0000gn/T/ipykernel_10318/1483438720.py:13: SyntaxWarning: invalid escape sequence '\z'
  ax[0].set(title='$\zeta\ \mathrm{[s^{-1}]}$', ylabel='y [km]')
/var/folders/p1/6grcm4fx3tx_dyzph3cbwk4c0000gn/T/ipykernel_10318/1483438720.py:18: SyntaxWarning: invalid escape sequence '\z'
  ax[1].set(title='$d\zeta/dt\ \mathrm{[s^{-2}]}$', xlabel='x [km]')
/var/folders/p1/6grcm4fx3tx_dyzph3cbwk4c0000gn/T/ipykernel_10318/1483438720.py:1: FutureWarning: In a future version of xarray decode_timedelta will default to False rather than None. To silence this warning, set decode_timedelta to True, False, or a 'CFTimedeltaCoder' instance.
  with xr.open_dataset('Rossby400.snapshot.nc') as ds:

We can see from this that the null in vorticity at x=7000 km has decreasing vorticity maximum. So we expect the negative vorticity anomaly to move to the left.

Q2 Overturning Circulation: Qualitative

These simulations take a few hours to run on my machine, and I want you to look at a few, so I’ll provide the netcdf files. The simulations are in Overturning.py, and then Hosing.py and StopHosing.py. The files are provided in OverturningRuns.zip (Note 776 Mb!).

Physically, the models are solved on a spherical grid so the coordinates are in degrees longitude (xt, xu) and latitude (yt,yu). The “wedge” is 48 degrees in longitude, and from the equator to 72 degrees N in latitude. There are 18 vertical levels used for these simulations, with the vertical grid spacing between the levels exponentially growing from 113 m to 318 m at the seafloor.

All the models use the same external forcing of a sin-wave temperature at the surface, going from 28 degrees at the equator to 4 degress at 72 degrees N.

The models are run with isopycnal lateral diffusion in an attempt to reduce the “Veronis effect”.

Q2.1 Plot the forcing

Plot the surface temperature for OverturningBeta4e6Coarse.snapshot.nc. How close is it to the forcing described above?

td = 'OverturningBeta10e5Coarse'
with xr.open_dataset(td+'.snapshot.nc') as ds:
    display(ds)
    fig, ax = plt.subplots()
    pc = ax.pcolormesh(ds.xt, ds.yt, ds.temp.isel(zt=-1, Time=-1))
    fig.colorbar(pc, label='T $[^oC]$')
    ax.set_xlabel('Lon [deg]')
    ax.set_ylabel('Lon [deg]')
    ax.set_title('Surface temp')

/var/folders/p1/6grcm4fx3tx_dyzph3cbwk4c0000gn/T/ipykernel_10318/4043134810.py:2: FutureWarning: In a future version of xarray decode_timedelta will default to False rather than None. To silence this warning, set decode_timedelta to True, False, or a 'CFTimedeltaCoder' instance.
  with xr.open_dataset(td+'.snapshot.nc') as ds:

<xarray.Dataset> Size: 178MB
Dimensions:            (xt: 12, xu: 12, yt: 18, yu: 18, zt: 18, zw: 18,
                        tensor1: 2, tensor2: 2, Time: 400)
Coordinates:
  * xt                 (xt) float64 96B -2.0 2.0 6.0 10.0 ... 34.0 38.0 42.0
  * xu                 (xu) float64 96B 0.0 4.0 8.0 12.0 ... 32.0 36.0 40.0 44.0
  * yt                 (yt) float64 144B -2.0 2.0 6.0 10.0 ... 58.0 62.0 66.0
  * yu                 (yu) float64 144B 0.0 4.0 8.0 12.0 ... 60.0 64.0 68.0
  * zt                 (zt) float64 144B -3.809e+03 -3.428e+03 ... -178.6 -47.76
  * zw                 (zw) float64 144B -3.619e+03 -3.264e+03 ... -113.2 0.0
  * tensor1            (tensor1) float64 16B 0.0 1.0
  * tensor2            (tensor2) float64 16B 0.0 1.0
  * Time               (Time) timedelta64[ns] 3kB 360 days ... -57768 days +0...
Data variables: (12/33)
    dxt                (xt) float64 96B ...
    dxu                (xu) float64 96B ...
    dyt                (yt) float64 144B ...
    dyu                (yu) float64 144B ...
    dzt                (zt) float64 144B ...
    dzw                (zw) float64 144B ...
    ...                 ...
    kappaM             (Time, zt, yt, xt) float64 12MB ...
    kappaH             (Time, zw, yt, xt) float64 12MB ...
    surface_taux       (Time, yt, xu) float64 691kB ...
    surface_tauy       (Time, yu, xt) float64 691kB ...
    forc_rho_surface   (Time, yt, xt) float64 691kB ...
    psi                (Time, yt, xt) float64 691kB ...
Attributes:
    date_created:       2023-09-07T15:42:47.965547
    veros_version:      1.5.0
    setup_identifier:   OverturningBeta10e5Coarse
    setup_description:  
    setup_settings:     {"identifier": "OverturningBeta10e5Coarse", "descript...
    setup_file:         /Users/jklymak/Dropbox/Teaching/Eos431Phy441/source/A...
    setup_code:         #!/usr/bin/env python\n\n"""\n"""\n\n__VEROS_VERSION_...

xarray.Dataset

• Dimensions:
  • xt: 12
  • xu: 12
  • yt: 18
  • yu: 18
  • zt: 18
  • zw: 18
  • tensor1: 2
  • tensor2: 2
  • Time: 400
• Coordinates: (9)
  • xt
    (xt)
    float64
    -2.0 2.0 6.0 ... 34.0 38.0 42.0

    long_name :

    Zonal coordinate (T)

    units :

    degrees_east

    array([-2.,  2.,  6., 10., 14., 18., 22., 26., 30., 34., 38., 42.])

  • xu
    (xu)
    float64
    0.0 4.0 8.0 12.0 ... 36.0 40.0 44.0

    long_name :

    Zonal coordinate (U)

    units :

    degrees_east

    array([ 0.,  4.,  8., 12., 16., 20., 24., 28., 32., 36., 40., 44.])

  • yt
    (yt)
    float64
    -2.0 2.0 6.0 ... 58.0 62.0 66.0

    long_name :

    Meridional coordinate (T)

    units :

    degrees_north

    array([-2.,  2.,  6., 10., 14., 18., 22., 26., 30., 34., 38., 42., 46., 50.,
           54., 58., 62., 66.])

  • yu
    (yu)
    float64
    0.0 4.0 8.0 12.0 ... 60.0 64.0 68.0

    long_name :

    Meridional coordinate (U)

    units :

    degrees_north

    array([ 0.,  4.,  8., 12., 16., 20., 24., 28., 32., 36., 40., 44., 48., 52.,
           56., 60., 64., 68.])

  • zt
    (zt)
    float64
    -3.809e+03 -3.428e+03 ... -47.76

    long_name :

    Vertical coordinate (T)

    units :

    m

    positive :

    up

    array([-3809.44645 , -3428.33935 , -3099.777178, -2767.592705, -2484.58057 ,
           -2194.806041, -1951.280227, -1698.269941, -1488.973945, -1267.833769,
           -1088.210846,  -894.698165,  -740.798173,  -571.235155,  -439.633805,
            -290.832221,  -178.56107 ,   -47.757115])

  • zw
    (zw)
    float64
    -3.619e+03 -3.264e+03 ... 0.0

    long_name :

    Vertical coordinate (W)

    units :

    m

    positive :

    up

    array([-3618.8929  , -3264.058264, -2933.684941, -2626.086637, -2339.693305,
           -2073.043134, -1824.775084, -1593.621943, -1378.403857, -1178.022308,
            -991.454505,  -817.748169,  -656.016664,  -505.43448 ,  -365.233013,
            -234.696645,  -113.159092,     0.      ])

  • tensor1
    (tensor1)
    float64
    0.0 1.0

    long_name :

    tensor1

    units :

    array([0., 1.])

  • tensor2
    (tensor2)
    float64
    0.0 1.0

    long_name :

    tensor2

    units :

    array([0., 1.])

  • Time
    (Time)
    timedelta64[ns]
    360 days ... -57768 days +02:58:...

    long_name :

    Time

    time_origin :

    01-JAN-1900 00:00:00

    array([   31104000000000000,   342144000000000000,   653184000000000000, ...,
           -5613224515966860288, -5302184515966860288, -4991144515966860288],
          shape=(400,), dtype='timedelta64[ns]')

• Data variables: (33)
  • dxt
    (xt)
    float64
    ...

    long_name :

    Zonal T-grid spacing

    units :

    m

    [12 values with dtype=float64]

  • dxu
    (xu)
    float64
    ...

    long_name :

    Zonal U-grid spacing

    units :

    m

    [12 values with dtype=float64]

  • dyt
    (yt)
    float64
    ...

    long_name :

    Meridional T-grid spacing

    units :

    m

    [18 values with dtype=float64]

  • dyu
    (yu)
    float64
    ...

    long_name :

    Meridional U-grid spacing

    units :

    m

    [18 values with dtype=float64]

  • dzt
    (zt)
    float64
    ...

    long_name :

    Vertical spacing (T)

    units :

    m

    [18 values with dtype=float64]

  • dzw
    (zw)
    float64
    ...

    long_name :

    Vertical spacing (W)

    units :

    m

    [18 values with dtype=float64]

  • ht
    (yt, xt)
    float64
    ...

    long_name :

    Total depth (T)

    units :

    m

    [216 values with dtype=float64]

  • hu
    (yt, xu)
    float64
    ...

    long_name :

    Total depth (U)

    units :

    m

    [216 values with dtype=float64]

  • hv
    (yu, xt)
    float64
    ...

    long_name :

    Total depth (V)

    units :

    m

    [216 values with dtype=float64]

  • beta
    (yt, xt)
    float64
    ...

    long_name :

    Change of Coriolis freq.

    units :

    1/(ms)

    [216 values with dtype=float64]

  • area_t
    (yt, xt)
    float64
    ...

    long_name :

    Area of T-box

    units :

    m^2

    [216 values with dtype=float64]

  • area_u
    (yt, xu)
    float64
    ...

    long_name :

    Area of U-box

    units :

    m^2

    [216 values with dtype=float64]

  • area_v
    (yu, xt)
    float64
    ...

    long_name :

    Area of V-box

    units :

    m^2

    [216 values with dtype=float64]

  • rho
    (Time, zt, yt, xt)
    float64
    ...

    long_name :

    Density

    units :

    kg/m^3

    [1555200 values with dtype=float64]

  • prho
    (Time, zt, yt, xt)
    float64
    ...

    long_name :

    Potential density

    units :

    kg/m^3

    [1555200 values with dtype=float64]

  • int_drhodT
    (Time, zt, yt, xt)
    float64
    ...

    long_name :

    Der. of dyn. enthalpy by temperature

    units :

    kg / (m^2 deg C)

    [1555200 values with dtype=float64]

  • int_drhodS
    (Time, zt, yt, xt)
    float64
    ...

    long_name :

    Der. of dyn. enthalpy by salinity

    units :

    kg / (m^2 g / kg)

    [1555200 values with dtype=float64]

  • Nsqr
    (Time, zw, yt, xt)
    float64
    ...

    long_name :

    Square of stability frequency

    units :

    1/s^2

    [1555200 values with dtype=float64]

  • Hd
    (Time, zt, yt, xt)
    float64
    ...

    long_name :

    Dynamic enthalpy

    units :

    m^2/s^2

    [1555200 values with dtype=float64]

  • temp
    (Time, zt, yt, xt)
    float64
    ...

    long_name :

    Temperature

    units :

    deg C

    [1555200 values with dtype=float64]

  • salt
    (Time, zt, yt, xt)
    float64
    ...

    long_name :

    Salinity

    units :

    g/kg

    [1555200 values with dtype=float64]

  • forc_temp_surface
    (Time, yt, xt)
    float64
    ...

    long_name :

    Surface temperature flux

    units :

    m deg C/s

    [86400 values with dtype=float64]

  • forc_salt_surface
    (Time, yt, xt)
    float64
    ...

    long_name :

    Surface salinity flux

    units :

    m g/s kg

    [86400 values with dtype=float64]

  • u
    (Time, zt, yt, xu)
    float64
    ...

    long_name :

    Zonal velocity

    units :

    m/s

    [1555200 values with dtype=float64]

  • v
    (Time, zt, yu, xt)
    float64
    ...

    long_name :

    Meridional velocity

    units :

    m/s

    [1555200 values with dtype=float64]

  • w
    (Time, zw, yt, xt)
    float64
    ...

    long_name :

    Vertical velocity

    units :

    m/s

    [1555200 values with dtype=float64]

  • p_hydro
    (Time, zt, yt, xt)
    float64
    ...

    long_name :

    Hydrostatic pressure

    units :

    m^2/s^2

    [1555200 values with dtype=float64]

  • kappaM
    (Time, zt, yt, xt)
    float64
    ...

    long_name :

    Vertical viscosity

    units :

    m^2/s

    [1555200 values with dtype=float64]

  • kappaH
    (Time, zw, yt, xt)
    float64
    ...

    long_name :

    Vertical diffusivity

    units :

    m^2/s

    [1555200 values with dtype=float64]

  • surface_taux
    (Time, yt, xu)
    float64
    ...

    long_name :

    Surface wind stress

    units :

    N/m^2

    [86400 values with dtype=float64]

  • surface_tauy
    (Time, yu, xt)
    float64
    ...

    long_name :

    Surface wind stress

    units :

    N/m^2

    [86400 values with dtype=float64]

  • forc_rho_surface
    (Time, yt, xt)
    float64
    ...

    long_name :

    Surface density flux

    units :

    kg / (m^2 s)

    [86400 values with dtype=float64]

  • psi
    (Time, yt, xt)
    float64
    ...

    long_name :

    Surface pressure

    units :

    m^2/s^2

    [86400 values with dtype=float64]

• Indexes: (9)
  • xt
    PandasIndex

    PandasIndex(Index([-2.0, 2.0, 6.0, 10.0, 14.0, 18.0, 22.0, 26.0, 30.0, 34.0, 38.0, 42.0], dtype='float64', name='xt'))

  • xu
    PandasIndex

    PandasIndex(Index([0.0, 4.0, 8.0, 12.0, 16.0, 20.0, 24.0, 28.0, 32.0, 36.0, 40.0, 44.0], dtype='float64', name='xu'))

  • yt
    PandasIndex

    PandasIndex(Index([-2.0,  2.0,  6.0, 10.0, 14.0, 18.0, 22.0, 26.0, 30.0, 34.0, 38.0, 42.0,
           46.0, 50.0, 54.0, 58.0, 62.0, 66.0],
          dtype='float64', name='yt'))

  • yu
    PandasIndex

    PandasIndex(Index([ 0.0,  4.0,  8.0, 12.0, 16.0, 20.0, 24.0, 28.0, 32.0, 36.0, 40.0, 44.0,
           48.0, 52.0, 56.0, 60.0, 64.0, 68.0],
          dtype='float64', name='yu'))

  • zt
    PandasIndex

    PandasIndex(Index([  -3809.44644995089,   -3428.33934985267, -3099.7771779861137,
            -2767.592704826615,  -2484.580569988048, -2194.8060409194495,
           -1951.2802265132486,  -1698.269940706693, -1488.9739448257642,
           -1267.8337690041126, -1088.2108463152017,  -894.6981645020714,
            -740.7981731809905,  -571.2351553493272, -439.63380474780934,
            -290.8322213281226, -178.56106956086296,  -47.75711480003247],
          dtype='float64', name='zt'))

  • zw
    PandasIndex

    PandasIndex(Index([  -3618.89289990178, -3264.0582639193917, -2933.6849414063645,
            -2626.086637407331, -2339.6933054537485,  -2073.043133716349,
           -1824.7750836099708, -1593.6219427662286, -1378.4038569149384,
           -1178.0223076596571,  -991.4545054086366,   -817.748168841531,
            -656.0166642651589,  -505.4344800485683,   -365.233013037966,
            -234.6966454444928, -113.15909218044771,                 0.0],
          dtype='float64', name='zw'))

  • tensor1
    PandasIndex

    PandasIndex(Index([0.0, 1.0], dtype='float64', name='tensor1'))

  • tensor2
    PandasIndex

    PandasIndex(Index([0.0, 1.0], dtype='float64', name='tensor2'))

  • Time
    PandasIndex

    PandasIndex(TimedeltaIndex([              '360 days 00:00:00',
                                 '3960 days 00:00:00',
                                 '7560 days 00:00:00',
                                '11160 days 00:00:00',
                                '14760 days 00:00:00',
                                '18360 days 00:00:00',
                                '21960 days 00:00:00',
                                '25560 days 00:00:00',
                                '29160 days 00:00:00',
                                '32760 days 00:00:00',
                    ...
                    '-90168 days +02:58:04.033139712',
                    '-86568 days +02:58:04.033139712',
                    '-82968 days +02:58:04.033139712',
                    '-79368 days +02:58:04.033139712',
                    '-75768 days +02:58:04.033139712',
                    '-72168 days +02:58:04.033139712',
                    '-68568 days +02:58:04.033139712',
                    '-64968 days +02:58:04.033139712',
                    '-61368 days +02:58:04.033139712',
                    '-57768 days +02:58:04.033139712'],
                   dtype='timedelta64[ns]', name='Time', length=400, freq=None))

• Attributes: (7)

  date_created :

  2023-09-07T15:42:47.965547

  veros_version :

  1.5.0

  setup_identifier :

  OverturningBeta10e5Coarse

  setup_description :

  setup_settings :

  {"identifier": "OverturningBeta10e5Coarse", "description": "", "nx": 12, "ny": 18, "nz": 18, "dt_mom": 3000.0, "dt_tracer": 86400.0, "runlen": 124416000000.0, "AB_eps": 0.1, "x_origin": 0.0, "y_origin": 0.0, "pi": 3.141592653589793, "radius": 6370000.0, "degtom": 111177.4733520388, "omega": 7.292123516990375e-05, "rho_0": 1024.0, "grav": 9.81, "coord_degree": true, "enable_cyclic_x": false, "eq_of_state_type": 1, "enable_implicit_vert_friction": true, "enable_explicit_vert_friction": false, "enable_hor_friction": true, "enable_hor_diffusion": false, "enable_biharmonic_friction": false, "enable_biharmonic_mixing": false, "enable_hor_friction_cos_scaling": true, "enable_ray_friction": false, "enable_bottom_friction": false, "enable_bottom_friction_var": false, "enable_quadratic_bottom_friction": false, "enable_tempsalt_sources": false, "enable_momentum_sources": false, "enable_superbee_advection": false, "enable_conserve_energy": true, "enable_store_bottom_friction_tke": false, "enable_store_cabbeling_heat": false, "enable_noslip_lateral": false, "enable_streamfunction": false, "A_h": 1758977.8451616622, "K_h": 1758977.8451616622, "r_ray": 0.0, "r_bot": 0.0, "r_quad_bot": 0.0, "hor_friction_cosPower": 1.0, "A_hbi": 1.313046875e+16, "K_hbi": 1.313046875e+16, "biharmonic_friction_cosPower": 0.0, "kappaH_0": 0.0001, "kappaM_0": 0.0001, "enable_neutral_diffusion": true, "enable_skew_diffusion": true, "enable_TEM_friction": false, "K_iso_0": 1000.0, "K_iso_steep": 500.0, "K_gm_0": 0.0, "iso_dslope": 0.005, "iso_slopec": 0.01, "enable_idemix": false, "tau_v": 172800.0, "tau_h": 1296000.0, "gamma": 1.57, "jstar": 5.0, "mu0": 0.3333333333333333, "enable_idemix_hor_diffusion": false, "enable_eke_diss_bottom": false, "enable_eke_diss_surfbot": false, "eke_diss_surfbot_frac": 1.0, "enable_idemix_superbee_advection": false, "enable_idemix_upwind_advection": false, "enable_tke": false, "c_k": 0.1, "c_eps": 0.7, "alpha_tke": 1.0, "mxl_min": 1e-12, "kappaM_min": 0.0, "kappaM_max": 100.0, "tke_mxl_choice": 1, "enable_tke_superbee_advection": false, "enable_tke_upwind_advection": false, "enable_tke_hor_diffusion": false, "K_h_tke": 2000.0, "enable_eke": false, "eke_lmin": 100.0, "eke_c_k": 1.0, "eke_cross": 1.0, "eke_crhin": 1.0, "eke_c_eps": 1.0, "eke_k_max": 10000.0, "alpha_eke": 1.0, "enable_eke_superbee_advection": false, "enable_eke_upwind_advection": false, "enable_eke_isopycnal_diffusion": false, "restart_input_filename": null, "restart_output_filename": "{identifier}_{itt:0>4d}.restart.h5", "restart_frequency": 311040000.0, "kappaH_min": 0.0, "enable_kappaH_profile": false, "enable_Prandtl_tke": true, "Prandtl_tke0": 10.0}

  setup_file :

  /Users/jklymak/Dropbox/Teaching/Eos431Phy441/source/Assign4/Overturning.py

  setup_code :

  #!/usr/bin/env python """ """ __VEROS_VERSION__ = '1.5.0' if __name__ == "__main__": raise RuntimeError( "Veros setups cannot be executed directly. " f"Try `veros run {__file__}` instead." ) # -- end of auto-generated header, original file below -- from veros import VerosSetup, veros_routine, veros_kernel, KernelOutput from veros.variables import allocate, Variable from veros.distributed import global_min, global_max from veros.core.operators import numpy as npx, update, at import matplotlib.pyplot as plt maxx = 0 maxy = 0 maxdepth = 4000 class BasinSetup(VerosSetup): """A model using cartesian coordinates in basin no beta effect Wind forcing over a channel; no buoyancy """ @veros_routine def set_parameter(self, state): settings = state.settings settings.identifier = "OverturningBeta" settings.nx, settings.ny, settings.nz = 12, 18, 30 settings.dt_mom = 4800 / 2 settings.dt_tracer = 86400 settings.runlen = 3600 * 24 * 30 * 12 * 1000 settings.x_origin = 0.0 settings.y_origin = 0.0 settings.coord_degree = True settings.enable_cyclic_x = False settings.enable_neutral_diffusion = True settings.K_iso_0 = 1000.0 settings.K_iso_steep = 500.0 settings.iso_dslope = 0.005 settings.iso_slopec = 0.01 settings.enable_skew_diffusion = True settings.enable_hor_friction = True # settings.enable_hor_diffusion = True dx = settings.degtom * 4 settings.A_h = (dx) ** 3 * 2e-11 settings.K_h = settings.A_h settings.enable_hor_friction_cos_scaling = True settings.hor_friction_cosPower = 1 settings.A_hbi = (3.5e5) ** 3 * 2e-11 * (3.5e5)**2 / 8 settings.K_hbi = (3.5e5) ** 3 * 2e-11 * (3.5e5)**2 / 8 # check dt: dt = dx**4 / 32 / settings.A_hbi print('dt', dt, 'dt_mom', settings.dt_mom) # assert dt > 2 * settings.dt_mom dt = dx**2 / 4 / settings.K_h print('dt', dt, 'dt_tracer', settings.dt_tracer) # assert dt > 2 * settings.dt_tracer settings.enable_superbee_advection = False #settings.enable_quadratic_bottom_friction = True #settings.r_quad_bot = 1e-1 settings.enable_implicit_vert_friction = True if False: settings.enable_tke = True settings.c_k = 0.1 settings.c_eps = 0.7 settings.alpha_tke = 30.0 settings.mxl_min = 1e-8 settings.tke_mxl_choice = 2 settings.kappaM_min = 1e-5 settings.kappaH_min = 1e-5 settings.enable_kappaH_profile = True settings.K_gm_0 = 1000.0 settings.enable_eke = True settings.eke_k_max = 1e4 settings.eke_c_k = 0.4 settings.eke_c_eps = 0.5 settings.eke_cross = 2.0 settings.eke_crhin = 1.0 settings.eke_lmin = 100.0 settings.enable_eke_superbee_advection = True settings.enable_eke_isopycnal_diffusion = True settings.kappaM_0 = 1e-4 settings.kappaH_0 = 1e-4 settings.enable_streamfunction = False settings.enable_idemix = False settings.enable_noslip_lateral = False settings.eq_of_state_type = 1 settings.restart_frequency = 3600 * 24 * 30 * 12 * 10 @veros_routine def set_grid(self, state): global maxx global maxy global maxdepth vs = state.variables settings = state.settings # ddz = 1000 / settings.nz ddz = npx.exp(-npx.arange(settings.nz) / 14) ddz = ddz / npx.sum(ddz) * maxdepth ddx = 4 # degrees maxx = settings.nx * ddx maxy = settings.ny * ddx vs.dxt = update(vs.dxt, at[...], ddx) vs.dyt = update(vs.dyt, at[...], ddx) vs.dzt = update(vs.dzt, at[...], ddz) @veros_routine def set_coriolis(self, state): vs = state.variables settings = state.settings vs.coriolis_t = update( vs.coriolis_t, at[:, :], 2 * settings.omega * npx.sin(vs.yt[None, :] / 180.0 * settings.pi) ) @veros_routine def set_topography(self, state): vs = state.variables settings = state.settings depths = -maxdepth * npx.ones((settings.nx, settings.ny)) if False: global maxx global maxy print('MAX', maxx) depths = (npx.interp(vs.xt[2:-2], [0, 150e-10, maxx-150e-10, maxx], [0, 1, 1, 0])[:, npx.newaxis] * npx.interp(vs.yt[2:-2], [0, 150e-10, maxy-150e-10, maxy], [0, 1, 1, 0])[npx.newaxis, :]) depths = - depths * 1000 depth_levels = 1 + npx.argmin(npx.abs(depths[:, :, npx.newaxis] - vs.zt[npx.newaxis, npx.newaxis, :]), axis=2) vs.kbot = update(vs.kbot, at[2:-2, 2:-2], npx.where(depth_levels > 0.0, depth_levels, 0)) vs.kbot = npx.where(vs.kbot < settings.nz, vs.kbot, 0) # depth of deepest cell. 0 is land, 1 is deepest # nz is shallowest.. if False: dd = npx.array([settings.nz+1, 1, 1, settings.nz+1]) yy = npx.array([0, 500, 20000-500, 20000]) k = npx.interp(vs.yt, yy, dd) k = npx.where(k>settings.nz, 0, k) # land k = npx.floor(k) # vs.kbot = update(vs.kbot, at[...], k) @veros_routine def set_initial_conditions(self, state): vs = state.variables settings = state.settings temp = npx.ones(settings.nz) * 4 # 4 degrees temp = npx.linspace(4, 4.04, settings.nz) vs.temp = update(vs.temp, at[...], temp[npx.newaxis, npx.newaxis, :, npx.newaxis]) vs.salt = update(vs.salt, at[...], 35) @veros_routine def set_forcing(self, state): vs = state.variables vs.update(set_forcing_kernel(state)) @veros_routine def set_diagnostics(self, state): settings = state.settings diagnostics = state.diagnostics diagnostics["snapshot"].output_frequency = 24 * 3600 * 30 * 12 # yearly diagnostics["overturning"].output_frequency = 24 * 3600 * 30 * 12 * 20 # 20 years diagnostics["overturning"].sampling_frequency = settings.dt_tracer * 10 @veros_routine def after_timestep(self, state): pass @veros_kernel def set_forcing_kernel(state): vs = state.variables global maxy sst = 16 - 12 * npx.cos((maxy - vs.yt) * npx.pi / maxy) sst = npx.clip(sst, 4, 30) # Model wants heat flux is in K m/s. Measured heat flux is order 100 W /m^2. # If dT = 10 C, Cp = 4000 J / C / kg, then wpiston = 2.5e-6 m/s # Bryan 87 uses 20 / 25 days = 9.3e-6 m/s. So lets use 10e-6 vs.forc_temp_surface = (sst - vs.temp[:, :, -1, vs.tau]) * vs.maskT[:, :, -1] * 100e-6 return KernelOutput( forc_temp_surface=vs.forc_temp_surface, )

The maximum temperature is very close to 28, and the minimum very close to 4 (though a bit warmer). The temperature is quite consistent across the basin, though close inspection reveals some small deviations.

Q2.2 Overturning as a function of mixing strength

There are 6 runs, each run at a different value of a constant vertical mixing co-efficient.

\(\kappa\ [m^2\,s^{-1}]\)	Runname
\(4\times10^{-6}\)	OverturningBeta4e6Coarse
\(1 \times 10^{-5}\)	OverturningBeta1e5Coarse
\(4\times 10^{-5}\)	OverturningBeta1e5Coarse
\(8\times 10^{-5}\)	OverturningBeta8e5Coarse
\(1\times 10^{-4}\)	OverturningBeta10e5Coarse
\(4\times 10^{-4}\)	OverturningBeta4e4Coarse

Plot the a cross-section at the end of the run from south to north with the mean north-south velocity and temprature contoured for at least three of these runs. Discuss the basic sense of the flow, and the differences observed as the mixing is increased.

fig, axs = plt.subplots(3, 1, layout='constrained', sharex=True, sharey=True)
for nn, td in enumerate(['OverturningBeta1e5Coarse', 'OverturningBeta4e5Coarse', 'OverturningBeta10e5Coarse']):
    with xr.open_dataset(td+'.snapshot.nc') as ds:
        ax = axs[nn]
        levels = np.arange(4.01, 30, 1)
        pc = ax.pcolormesh(ds.yu, ds.zt, ds.isel(Time=-1).v.integrate(coord='xt'), vmin=-2e-1, vmax=2e-1, cmap='RdBu_r')
        cl = ax.contour(ds.yu, ds.zt, ds.isel(Time=-1).temp.mean(dim='xt'), levels=levels, 
                        linewidths=0.4, colors='0.5')
        ax.clabel(cl, levels[::3], fontsize='small')
        ax.set_title(td, fontsize='small')
        ax.set_ylim(-4000, 0)
fig.colorbar(pc, ax=axs, shrink=0.6, extend='both', label='Mean velocity [m/s]')
ax.set_xlabel('Latitude')
ax.set_ylabel('Depth [m]')

/var/folders/p1/6grcm4fx3tx_dyzph3cbwk4c0000gn/T/ipykernel_10318/27345152.py:3: FutureWarning: In a future version of xarray decode_timedelta will default to False rather than None. To silence this warning, set decode_timedelta to True, False, or a 'CFTimedeltaCoder' instance.
  with xr.open_dataset(td+'.snapshot.nc') as ds:
/var/folders/p1/6grcm4fx3tx_dyzph3cbwk4c0000gn/T/ipykernel_10318/27345152.py:3: FutureWarning: In a future version of xarray decode_timedelta will default to False rather than None. To silence this warning, set decode_timedelta to True, False, or a 'CFTimedeltaCoder' instance.
  with xr.open_dataset(td+'.snapshot.nc') as ds:
/var/folders/p1/6grcm4fx3tx_dyzph3cbwk4c0000gn/T/ipykernel_10318/27345152.py:3: FutureWarning: In a future version of xarray decode_timedelta will default to False rather than None. To silence this warning, set decode_timedelta to True, False, or a 'CFTimedeltaCoder' instance.
  with xr.open_dataset(td+'.snapshot.nc') as ds:

Text(0, 0.5, 'Depth [m]')

fig, axs = plt.subplots(1, 3, layout='constrained', sharex=True, sharey=True, figsize=(8,4))
for nn, td in enumerate(['OverturningBeta1e5Coarse', 'OverturningBeta4e5Coarse', 'OverturningBeta10e5Coarse']):
    with xr.open_dataset(td+'.snapshot.nc') as ds:
        ds = ds.isel(Time=-1).sel(zt=-2000, method='nearest')
        ax = axs[nn]
        levels = np.arange(4.01, 30, 1)
        pc = ax.pcolormesh(ds.xt, ds.yu, ds.v, vmin=-2e-3, vmax=2e-3, cmap='RdBu_r')
        ax.clabel(cl, levels[::3], fontsize='small')
        ax.set_title(td, fontsize='small')
        ax.set_aspect(1)
        #ax.set_ylim(-4000, 0)
fig.colorbar(pc, ax=axs, shrink=0.6, extend='both', label='Mean velocity [m/s]')
ax.set_xlabel('Latitude')
ax.set_ylabel('Depth [m]')

/var/folders/p1/6grcm4fx3tx_dyzph3cbwk4c0000gn/T/ipykernel_10318/2472369364.py:3: FutureWarning: In a future version of xarray decode_timedelta will default to False rather than None. To silence this warning, set decode_timedelta to True, False, or a 'CFTimedeltaCoder' instance.
  with xr.open_dataset(td+'.snapshot.nc') as ds:
/var/folders/p1/6grcm4fx3tx_dyzph3cbwk4c0000gn/T/ipykernel_10318/2472369364.py:3: FutureWarning: In a future version of xarray decode_timedelta will default to False rather than None. To silence this warning, set decode_timedelta to True, False, or a 'CFTimedeltaCoder' instance.
  with xr.open_dataset(td+'.snapshot.nc') as ds:
/var/folders/p1/6grcm4fx3tx_dyzph3cbwk4c0000gn/T/ipykernel_10318/2472369364.py:3: FutureWarning: In a future version of xarray decode_timedelta will default to False rather than None. To silence this warning, set decode_timedelta to True, False, or a 'CFTimedeltaCoder' instance.
  with xr.open_dataset(td+'.snapshot.nc') as ds:

Text(0, 0.5, 'Depth [m]')

So, in general there is a poleward flow at the surface as warm water moves north to replace sinking cold water. There is a weaker return flow at depth to balance it.

As the mixing gets stronger the lateral gradients increase, and the warm water is mixed deeper. This drives a stronger overturning strength.

Q2.3 Overturning as function of time

We can calculate the strength of the overturning circulation ourselves. However, this is done so often in physical oceanography that many models have an overturning diagnostic.

Using the files like OverturningBeta4e6Coarse.overturning.nc, plot the minimum of the overturning strength south of 60 N (the reuslts are less clear north of 60 N) as a function of time. The correct field in these files to look at is vsf_depth.

How long do these runs take to reach “steady state”? Which runs are the fastest to get to “steady state”?

fig, ax = plt.subplots(1, layout='constrained', sharex=True, sharey=True)
for nn, td in enumerate([ 'OverturningBeta4e6Coarse', 'OverturningBeta1e5Coarse', 
                         'OverturningBeta4e5Coarse', 'OverturningBeta8e5Coarse',
                         'OverturningBeta10e5Coarse', 'OverturningBeta4e4Coarse', ]):
    with xr.open_dataset(td+'.overturning.nc', decode_timedelta=False) as ds:
        #display(ds)
        ds['Time'] = ds.Time / 360. # model years
        ds['vsf_depth'] = ds.vsf_depth / 1e6
        ds.vsf_depth.sel(yu=slice(20, 60)).min(dim=("yu", "zw")).plot(x='Time', ax=ax, label=td)
        ax.set_ylabel('Overturning [Sv]')
        ax.set_xlabel('Time [years]')
ax.legend()

<matplotlib.legend.Legend at 0x1149fda90>

The weaker mixing runs take longer to reach steady state, on the order of 1500 years.

Q2.4 Hosing experiment and stop…

Thre more runs were made, based on OverturningBeta4e5Coarse. The first was Hosing, which adds fresh water north of 50 N. This is supposed to be like glaciers and sea ice melting. The rate was approximately 0.1 Sv of fresh water added to the upper ocean.

The next two were StopHosing and StopHosing3. These stop the fresh water flux (all the glaciers have melted. The first stops after 100 years, and the second stops after 300 years.

As above, plot a time series of the strength of the overturning circulation for all 4 of these runs versus time. I found it really useful to also plot the difference between the minimum and maximum salinities in the model.

Describe what happens to the overturning circulation in the three new runs.

fig, axs = plt.subplots(2, 1, layout='constrained', sharex=True, sharey=False)
for nn, td in enumerate(['OverturningBeta4e5Coarse', 'Hosing', 'StopHosing', 'StopHosing3']):
    with xr.open_dataset(td+'.overturning.nc', decode_timedelta=False) as ds, xr.open_dataset(td+'.snapshot.nc', decode_timedelta=False) as ds0:
        ds['Time'] = ds.Time / 360. # model years
        ds0['Time'] = np.arange(0, len(ds0.Time)*10, 10) 
        ds0['Time'] = ds0['Time'] + ds['Time'][0]
        ds0 = ds0.isel(xt=4, xu=4)
        ds['vsf_depth'] = ds.vsf_depth / 1e6
        ds.vsf_depth.sel(yu=slice(20, 60), zw=slice(None, -100)).min(dim=("yu", "zw")).plot(x='Time', ax=axs[0], label=td)
        axs[0].set_ylabel('Overturning [Sv]')
        axs[0].set_xlabel('')
        
        dS = (ds0.salt.sel(yt=slice(20, 60)).max(dim=("yt", "zt")) - 
        ds0.salt.sel(yt=slice(20, 60)).min(dim=("yt", "zt")))
        axs[1].plot(dS.Time, dS, label=td)
        axs[1].set_ylabel('Max Salinity - Min salinity [g/kg]')
        axs[1].set_title('')
        ax.set_xlabel('Time [years]')
axs[1].legend(fontsize='small')

<matplotlib.legend.Legend at 0x116592fd0>

As the fresh water is added to the model, the overturning circulation rapidly stops; within 10 years or so. When the fresh water forcing is stopped, the overturning circulation recovers, but takes longer when more fresh water is added.

fig, ax = plt.subplots()
dt = np.diff(ds0.Time)
#ax.hist(dt, np.arange(-100.05, 100, 0.1))
ax.plot(ds0.Time)

[<matplotlib.lines.Line2D at 0x116795950>]