ECCOv4 Global Volume Budget Closure

ECCOv4 Global Volume Budget Closure#

Here we demonstrate the closure of volume budgets in ECCOv4 configurations. This notebook is draws heavily from evaluating_budgets_in_eccov4r3.pdf which explains the procedure with Matlab code examples by Christopher G. Piecuch. See ECCO Version 4 release documents: /doc/evaluating_budgets_in_eccov4r3.pdf).

Objectives#

Illustrate how volume budgets are closed globally.

Introduction#

ECCOv4 uses the $z^*$ coordinate system in which the depth of the vertical coordinate, $z^*$ varies with time as:

\[ z^* = \frac{z - \eta(x,y,t)}{H(x,y) + \eta(x,y,t)} H(x,y) \]

With $H$ being the model depth, $\eta$ being the model sea level anomaly, and $z$ being depth.

If the vertical coordinate didn’t change through time then volume fluxes across the ‘u’ and ‘v’ grid cell faces of a tracer cell could be calculated by multiplying the velocities at the face with the face area:

volume flux across ‘u’ face in the +x direction = $\mathit{UVEL}(x,y,k) \times \mathit{drF}(k) \times \mathit{dyG}(x,y) \times \mathit{hFacW}(x,y,k)$

volume flux across ‘v’ face in the +y direction = $\mathit{VVEL}(x,y,k) \times \mathit{drF}(k) \times \mathit{dxG}(x,y) \times \mathit{hFacS}(x,y,k)$

With dyG and dxG being the lengths of the ‘u’ and ‘v’ faces, drF being the grid cell height and hFacW and hFacS being the vertical fractions of the ‘u’ and ‘v’ grid cell faces that are open water (ECCOv4 uses partial cells to better represent bathymetry which can allows 0 < hfac $\le$ 1).

However, because the vertical coordiate varies with time in the $z^*$ system, the grid cell height drF varies with time as drF$\times s^*(t)$, with

\[ s^*(x,y,k,t) = 1 + \frac{\eta(x,y,t)}{H} \]

with $s^* > 1$ when $\eta > 0$

Thus, to calculate the volume fluxes grid cell through horizontal faces we must account for the time-varying grid cell face areas:

volume flux across ‘u’ in the +x direction face with $z^*$ coordinates = $\mathit{UVEL}(x,y,k) \times \mathit{drF}(k) \times \mathit{dyG}(x,y) \times \mathit{hFacW}(x,y,k) \times s^*(x,y,k,t)$

volume flux across ‘v’ in the +y direction face with $z^*$ coordinates = $\mathit{VVEL}(x,y,k) \times \mathit{drF}(k) \times \mathit{dxG}(x,y) \times \mathit{hFacS}(x,y,k) \times s^*(x,y,k,t)$

To make budget calculations easier we provide the scaled velocities quantities UVELMASS and VVELMASS,

\[ \mathit{UVELMASS}(x,y,k) = \mathit{UVEL}(x,y,k) \times \mathit{hFacW}(x,y,k) \times s^{*}(x,y,k,t) \]

and $$ \mathit{VVELMASS}(x,y,k) = \mathit{VVEL}(x,y,k) \times \mathit{hFacS}(x,y,k) \times s^{*}(x,y,k,t) $$

It is worth noting that the word mass in UVELMASS and VVELMASS is confusing since there is no mass involved here. Think of these terms as simply being UVEL and VVEL multiplied by the fraction of the grid cell height that is open grid cell face across which the volume transport occurs. Partial cell bathymetry can make this fraction (hFacW, hFacS) less than one, and the $s^*$ scaling factor further adjusts this fraction higher or lower through time.

Fully closing the budget requires the vertical volume fluxes across the top and bottom ‘w’ faces of the grid cell and surface freshwater fluxes. Regarding vertical volume fluxes, there are no $s^*$ or hFac equivalent scaling factors that modify our top and bottom grid cell areas. Therefore, vertical volume fluxes through ‘w’ faces are simply:

volume flux across ‘w’ face in the +z direction = $\mathit{WVEL}(x,y,k) \times \mathit{rA}(x,y)$

Note: Inexplicably, the term WVEL is provided with the silly name WVELMASS. Sometimes it’s difficult to ignore other people’s poor life choices, but please try to do so here. Ignore the confusing name, WVELMASS is identical to WVEL.

In the $z^*$ coordinate system the depth of the surface grid cell is always $z^* = 0$. In the MITgcm, WELMASS at the top of the surface grid cell is the liquid volume flux out of the ocean surface and is proportional to the vertical ocean mass flux, oceFWflx

`ETAN` in a Boussinesq Model#

ECCOv4 uses a volume-conserving Boussinesq formulation of the MITgcm. Because volume is conserved in Boussinesq formulations, seawater density changes do not change model sea level anomaly, ETAN. The following demonstration ofETAN budget closure considers volumetric fluxes but does not take into consideration expansion/contraction due to changes in density. Furthermore, model ETAN changes with the exchange of water between ocean and sea-ice while in reality the Archimedes principle holds that sea level should not change following the growth or melting of sea-ice because floating sea-ice displaces a volume of seawater equal to its weight. Thus, ETAN is not comparable to observed sea level.

We correct ETAN to make a sea surface height field that is comaprable with observations by making three corrections: with a) the “Greatbatch correction”, a time varying, globally-uniform correction to ocean volume due to changes in global mean density, b) the inverted barometer (IB) correction (see SSHIBC) and c) the ‘sea ice load’ correction to account for the displacement of seawater due to submerged sea-ice and snow (see sIceLoad). A demonstration of these corrections is outside the scope of this tutorial. Here we focus on closing the model budget keeping in mind that we are neglecting sea level changes from changes in global mean density and the fact that ETAN does not account for volume displacement due to submerged sea-ice.

Greatbatch, 1994. J. of Geophys. Res. Oceans, https://doi.org/10.1029/94JC00847

Calculate LHS: $\eta$ time tendency: $G_{\text{total tendency}}$#

We calculate the monthly-averaged time tendency of ETAN by differencing monthly ETAN snapshots. Subtract the numpy arrays $\eta(t+1)$ - $\eta(t)$. This operation gives us $\Delta$ ETAN $/ \Delta$ t (month) records.

num_months = len(ecco_monthly_snaps.time)
G_total_tendency_month = \
    ecco_monthly_snaps.ETAN.isel(time=range(1,num_months)).values - \
    ecco_monthly_snaps.ETAN.isel(time=range(0,num_months-1)).values

# The result is a numpy array of 264 months
print('shape of G_total_tendency_month: ', G_total_tendency_month.shape)

shape of G_total_tendency_month:  (288, 13, 90, 90)

ecco_monthly_mean.oceFWflx.shape

(288, 13, 90, 90)

# Convert this numpy array to an xarray Dataarray to take advantage 
# of xarray time indexing.
# The easiest way is to copy an existing DataArray that has the 
# dimensions and time indexes that we want,
# replace its values, and change its name.  

tmp = ecco_monthly_mean.oceFWflx.copy(deep=True)
tmp.values = G_total_tendency_month
tmp.name = 'G_total_tendency_month'
G_total_tendency_month = tmp

# the nice thing is that now the time values of G_total_tendency_month now line 
# up with the time values of the time-mean fields (middle of the month)
print('\ntime of first array of G_total_tendency_month');
print(G_total_tendency_month.time[0].values)

time of first array of G_total_tendency_month
1993-01-16T12:00:00.000000000

Now convert $\Delta$ ETAN $/ \Delta$ t (month) to $\Delta$ ETAN $/ \Delta$ t (seconds) by dividing by the number of seconds in each month.

secs_per_month = (ecco_monthly_snaps.time[1:].values - \
ecco_monthly_snaps.time[0:-1].values).astype('float')/(1.e9)   #divide by 1.e9 to convert ns to seconds


# Make a DataArray with the number of seconds in each month, 
#time indexed to the times in dETAN_dT_perMonth (middle of each month)
secs_per_month = xr.DataArray(secs_per_month, \
                              coords={'time': G_total_tendency_month.time.values}, \
                              dims='time')

# show number of seconds in the first two months:
print('# of seconds in Jan and Feb 1993 ', secs_per_month[0:2].values)

# sanity check: show number of days in the first two months:
print('# of days in Jan and Feb 1993 ', secs_per_month[0:2].values/3600/24)

# of seconds in Jan and Feb 1993  [2678400. 2419200.]
# of days in Jan and Feb 1993  [31. 28.]

# convert dETAN_dT_perMonth to perSeconds
G_total_tendency = G_total_tendency_month / secs_per_month

Plot the time-mean $\partial \eta / \partial t$, total $\Delta \eta$, and one example $\partial \eta / \partial t$ field#

Time-mean $\partial \eta / \partial t$#

To calculate the time averaged $G_{total_tendency}$ one might be tempted to do the following

> G_total_tendency_mean = G_total_tendency.mean('time')

But that would be folly because the ‘mean’ function does not know that the number of days in each month is different! The result would downweight Februarys and upweight Julys. We have to weight the tendency records by the length of each month. A clever way of doing that is provided in the xarray documents: https://xarray-test.readthedocs.io/en/latest/examples/monthly-means.html

In our case we know the length of each month, we just calculated it above in secs_per_month. We will weight each month by the relative # of seconds in each month and sum to get a weighted average.

# the weights are just the # of seconds per month divided by total seconds
month_length_weights = secs_per_month / secs_per_month.sum()

The time mean of the ETAN tendency, $\overline{G_{\text{total tendency}}}$, is given by

$\overline{G_{\text{total tendency}}} = \sum_{i=1}^{nm} w_i G_{\text{total tendency}}$

with $\sum_{i=1}^{nm} w_i = 1$ and nm=number of months

# the weights sum to 1
print(month_length_weights.sum())

<xarray.DataArray ()> Size: 8B
array(1.)

# the weighted mean weights by the length of each month (in seconds)
G_total_tendency_mean = (G_total_tendency*month_length_weights).sum('time')

plt.figure(figsize=(20,8));
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC,\
                              G_total_tendency_mean,show_colorbar=True,\
                              cmin=-1e-9, cmax=1e-9, \
                              cmap='RdBu_r', user_lon_0=-67,\
                              dx=map_dx,dy=map_dy);
plt.title('Average $\partial \eta / \partial t$ [m/s]', fontsize=20);

../../_images/f0254b6fe9096708255fa6b79eb770c0e4b72eadd867a490ee36af511b32a8e4.png

Total $\Delta \eta$#

The time average eta tendency is small, about 1 billionth of a meter per second. The ECCO period is coming up to a billion seconds though… How much did ETAN change over the analysis period?

# the number of seconds in the entire period 
seconds_in_entire_period = \
    float(ecco_monthly_snaps.time[-1] - ecco_monthly_snaps.time[0])/1e9
print ('seconds in analysis period: ', seconds_in_entire_period)

# which is also the sum of the number of seconds in each month
print('sum of seconds in each month ', secs_per_month.sum().values)

seconds in analysis period:  757382400.0
sum of seconds in each month  757382400.0

ETAN_delta = G_total_tendency_mean*seconds_in_entire_period

plt.figure(figsize=(20,8));
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, \
                              ETAN_delta,show_colorbar=True,\
                              cmin=-1, cmax=1, \
                              cmap='RdBu_r', user_lon_0=-67, \
                              dx=map_dx,dy=map_dy);
plt.title('Predicted $\Delta \eta$ [m] from $\partial \eta / \partial t$', \
          fontsize=20);

../../_images/66402b77bf20a792be6a43de89268b3d31b396ac4af46a80856b2fae2f81066f.png

We can sanity check the total ETAN change that we found by multipling the time-mean ETAN tendency with the number of seconds in the simulation by comparing that with the difference in ETAN between the end of the last month and start of the first month.

ETAN_delta_method_2 = ecco_monthly_snaps.ETAN.isel(time=-1).values - \
    ecco_monthly_snaps.ETAN.isel(time=0).values

plt.figure(figsize=(20,8));
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, \
                              ETAN_delta_method_2,
                              show_colorbar=True,\
                              cmin=-1, cmax=1, \
                              cmap='RdBu_r', user_lon_0=-67,\
                              dx=map_dx,dy=map_dy);

plt.title('Actual $\Delta \eta$ [m]', fontsize=20);

../../_images/661bbf3d1dcdb748bd2104fbb0dfb6ec4f9f19494c48132bb8670c0388719dcf.png

plt.figure(figsize=(20,8));
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, \
                              ETAN_delta_method_2-ETAN_delta,\
                              show_colorbar=True,\
                              cmin=-1e-6, cmax=1e-6, \
                              cmap='RdBu_r', user_lon_0=-67,\
                              dx=map_dx,dy=map_dy);
plt.title('Difference between actual and predicted $\Delta \eta$ [m]', \
          fontsize=20);

../../_images/1494b24d733b0bb897988a7e30136aa3db5f8882c90876e951ce5e4d683dc37d.png

That’s a big woo, these are the same to within 10^-6 meters!

Example $\partial \eta / \partial t$ field#

plt.figure(figsize=(20,8));

curr_t_ind = 100
curr_time = G_total_tendency.time[curr_t_ind].values
print(curr_time)
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, \
                              G_total_tendency.isel(time=curr_t_ind),\
                              show_colorbar=True,\
                              cmin=-1e-7, cmax=1e-7,\
                              cmap='RdBu_r', user_lon_0=-67,\
                              dx=map_dx,dy=map_dy);

plt.title('$\partial \eta / \partial t$ [m/s] during ' + 
          str(curr_time)[0:4] +'/' + str(curr_time)[5:7], fontsize=20);

2001-05-16T12:00:00.000000000

../../_images/a9fd9f054a8cfc54fc006bf016acd09eb3ba8d4a6f9a91c827f97e7fe0618fc0.png

For any given month the time rate of change of ETAN is two orders of magnitude smaller than the 1993-2015 mean. In the above we are looking at May 2001. We see positive ETAN tendency due sea ice melting in the northern hemisphere (e.g., Baffin Bay, Greenland Sea, and Chukchi Sea).

Calculate RHS: $\eta$ tendency due to surface fluxes, $G_\text{surface fluxes}$#

Surface mass fluxes are given in oceFWflx. Convert surface mass flux to a vertical velocity by dividing by the reference density rhoConst= 1029 kg m-3

# tendency of eta implied by surface volume fluxes (m/s)
G_surf_fluxes = ecco_monthly_mean.oceFWflx/rhoconst

Plot the time-mean, total, and one month average of $G_{\text{surface fluxes}}$#

Time-mean $G_{\text{surface fluxes}}$#

We calculate the time-mean surface flux $\eta$ tendency using the same weights as the total $\eta$ tendency.

G_surf_fluxes_mean= (G_surf_fluxes*month_length_weights).sum('time')
G_surf_fluxes_mean.shape

(13, 90, 90)

plt.figure(figsize=(20,8));
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, \
                              G_surf_fluxes_mean,
                              show_colorbar=True,\
                              cmin=-1e-7, cmax=1e-7, \
                              cmap='RdBu_r', user_lon_0=-67,
                              dx=map_dx,dy=map_dy)
plt.title('Average $\partial \eta / \partial t$ [m/s] implied by oceFWflx surface mass fluxes\n Negative = Water out',
          fontsize=20);

../../_images/06f65641b00a6f0fac37e984559ee377e9c3e964a6dd70d7506e52ca0c2b6b85.png

Total $\Delta \eta$ due to surface fluxes#

If there were no other terms on the RHS to balance surface fluxes, the total change in ETAN between 1993 and 2017 would be order of 10s of meters almost everwhere.

ETAN_delta_surf_fluxes = G_surf_fluxes_mean*seconds_in_entire_period

plt.figure(figsize=(20,8));
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, \
                              ETAN_delta_surf_fluxes,show_colorbar=True,\
                              cmin=-100, cmax=100, \
                              cmap='RdBu_r', user_lon_0=-67, \
                              dx=map_dx,dy=map_dy);
plt.title('$\Delta \eta$ [m] implied by oceFWflx surface mass fluxes', 
          fontsize=20);

../../_images/964a04a2fbd5fbffa685772283fecf556c5f4c20935434cf3c5a9ddc5d7ef81a.png

Example $\partial \eta / \partial t$ impliedy by surface fluxes#

plt.figure(figsize=(20,8));

curr_t_ind = 100
curr_time = G_surf_fluxes.time[curr_t_ind].values
print(curr_time)
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, \
                              G_surf_fluxes.isel(time=curr_t_ind),show_colorbar=True,\
                             cmin=-1e-7, cmax=1e-7, \
                              cmap='RdBu_r', user_lon_0=-67,\
                              dx=map_dx,dy=map_dy);
plt.title('$\partial \eta / \partial t$ [m/s] implied by oceFWflx surface mass fluxes ' + 
          str(curr_time)[0:4] +'-' + str(curr_time)[5:7], fontsize=20);

2001-05-16T12:00:00.000000000

../../_images/6d60aaac434471a5047ae1a71844cf364110295c9752d9271cdf8ef4c8cb6cca.png

For any given month the time rate of change of ETAN is almost the same as its 22 year mean. Differences are largest in the high latitudes where sea-ice melt and growth during any particular month induce large changes in ETAN.

Calculate RHS: $\eta$ tendency due to volumetric flux divergence, $G_\text{volumetric fluxes}$#

First we will look at vertical volumetric flux divergence, then horizontal volumetric flux divergence.

Vertical volumetric flux divergence#

It turns out that WVELMASS at k_l=0 (the top face of the top tracer cell) is proportional to the ocean surface mass flux oceFWflx. The vertical velocity of the ocean surface is equal to the rate at which water is being added or removed across the top surface of the uppermost grid cell. This is demonstrated by differencing the velocity at the top ‘w’ face of the uppermost tracr cell WVELMASS (k_l = 0) and the velocity equivalent of transporting the surface mass flux term oceFWFlx through this same face.

First, find the time-mean vertical velocity at the liquid ocean surface

WVELMASS_surf_mean = \
    (ecco_monthly_mean.WVELMASS.isel(k_l=0)*month_length_weights).sum('time')

Next, find the time-mean vertical velocity implied by the oceFWflx at k_l=0:

WVEL_from_oceFWflx_mean = \
    -(ecco_monthly_mean.oceFWflx*month_length_weights).sum('time')/rhoconst

plt.figure(figsize=(15,15))
#plt.sca(axs[0,0])
F=ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, \
                                WVELMASS_surf_mean,\
                                show_colorbar=True,\
                                cmin=-1e-7, cmax=1e-7, \
                                cmap='RdBu_r', user_lon_0=-67,\
                                dx=2,dy=2, subplot_grid=[3,1,1]);
F[1].set_title('A: surface velocity from WVEL [m/s]')

F=ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, \
                                WVEL_from_oceFWflx_mean,\
                                show_colorbar=True,\
                                cmin=-1e-7, cmax=1e-7, \
                                cmap='RdBu_r', user_lon_0=-67,\
                                dx=2,dy=2, subplot_grid=[3,1,2])
F[1].set_title("B: surface velocity implied by oceFWflx [m/s]")

F=ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, \
                                WVELMASS_surf_mean-WVEL_from_oceFWflx_mean,\
                                show_colorbar=True,\
                                cmin=-1e-7, cmax=1e-7, \
                                cmap='RdBu_r', user_lon_0=-67,\
                                dx=2,dy=2, subplot_grid=[3,1,3])
F[1].set_title("difference between A and B");

../../_images/6ec9bd2948c06447f18897c6242c2616df57f4b861efac8bb8ec7be61b803c36.png

WVELMASS at the surface evidently is the same as the surface velocity implied by the surface mass flux oceFWflx. Thus, we do not actually need oceFWflx to close the volume budget. However, to keep the surface forcing term explicitly represented, we will keep oceFWflx and instead zero out the values of WVELMASS at the surface so as to avoid double counting.

Calculate the vertical volume fluxes at each level: (velocity x grid cell area) [m3 s-1]

vol_transport_z = ecco_monthly_mean.WVELMASS * ecco_grid.rA

Set the volume transport at the surface level to be zero because we already counted the fluxes out of the domain with oceFWflx. We will also set the volume transport in land areas (currently NaN values) to zero so the divergence adjacent to the bottom is computed correctly.

# vol_transport_z.isel(k_l=0).values[:] = 0
# vol_transport_z.values[ecco_grid.hFacC.values == 0] = 0

# ## For some reason the above commands don't seem to work with dask arrays, 
# ## probably due to some complication resetting the values of an existing dask array.
## However, using where works with dask arrays
vol_transport_z = vol_transport_z.where(np.expand_dims(ecco_grid.k_l.values,axis=(0,-1,-2,-3)) != 0,\
                                        0)
vol_transport_z = vol_transport_z.where(ecco_grid.hFacC.values > 0,0)

Each grid cell has a top and bottom surface and therefore WVELMASS should have 51 vertical levels (one more than the number of tracer cells). For some reason we only have 50 vertical levels, with the bottom of the 50th tracer cell missing. To calculate vertical flux divergence we need to add this 51st WVELMASS which is everywhere zero (no volume flux from the seafloor). The xgcm library handles this in its diff routine by specifying the boundary=’fill’ and fill_value = 0.

# re-chunk xarray DataArray for better performance
vol_transport_z = vol_transport_z.chunk({'time':12,'k_l':vol_transport_z.sizes['k_l'],'tile':1})

# this computation can take a little while, so we'll loop through 12 time indices (1 year)
# at a time and show when each computation has been completed

import time

G_vertical_flux_divergence_depth_integrated_time_mean = xr.DataArray(\
    data=np.zeros(vol_transport_z.shape[-3:]),dims=['tile','j','i'])
size_superchunk = 12
for time_superchunk in range(int(np.ceil(vol_transport_z.sizes['time']/size_superchunk))):
    time_log = [time.time()]
    
    curr_time_ind = np.arange(size_superchunk*time_superchunk,\
                              np.fmin(size_superchunk*(time_superchunk+1),\
                                      vol_transport_z.sizes['time']))
    
    
    # volume flux divergence into each grid cell, m^3 / s 
    vol_vert_divergence = ecco_xgcm_grid.diff(vol_transport_z.isel(time=curr_time_ind),axis='Z',\
                                              boundary='fill',fill_value=0).squeeze()
    
    # change in eta per unit time due to volumetric vertical convergence
    # at each depth level: m/s
    G_vertical_flux_divergence = vol_vert_divergence / ecco_grid.rA
    
    # change in eta per unit time due to vertical integral of 
    # volumetric horizonal convergence: m/s
    G_vertical_flux_divergence_depth_integrated = G_vertical_flux_divergence.sum('k')
        
    # Calculate the time-mean surface flux $\eta$ tendency using 
    # the same weights as the total $\eta$ tendency.
    G_vertical_flux_divergence_depth_integrated_time_mean.values += \
        (G_vertical_flux_divergence_depth_integrated * month_length_weights[curr_time_ind]).sum('time').values

    time_log = np.append(time_log,time.time())
    
    print(f'Completed calculation for time ind = {curr_time_ind[[0,-1]]}, '\
            +f'time elapsed = {np.diff(time_log)[0]} s')

Completed calculation for time ind = [ 0 11], time elapsed = 14.48197865486145 s
Completed calculation for time ind = [12 23], time elapsed = 2.9820005893707275 s
Completed calculation for time ind = [24 35], time elapsed = 2.9535562992095947 s
Completed calculation for time ind = [36 47], time elapsed = 2.851987838745117 s
Completed calculation for time ind = [48 59], time elapsed = 2.9231979846954346 s
Completed calculation for time ind = [60 71], time elapsed = 3.1725475788116455 s
Completed calculation for time ind = [72 83], time elapsed = 2.911944627761841 s
Completed calculation for time ind = [84 95], time elapsed = 2.9191393852233887 s
Completed calculation for time ind = [ 96 107], time elapsed = 3.018028974533081 s
Completed calculation for time ind = [108 119], time elapsed = 3.019850015640259 s
Completed calculation for time ind = [120 131], time elapsed = 2.9223062992095947 s
Completed calculation for time ind = [132 143], time elapsed = 3.1544744968414307 s
Completed calculation for time ind = [144 155], time elapsed = 2.9720146656036377 s
Completed calculation for time ind = [156 167], time elapsed = 2.958627462387085 s
Completed calculation for time ind = [168 179], time elapsed = 2.9597220420837402 s
Completed calculation for time ind = [180 191], time elapsed = 3.2003564834594727 s
Completed calculation for time ind = [192 203], time elapsed = 2.9414868354797363 s
Completed calculation for time ind = [204 215], time elapsed = 2.937130928039551 s
Completed calculation for time ind = [216 227], time elapsed = 2.928783893585205 s
Completed calculation for time ind = [228 239], time elapsed = 2.9334447383880615 s
Completed calculation for time ind = [240 251], time elapsed = 2.9076087474823 s
Completed calculation for time ind = [252 263], time elapsed = 2.955669641494751 s
Completed calculation for time ind = [264 275], time elapsed = 2.9211392402648926 s
Completed calculation for time ind = [276 287], time elapsed = 2.9150209426879883 s

Plot the time-mean $G_{\text{vertical flux divergence}}$#

plt.figure(figsize=(20,8));
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC,\
                              G_vertical_flux_divergence_depth_integrated_time_mean,\
                              show_colorbar=True,\
                              cmap='RdBu_r', \
                              cmin=-1e-9, cmax=1e-9,
                              dx=map_dx,dy=map_dy);
plt.title('Average $\partial \eta / \partial t$ [m/s] due to vertical flux divergence\n Negative = Water out',
          fontsize=20);

../../_images/b9f821a1b91fa0b9179510151f255ec15ce7b0bc77c29d8ad972a2c3fdc0ba65.png

These values are everywhere essentially zero (numerical noise). On average, the only vertical flux divergence in the column is across the ocean surface. Below the surface, the sum of vertical flux divergence in all tracer cells in the column must be zero because any divergence in any one particular cell is exactly offset by convergence in another cell. Net convergence into the column manifests as a positive vertical velocity at the surface which is equal to oceFWflux in the time-mean. Thanks to Hong Zhang for comments that improved this explanation.

Horizontal Volume Flux Divergence#

# Volumetric transports in x and y(m^3/s)
vol_transport_x = ecco_monthly_mean.UVELMASS * ecco_grid.dyG * ecco_grid.drF
vol_transport_y = ecco_monthly_mean.VVELMASS * ecco_grid.dxG * ecco_grid.drF

# Set fluxes on land to zero (instead of NaN)
vol_transport_x = vol_transport_x.where(ecco_grid.hFacW.values > 0,0)
vol_transport_y = vol_transport_y.where(ecco_grid.hFacS.values > 0,0)

# re-chunk arrays for better performance
vol_transport_x = vol_transport_x.chunk({'time':1,'k':-1,'tile':-1,'j':-1,'i_g':-1})
vol_transport_y = vol_transport_y.chunk({'time':1,'k':-1,'tile':-1,'j_g':-1,'i':-1})

def diff_2d_flux_llc90(flux_vector_dict):
    """
    A function that differences flux variables on the llc90 grid.
    Can be used in place of xgcm's diff_2d_vector.
    """

    u_flux = flux_vector_dict['X']
    v_flux = flux_vector_dict['Y']
    
    u_flux_padded = u_flux.pad(pad_width={'i_g':(0,1)},mode='constant',constant_values=np.nan)\
                            .chunk({'i_g':u_flux.sizes['i_g']+1})
    v_flux_padded = v_flux.pad(pad_width={'j_g':(0,1)},mode='constant',constant_values=np.nan)\
                            .chunk({'j_g':v_flux.sizes['j_g']+1})
    
    def da_replace_at_indices(da,indexing_dict,replace_values):
        # replace values in xarray DataArray using locations specified by indexing_dict
        array_data = da.data
        indexing_dict_bynum = {}
        for axis,dim in enumerate(da.dims):
            if dim in indexing_dict.keys():
                indexing_dict_bynum = {**indexing_dict_bynum,**{axis:indexing_dict[dim]}}
        ndims = len(array_data.shape)
        indexing_list = [':']*ndims
        for axis in indexing_dict_bynum.keys():
            indexing_list[axis] = indexing_dict_bynum[axis]
        indexing_str = ",".join(indexing_list)

        # using exec isn't ideal, but this works for both NumPy and Dask arrays
        exec('array_data['+indexing_str+'] = replace_values')        
        
        return da
    
    # u flux padding
    for tile in range(0,3):
        u_flux_padded = da_replace_at_indices(u_flux_padded,{'tile':str(tile),'i_g':'-1'},\
                                              u_flux.isel(tile=tile+3,i_g=0).data)
    for tile in range(3,6):
        u_flux_padded = da_replace_at_indices(u_flux_padded,{'tile':str(tile),'i_g':'-1'},\
                                              v_flux.isel(tile=12-tile,j_g=0,i=slice(None,None,-1)).data)
    u_flux_padded = da_replace_at_indices(u_flux_padded,{'tile':'6','i_g':'-1'},\
                                          u_flux.isel(tile=7,i_g=0).data)
    for tile in range(7,9):
        u_flux_padded = da_replace_at_indices(u_flux_padded,{'tile':str(tile),'i_g':'-1'},\
                                              u_flux.isel(tile=tile+1,i_g=0).data)
    for tile in range(10,12):
        u_flux_padded = da_replace_at_indices(u_flux_padded,{'tile':str(tile),'i_g':'-1'},\
                                              u_flux.isel(tile=tile+1,i_g=0).data)
        
    # v flux padding
    for tile in range(0,2):
        v_flux_padded = da_replace_at_indices(v_flux_padded,{'tile':str(tile),'j_g':'-1'},\
                                              v_flux.isel(tile=tile+1,j_g=0).data)
    v_flux_padded = da_replace_at_indices(v_flux_padded,{'tile':'2','j_g':'-1'},\
                                          u_flux.isel(tile=6,j=slice(None,None,-1),i_g=0).data)
    for tile in range(3,6):
        v_flux_padded = da_replace_at_indices(v_flux_padded,{'tile':str(tile),'j_g':'-1'},\
                                              v_flux.isel(tile=tile+1,j_g=0).data)
    v_flux_padded = da_replace_at_indices(v_flux_padded,{'tile':'6','j_g':'-1'},\
                                          u_flux.isel(tile=10,j=slice(None,None,-1),i_g=0).data)
    for tile in range(7,10):
        v_flux_padded = da_replace_at_indices(v_flux_padded,{'tile':str(tile),'j_g':'-1'},\
                                              v_flux.isel(tile=tile+3,j_g=0).data)
    for tile in range(10,13):
        v_flux_padded = da_replace_at_indices(v_flux_padded,{'tile':str(tile),'j_g':'-1'},\
                                              u_flux.isel(tile=12-tile,j=slice(None,None,-1),i_g=0).data)
    
    # take differences
    diff_u_flux = u_flux_padded.diff('i_g')
    diff_v_flux = v_flux_padded.diff('j_g')
    
    # include coordinates of input DataArrays and correct dimension/coordinate names
    diff_u_flux = diff_u_flux.assign_coords(u_flux.coords).rename({'i_g':'i'})
    diff_v_flux = diff_v_flux.assign_coords(v_flux.coords).rename({'j_g':'j'})

    diff_flux_vector_dict = {'X':diff_u_flux,'Y':diff_v_flux}
    
    return diff_flux_vector_dict

# # This way of computing horizontal volume divergence has more concise code, 
# # but it may not run optimally in some computing environments.
# # In that case, try doing the computation in "superchunks" of 12 months (1 year)
# # at a time, as coded in the next cell.


# # Difference of horizontal transports in x and y directions
# vol_flux_diff = ecco_xgcm_grid.diff_2d_vector({'X': vol_transport_x, \
#                                                'Y': vol_transport_y},\
#                                               boundary='fill')
# 
# # volume flux divergence into each grid cell, m^3 / s 
# vol_horiz_divergence = (vol_flux_diff['X'] + vol_flux_diff['Y'])
# # restore time coordinate to DataArray (lost in xgcm.diff_2d_vector operation)
# vol_horiz_divergence = vol_horiz_divergence.assign_coords({'time':vol_transport_x.time.data})
# 
# # change in eta per unit time due to volumetric horizonal
# # convergence at each depth level: m/s
# # a positive DIVERGENCE leads to negative eta tendency
# G_vol_horiz_divergence = -vol_horiz_divergence / ecco_grid.rA
# 
# # change in eta in each grid cell per unit time due to horiz. divergence: m/s
# G_vol_horiz_divergence_depth_integrated = G_vol_horiz_divergence.sum('k').compute()
# 
# # calculate time-mean using the month length weights
# G_vol_horiz_divergence_depth_integrated_mean = \
#     (G_vol_horiz_divergence_depth_integrated * month_length_weights).sum('time')

# complete horizontal divergence calculation, 12 time indices (1 year) at a time

import copy

area_masked = ecco_grid.rA.where(ecco_grid.hFacC.isel(k=0) > 0).compute()
area_masked_sum = area_masked.sum(dim=('i','j','tile')).values

G_vol_horiz_divergence_depth_integrated_globalmean_array = np.empty((vol_transport_x.sizes['time'],))
G_vol_horiz_divergence_depth_integrated_mean = xr.DataArray(\
    data=np.zeros(vol_transport_x.shape[-3:]),dims=['tile','j','i'])
size_superchunk = 12
for time_superchunk in range(int(np.ceil(vol_transport_x.sizes['time']/size_superchunk))):
    time_log = np.array([time.time()])
    
    curr_time_ind = np.arange(size_superchunk*time_superchunk,\
                              np.fmin(size_superchunk*(time_superchunk+1),\
                                      vol_transport_x.sizes['time']))    
    
    # # Difference of horizontal transports in x and y directions
    # # (can use either xgcm's diff_2d_vector or the diff_2d_flux_llc90
    # # function defined above)
    # vol_flux_diff = ecco_xgcm_grid.diff_2d_vector({'X': vol_transport_x.isel(time=curr_time_ind), \
    #                                                'Y': vol_transport_y.isel(time=curr_time_ind)},\
    #                                               boundary='fill')
    vol_flux_diff = diff_2d_flux_llc90({'X': vol_transport_x.isel(time=curr_time_ind),\
                                         'Y': vol_transport_y.isel(time=curr_time_ind)})
    
    # volume flux divergence into each grid cell, m^3 / s 
    vol_horiz_divergence = (vol_flux_diff['X'] + vol_flux_diff['Y'])
    
    # restore time coordinate to DataArray (lost in xgcm.diff_2d_vector operation)
    vol_horiz_divergence = vol_horiz_divergence.assign_coords({'time':vol_transport_x.time[curr_time_ind].data})
    
    # change in eta per unit time due to volumetric horizonal
    # convergence at each depth level: m/s
    # a positive DIVERGENCE leads to negative eta tendency
    curr_G_vol_horiz_divergence = -vol_horiz_divergence / ecco_grid.rA
    if curr_time_ind[0] == 0:
        G_vol_horiz_divergence = copy.deepcopy(curr_G_vol_horiz_divergence)
    else:
        G_vol_horiz_divergence = xr.concat((G_vol_horiz_divergence,\
                                            curr_G_vol_horiz_divergence),\
                                           dim='time',\
                                           coords='minimal',compat='override')
    
    # change in eta in each grid cell per unit time due to horiz. divergence: m/s
    curr_G_vol_horiz_divergence_depth_integrated = xr.DataArray(\
        data=curr_G_vol_horiz_divergence.sum('k'),dims=['time','tile','j','i'])
    if curr_time_ind[0] == 0:
        G_vol_horiz_divergence_depth_integrated = copy.deepcopy(curr_G_vol_horiz_divergence_depth_integrated)
    else:
        G_vol_horiz_divergence_depth_integrated = xr.concat((G_vol_horiz_divergence_depth_integrated,\
                                                             curr_G_vol_horiz_divergence_depth_integrated),\
                                                            dim='time',\
                                                            coords='minimal',compat='override')

    # compute global mean at each time using area weighting
    G_vol_horiz_divergence_depth_integrated_globalmean_array[curr_time_ind] = \
        (area_masked*curr_G_vol_horiz_divergence_depth_integrated).sum(dim=('i','j','tile')).values/area_masked_sum
    
    # calculate time-mean using the month length weights
    G_vol_horiz_divergence_depth_integrated_mean += \
        (curr_G_vol_horiz_divergence_depth_integrated*month_length_weights[curr_time_ind]).sum('time').values
    
    time_log = np.append(time_log,time.time())
    
    print(f'Completed calculation for time ind = {curr_time_ind[[0,-1]]}, '\
            +f'time elapsed = {np.diff(time_log)[0]} s')

# # create DataArray of global mean values
G_vol_horiz_divergence_depth_integrated_globalmean = xr.DataArray(\
    data=G_vol_horiz_divergence_depth_integrated_globalmean_array,dims=['time'],coords={'time':vol_transport_x.time})

# # calculate time-mean using the month length weights
# G_vol_horiz_divergence_depth_integrated_mean = (G_vol_horiz_divergence_depth_integrated*month_length_weights).sum('time')

Completed calculation for time ind = [ 0 11], time elapsed = 5.920902490615845 s
Completed calculation for time ind = [12 23], time elapsed = 5.610679864883423 s
Completed calculation for time ind = [24 35], time elapsed = 5.947405576705933 s
Completed calculation for time ind = [36 47], time elapsed = 5.594774961471558 s
Completed calculation for time ind = [48 59], time elapsed = 5.660876989364624 s
Completed calculation for time ind = [60 71], time elapsed = 5.619669675827026 s
Completed calculation for time ind = [72 83], time elapsed = 6.140012502670288 s
Completed calculation for time ind = [84 95], time elapsed = 5.650403738021851 s
Completed calculation for time ind = [ 96 107], time elapsed = 5.640556812286377 s
Completed calculation for time ind = [108 119], time elapsed = 5.6404054164886475 s
Completed calculation for time ind = [120 131], time elapsed = 6.055705785751343 s
Completed calculation for time ind = [132 143], time elapsed = 5.766879081726074 s
Completed calculation for time ind = [144 155], time elapsed = 5.683758735656738 s
Completed calculation for time ind = [156 167], time elapsed = 5.638906717300415 s
Completed calculation for time ind = [168 179], time elapsed = 6.075970411300659 s
Completed calculation for time ind = [180 191], time elapsed = 5.683656692504883 s
Completed calculation for time ind = [192 203], time elapsed = 5.779720783233643 s
Completed calculation for time ind = [204 215], time elapsed = 5.735966920852661 s
Completed calculation for time ind = [216 227], time elapsed = 6.16389012336731 s
Completed calculation for time ind = [228 239], time elapsed = 5.760379314422607 s
Completed calculation for time ind = [240 251], time elapsed = 5.806980848312378 s
Completed calculation for time ind = [252 263], time elapsed = 5.7426581382751465 s
Completed calculation for time ind = [264 275], time elapsed = 6.168529033660889 s
Completed calculation for time ind = [276 287], time elapsed = 5.6668009757995605 s

Plot the time-mean $G_{\text{horizontal flux divergence}}$#

plt.figure(figsize=(20,8));
obj = ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, \
                              G_vol_horiz_divergence_depth_integrated_mean,\
                              show_colorbar=True,\
                              cmin=-1e-7, cmax=1e-7, \
                              cmap='RdBu_r', user_lon_0=-67,\
                              dx=map_dx,dy=map_dy);
plt.title('Average $\partial \eta / \partial t$ [m/s] implied by horizontal flux divergence\n Positive = Water out',
          fontsize=20);

../../_images/0b53a30e2c3779d366780e16f0fa23282c30e979a87545da1bd511945f046681.png

Plot one example $G_{\text{horizontal flux divergence}}$#

G_vol_horiz_divergence_depth_integrated.shape

(288, 13, 90, 90)

curr_t_ind = 100
curr_time = ecco_monthly_mean.time.values[curr_t_ind]
curr_timestr = str(curr_time)[0:4]+'-'+str(curr_time)[5:7]

plt.figure(figsize=(20,8));
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, \
                              G_vol_horiz_divergence_depth_integrated.isel(time=curr_t_ind),
                              show_colorbar=True,\
                              cmin=-2e-7, cmax=2e-7, \
                              cmap='RdBu_r', user_lon_0=-67,\
                              dx=map_dx,dy=map_dy);
plt.title('$\partial \eta / \partial t$ [m/s] implied by horizontal flux divergence for '\
          +curr_timestr+'\nPositive = Water out',
          fontsize=20);

../../_images/ddb62c7fba8548cb198f930d5f737b353e4fa30277cb81760486fed557690e66.png

Plot $G_{\text{horizontal flux divergence}}$ at 100m depth#

plt.figure(figsize=(20,8));
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, \
                              G_vol_horiz_divergence.isel(k=10,time=curr_t_ind),
                              show_colorbar=True,\
                              cmin=-2e-6, cmax=2e-6, \
                              cmap='RdBu_r', user_lon_0=-67,\
                              dx=map_dx,dy=map_dy);
plt.title('$\partial \eta / \partial t$ [m/s] at 100 m implied by horizontal flux divergence for '\
          +curr_timestr+'\nPositive = Water out',
          fontsize=20);

../../_images/58768a6788d680a5e98a9b41d2620e075d8c659f9f71e3a515afd12fb60fd019.png

Comparison of LHS and RHS#

Time mean difference of LHR and RHS#

Plot the time-mean difference between the LHS and RHS of the volume budget equation.

#LHS ETA TENDENCY
a = G_total_tendency_mean
# RHS ETA TENDENCY FROM VOL DIVERGENCE AND SURFACE FLUXES
b = G_vol_horiz_divergence_depth_integrated_mean
c = G_surf_fluxes_mean

delta = a - b - c 

plt.figure(figsize=(20,8));
ecco.plot_proj_to_latlon_grid(ecco_grid.XC, ecco_grid.YC, delta, \
                              show_colorbar=True,\
                              cmin=-1e-9, cmax=1e-9, \
                              cmap='RdBu_r', \
                              dx=map_dx,dy=map_dy);
plt.title('Residual $\partial \eta / \partial t$ [m/s]: LHS - RHS ', 
          fontsize=20);

../../_images/4f3ca0c87745d62ce50cdfde64a5ea5efbee0302c52bc8d14214d4a43bb50269.png

The residual of the time-mean tendency terms is essentially zero. We can look at the distribution of residuals to get a little more confidence.

Histogram of residuals#

tmp = np.abs( a-b-c).values.ravel();
plt.figure(figsize=(10,3));

plt.hist(tmp[np.nonzero(tmp > 0)],np.linspace(0, .5e-10,1000));
plt.grid()

../../_images/56a08719b47caa7e5e250d2a51caa38df57bb546addd6b5dea1ab2a99d928114.png

Almost all residuals < $10^{-11}$ m/s. We can close the ETAN budget using UVELMASS, VVELMASS, WVELMASS and oceFWflx.

Note: As stated earlier, we don’t actually need oceFWflx because the surface velocity of WVELMASS is proportional to oceFWflx but we kept it so that the surface forcing term is explicit.

ETAN budget closure through time#

Global average ETAN budget closure#

Another way of demonstrating volume budget closure is to show the global spatially-averaged ETAN tendency terms through time

# LHR and RHS through time
a = G_total_tendency
b = G_vol_horiz_divergence_depth_integrated
c = G_surf_fluxes
# residuals
d = a-b-c

# area_masked = ecco_grid.rA.where(ecco_grid.hFacC.isel(k=0) > 0).compute()
# area_masked_sum = area_masked.sum(dim=('i','j','tile')).values

# take area-weighted mean of these terms
tmp_a=(a*area_masked).sum(dim=('i','j','tile')).compute()/area_masked_sum
# tmp_b=(b*area_masked).sum(dim=('i','j','tile')).compute()/area_masked_sum     # can compute global mean earlier when vol horiz div is computed
tmp_b=G_vol_horiz_divergence_depth_integrated_globalmean
tmp_c=(c*area_masked).sum(dim=('i','j','tile')).compute()/area_masked_sum
# tmp_d=(d*area_masked).sum(dim=('i','j','tile'))/area_masked_sum
tmp_d = tmp_a - tmp_b - tmp_c

# result is four time series
tmp_a.dims

('time',)

fig, axs = plt.subplots(2, 2, figsize=(12,8))

plt.sca(axs[0,0])
tmp_a.plot(color='orange')
axs[0,0].set_title("d(eta)/d(t) total [m/s]")
plt.grid()
plt.ylim([-1e-8, 1e-8]);

plt.sca(axs[0,1])
tmp_b.plot(color='g')
axs[0,1].set_title("d(eta)/d(t) vol divergence [m/s]")
plt.grid()
plt.ylim([-1e-8, 1e-8]);

plt.sca(axs[1,0])
tmp_c.plot(color='r')
axs[1,0].set_title("d(eta)/d(t) surf fluxes [m/s]")
plt.grid()
plt.ylim([-1e-8, 1e-8]);

plt.sca(axs[1,1])
tmp_d.plot(color='b')
axs[1,1].set_title("d(eta)/d(t) LHS - RHS [m/s]")
plt.grid()
plt.subplots_adjust(hspace = .5, wspace=.2)
plt.ylim([-1e-8, 1e-8])
plt.suptitle('Global Volume Budget',fontsize=20);

../../_images/dd27c6eaacd614ddb74ef3384f8535fa387de55f1a8ddcc2acdd4a6f8514b1d5.png

When averaged over the entire ocean surface the volumetric divergence has no net impact on $\partial \eta / \partial t$. This makes sence because horizontal flux divergence can only redistribute volume. Globally, $\eta$ can only change via net surface fluxes.

Local ETAN budget closure#

Locally we expect that volume divergence can impact $\eta$. This is demonstrated for a single point the Arctic.

# close distributed client (removes overhead of distributed client with numerous small chunks)
client.close()

%%time

# Recall, from above...
a = G_total_tendency
b = G_vol_horiz_divergence_depth_integrated
c = G_surf_fluxes
d = a-b-c

plot_tile = 6
plot_j = 40
plot_i = 29
# tmp_aa = a.isel(tile=plot_tile,j=plot_j,i=plot_i)
# tmp_bb = b.isel(tile=plot_tile,j=plot_j,i=plot_i)
# tmp_cc = c.isel(tile=plot_tile,j=plot_j,i=plot_i)
# tmp_dd = d.isel(tile=plot_tile,j=plot_j,i=plot_i)

# use superchunks again (1 year = 1 superchunk) to partition computation
# and improve performance in limited-memory 
tmp_aa = np.empty((a.sizes['time'],))
tmp_bb = np.empty((b.sizes['time'],))
tmp_cc = np.empty((c.sizes['time'],))
size_superchunk = 12
for time_superchunk in range(int(np.ceil(b.sizes['time']/size_superchunk))):
    curr_time_ind = np.arange(size_superchunk*time_superchunk,\
                              np.fmin(size_superchunk*(time_superchunk+1),\
                                      b.sizes['time']))
    tmp_aa[curr_time_ind] = G_total_tendency.isel(time=curr_time_ind,\
                                                  tile=plot_tile,j=plot_j,i=plot_i).compute()
    tmp_bb[curr_time_ind] = G_vol_horiz_divergence_depth_integrated.isel(time=curr_time_ind,\
                                                                         tile=plot_tile,j=plot_j,i=plot_i).compute()
    tmp_cc[curr_time_ind] = G_surf_fluxes.isel(time=curr_time_ind,\
                                               tile=plot_tile,j=plot_j,i=plot_i).compute()
    print(f'Completed time_superchunk = {time_superchunk}')
tmp_dd = tmp_aa-tmp_bb-tmp_cc

tmp_aa = xr.DataArray(data=tmp_aa,dims=['time'],coords=a.isel(tile=plot_tile,j=plot_j,i=plot_i).coords)
tmp_bb = xr.DataArray(data=tmp_bb,dims=['time'],coords=b.isel(tile=plot_tile,j=plot_j,i=plot_i).coords)
tmp_cc = xr.DataArray(data=tmp_cc,dims=['time'],coords=c.isel(tile=plot_tile,j=plot_j,i=plot_i).coords)
tmp_dd = xr.DataArray(data=tmp_dd,dims=['time'],coords=d.isel(tile=plot_tile,j=plot_j,i=plot_i).coords)

fig, axs = plt.subplots(2, 2, figsize=(12,8))
plt.sca(axs[0,0])
tmp_aa.plot(color='orange')
axs[0,0].set_title("d(eta)/d(t) total [m/s]")
plt.ticklabel_format(axis='y', style='sci', useMathText=True)
plt.ylim([-5e-7, 5e-7]);
plt.grid()

plt.sca(axs[0,1])
tmp_bb.plot(color='g')
axs[0,1].set_title("d(eta)/d(t) vol divergence [m/s]")
plt.ticklabel_format(axis='y', style='sci', useMathText=True)
plt.ylim([-5e-7, 5e-7]);
plt.grid()

plt.sca(axs[1,0])
tmp_cc.plot(color='r')
axs[1,0].set_title("d(eta)/d(t) surf fluxes [m/s]")
plt.ticklabel_format(axis='y', style='sci', useMathText=True)
plt.ylim([-5e-7, 5e-7]);
plt.grid()

plt.sca(axs[1,1])
tmp_dd.plot(color='b')
axs[1,1].set_title("d(eta)/d(t) LHS - RHS [m/s]")
plt.ticklabel_format(axis='y', style='sci', useMathText=True)
plt.ylim([-5e-7, 5e-7]);
plt.grid()

plt.subplots_adjust(hspace = .5, wspace=.2)
plt.suptitle('Volume Budget for one Arctic Ocean point',
                 fontsize=20);

Completed time_superchunk = 0
Completed time_superchunk = 1
Completed time_superchunk = 2
Completed time_superchunk = 3
Completed time_superchunk = 4
Completed time_superchunk = 5
Completed time_superchunk = 6
Completed time_superchunk = 7
Completed time_superchunk = 8
Completed time_superchunk = 9
Completed time_superchunk = 10
Completed time_superchunk = 11
Completed time_superchunk = 12
Completed time_superchunk = 13
Completed time_superchunk = 14
Completed time_superchunk = 15
Completed time_superchunk = 16
Completed time_superchunk = 17
Completed time_superchunk = 18
Completed time_superchunk = 19
Completed time_superchunk = 20
Completed time_superchunk = 21
Completed time_superchunk = 22
Completed time_superchunk = 23
CPU times: user 37.9 s, sys: 2.41 s, total: 40.3 s
Wall time: 13.7 s

../../_images/d0ff7e9f43971f62744065a93c788b255ccd914c6cf466e3c4d54fadea3677b3.png

Indeed, the volume divergence term does contribute to $\eta$ variations at this one point.

The seasonal cycle of surface fluxes from sea-ice growth and melt can be seen in the surface fluxes term (plotted below for just the year 1995)

tmp_cc.sel(time='1995').plot(color='r')
plt.ticklabel_format(axis='y', style='sci', useMathText=True)
plt.ylim([-5e-7, 5e-7]);
plt.grid()
plt.title('d(eta)/d(t) from surf. fluxes in 1995 [m/s]');

../../_images/c1e92c6f7dd3f6bbb4527cb6a3daca234516617ac76954605f23fbf53ece2a05.png

Predicted vs. actual $\eta$#

As we have shown, in our Boussinesq model the only term that can change global mean model sea level anomaly $\eta$ is net surface freshwater flux. Let us compare the time-evolution of $\eta$ implied by surface freshwater fluxes and the actual $\eta$ from the model output.

The predicted $\eta$ time series is calculated by time integrating $G_{surface fluxes}$.
This time series is compared against the actual $\eta$ time series anomaly relative to the $\eta(t=0)$.

area_masked = ecco_grid.rA.where(ecco_grid.maskC.isel(k=0) == 1)

dETA_per_month_predicted_from_surf_fluxes = \
    ((G_surf_fluxes * area_masked).sum(dim=('i','j','tile')) / \
     area_masked.sum())*secs_per_month

ETA_predicted_by_surf_fluxes = \
    np.cumsum(dETA_per_month_predicted_from_surf_fluxes.values)

ETA_from_ETAN = \
    (ecco_monthly_snaps.ETAN * area_masked).sum(dim=('i','j','tile')) / \
    area_masked.sum().compute()

# plotting
plt.figure(figsize=(14,5))

plt.plot(dETA_per_month_predicted_from_surf_fluxes.time, \
         ETA_predicted_by_surf_fluxes,'b.')
plt.plot(ETA_from_ETAN.time.values, ETA_from_ETAN-ETA_from_ETAN.values[0],'r-')
plt.grid()
plt.ylabel('global mean $\eta$')
plt.legend(('predicted', 'actual'))
plt.title('$\eta(t)$ as predicted from net surface fluxes and model ETAN [m]', 
          fontsize=20)

../../_images/94036e522d43de25a25fdfce9479d0242076d1fdb87dcae87f692d69692fca72.png

The first predicted $\eta$ occurs at the end of the first month (one month time integral of $\partial \eta / \partial t$). The first actual $\eta$ is set to be zero.

The above plot is another way of confirming that in our Boussinesq model the only term that can change global mean model sea level anomaly $\eta$ is net surface freshwater flux. To account for changes in global mean density we must apply the Greatbatch correction, inverse-barometer correction, and a correction term to account for the fact that sea-ice does not ‘float’ on top of the ocean but in fact displaces seawater upwards. All of these corrections are made for the term SSH, dynamic sea surface height anomaly (not shown here).

One can compare the sea level rise from mass fluxes in ECCO vs those estimated from GRACE and other datasets published by the WCRP Global Sea Level Budget Group: Global sea-level budget 1993-present, Earth Syst. Sci. Data, 10, 1551-1590, https://doi.org/10.5194/essd-10-1551-2018, 2018.

Available here. See Figure 16. https://www.earth-syst-sci-data.net/10/1551/2018/essd-10-1551-2018.pdf

# Annual mean SL calculation must account for different lengths of each month.

# step 1. weight predicted eta by seconds in each month
tmp1=ETA_predicted_by_surf_fluxes*secs_per_month
# step 2, group records by year and sum across each year
tmp2=tmp1.groupby('time.year').sum()
# step 3, group secs per month by year and sum across each year
secs_per_year = secs_per_month.groupby('time.year').sum()

# step 3, divide time-weighted ETA by seconds per year
annual_mean_GMSL_due_to_mass_fluxes =tmp2/secs_per_year
num_years = len(annual_mean_GMSL_due_to_mass_fluxes.year.values)

plt.figure(figsize=(8,5));
# the -0.13 is to make the starting value comparable with WCRP fig 16.
plt.bar(annual_mean_GMSL_due_to_mass_fluxes.year.values[12:num_years],\
        1000*(annual_mean_GMSL_due_to_mass_fluxes.values[12:num_years]-.017), \
        color='k')
plt.grid()
plt.xticks(np.arange(2005, 
                     annual_mean_GMSL_due_to_mass_fluxes.year.values[-1]+1,step=1));
plt.title('Sea level (mm) caused by mass fluxes');

../../_images/858c852a82294f04ed7ef80b0e49dcddfded06cdc7e910454e72de8675d26e3e.png

vs the WCRP data:

from IPython.display import Image
Image('../figures/cazenave_ocean_mass.png')

../../_images/88294898069e70059b9807e77a9e415b98033cda643e64cfff3162b39de4eaf0.png

Time-mean $\partial \eta / \partial t$ due to net freshwater fluxes#

We can easily calculate the time mean rate of global sea level rise due to net freshwater flux.

x=(G_surf_fluxes_mean * area_masked).sum() / area_masked.sum()
total_GMSLR_due_to_mass_fluxes = x*seconds_in_entire_period
print('total global mean sea level rise due to mass fluxes m: ', \
      total_GMSLR_due_to_mass_fluxes.values)

total global mean sea level rise due to mass fluxes m:  0.06982078709392553

Dividing this by the total number of years in the analysis period gives a average rate per year:

total_number_of_years = len(secs_per_month)/12
total_number_of_years

24.0

mean_rate_of_GMSLR_due_to_mass_fluxes = \
    total_GMSLR_due_to_mass_fluxes/total_number_of_years

print('mean rate of GMSLR due to mass fluxes [mm/yr] ',  \
      1000*mean_rate_of_GMSLR_due_to_mass_fluxes.values)

mean rate of GMSLR due to mass fluxes [mm/yr]  2.909199462246897

Compare with other estimates of GMSLR due to mass fluxes.

Remaining components of sea surface height#

# estimate "snapshots" of global mean steric height anomaly at month boundaries
ds_glo_steric = ea.ecco_podaac_to_xrdataset('ster',time_res='daily',\
                                            StartDate='1993-01',EndDate='2017-01',\
                                            mode='s3_open')
glo_steric_monthly_snaps = ds_glo_steric.global_mean_sterodynamic_sea_level_anomaly\
                            .interp(time=ecco_monthly_snaps.time)

ShortName Options for query "ster":
                  Variable Name     Description (units)

Option 1: ECCO_L4_GMSL_TIME_SERIES_DAILY_V4R4    *daily means*
                  global_mean_barystatic_sea_level_anomaly
                                    Global mean of barystatic sea level anomaly (m)
                  global_mean_sterodynamic_sea_level_anomaly
                                    Global mean of sterodynamic sea level anomaly (m)
                  global_mean_sea_level_anomaly
                                    Global mean of dynamic SSH (m)


Using dataset with ShortName: ECCO_L4_GMSL_TIME_SERIES_DAILY_V4R4
{'ShortName': 'ECCO_L4_GMSL_TIME_SERIES_DAILY_V4R4', 'temporal': '1993-01-02,2017-01-31'}

Total number of matching granules: 1

# retrieve sea ice load (mass) per unit area
ds_siceload = ea.ecco_podaac_to_xrdataset('siceload',time_res='snapshot',\
                                          StartDate='1993-01',EndDate='2016-12',\
                                          snapshot_interval='monthly',\
                                          mode='s3_open_fsspec',\
                                          jsons_root_dir=join('/efs_ecco','mzz-jsons'))

ShortName Options for query "siceload":
                  Variable Name     Description (units)

Option 1: ECCO_L4_SEA_ICE_CONC_THICKNESS_LLC0090GRID_SNAPSHOT_V4R4    *native grid,snapshots at daily intervals*
                  SIarea            Sea-ice concentration (fraction between 0 and 1)
                  SIheff            Area-averaged sea-ice thickness (m)
                  SIhsnow           Area-averaged snow thickness (m)
                  sIceLoad          Average sea-ice and snow mass per unit area
                                    (kg/m^2)


Using dataset with ShortName: ECCO_L4_SEA_ICE_CONC_THICKNESS_LLC0090GRID_SNAPSHOT_V4R4

# global mean sea ice load correction
rhoConst = 1029
sIceLoad_m = ds_siceload.sIceLoad/rhoConst
sIceLoad_m_globmean = \
    (sIceLoad_m * area_masked).sum(dim=('i','j','tile')) / \
    area_masked_sum

# adding components to ETAN
ETA_plus_glob_steric = ETA_from_ETAN + glo_steric_monthly_snaps
ETA_plus_glob_steric_sIceLoad = ETA_plus_glob_steric + sIceLoad_m_globmean

# global mean dynamic SSH
SSH_globmean = \
    (ds_dict[ShortNames_list[-2]]['SSH'] * area_masked).sum(dim=('i','j','tile')) / \
    area_masked_sum

# global mean SSH, with inverted barometer effect
SSHNOIBC_globmean = \
    (ds_dict[ShortNames_list[-2]]['SSHNOIBC'] * area_masked).sum(dim=('i','j','tile')) / \
    area_masked_sum



# plotting
fig,ax = plt.subplots(2,1,figsize=(14,10))

ax[0].plot(ETA_from_ETAN.time.values, (ETA_from_ETAN-ETA_from_ETAN[0]).values,'r-')
ax[0].plot(ETA_from_ETAN.time.values, (ETA_plus_glob_steric - ETA_plus_glob_steric[0]).values,'g-')
ax[0].plot(ETA_from_ETAN.time.values, (ETA_plus_glob_steric_sIceLoad \
                                     - ETA_plus_glob_steric_sIceLoad[0]).values,'b-')
# plt.plot(SSH_globmean.time.values, (SSH_globmean - SSH_globmean[0]).values,'k-')
ax[0].grid()
ax[0].set_ylabel('global mean sea level anomaly [m]')
ax[0].legend(('ETAN','ETAN+global steric','ETAN+global steric+sea ice load'))
# plt.legend(('ETAN','ETAN+glob steric','ETAN+glob steric+sea ice load','Dynamic SSH'))
ax[0].set_title('Components of dynamic sea surface height, global mean', 
          fontsize=16)

ax[1].plot(ETA_from_ETAN.time.values, (ETA_from_ETAN-ETA_from_ETAN[0]).values,'r-')
ax[1].plot(ETA_from_ETAN.time.values, (ETA_plus_glob_steric - ETA_plus_glob_steric[0]).values,'g-')
ax[1].plot(ETA_from_ETAN.time.values, (ETA_plus_glob_steric_sIceLoad \
                                     - ETA_plus_glob_steric_sIceLoad[0]).values,'b-')
ax[1].plot(SSH_globmean.time.values, (SSH_globmean - SSH_globmean[0]).values,'k-')
ax[1].plot(SSHNOIBC_globmean.time.values, (SSHNOIBC_globmean - SSHNOIBC_globmean[0]).values,\
           color=(1,.6,0),linestyle=(0,(5,3)))
ax[1].grid()
ax[1].set_ylabel('global mean sea level anomaly [m]')
ax[1].legend(('ETAN','ETAN+glob steric','ETAN+glob steric+sea ice load','Dynamic SSH','SSH w/ atm pressure displacement'))
ax[1].set_title('Components of sea surface height, global mean', 
          fontsize=16)
plt.show()

../../_images/b25bc5e6996a01e4b9ceb7308767169030fe89a6d1d459935cf2ea592f783cfe.png

Dashboard: http://127.0.0.1:8787/status	Workers: 4
Total threads: 4	Total memory: 30.67 GiB
Status: running	Using processes: True

Comm: tcp://127.0.0.1:43177	Workers: 4
Dashboard: http://127.0.0.1:8787/status	Total threads: 4
Started: Just now	Total memory: 30.67 GiB

Comm: tcp://127.0.0.1:40511	Total threads: 1
Dashboard: http://127.0.0.1:34353/status	Memory: 7.67 GiB
Nanny: tcp://127.0.0.1:46295
Local directory: /tmp/dask-scratch-space/worker-cnp08wqq

Comm: tcp://127.0.0.1:39179	Total threads: 1
Dashboard: http://127.0.0.1:35527/status	Memory: 7.67 GiB
Nanny: tcp://127.0.0.1:35897
Local directory: /tmp/dask-scratch-space/worker-8wt5d6fv

Comm: tcp://127.0.0.1:35659	Total threads: 1
Dashboard: http://127.0.0.1:41823/status	Memory: 7.67 GiB
Nanny: tcp://127.0.0.1:45545
Local directory: /tmp/dask-scratch-space/worker-ln82bczw

Connection method: Cluster object	Cluster type: distributed.LocalCluster
Dashboard: http://127.0.0.1:8787/status

Comm: tcp://127.0.0.1:41919	Total threads: 1
Dashboard: http://127.0.0.1:36037/status	Memory: 7.67 GiB
Nanny: tcp://127.0.0.1:39659
Local directory: /tmp/dask-scratch-space/worker-j449fv1k

ECCOv4 Global Volume Budget Closure

Contents

ECCOv4 Global Volume Budget Closure#

Objectives#

Introduction#

ETAN in a Boussinesq Model#

Evaluating the model sea level anomaly ETAN volume budget#

Datasets to download#

Prepare environment and loading the relevant model variables#

Client

Cluster Info

LocalCluster

Scheduler Info

Scheduler

Workers

Worker: 0

Worker: 1

Worker: 2

Worker: 3

Load ecco_grid#

Open 2D MONTHLY \(\eta\) snapshots#

Open MONTHLY mean data#

Create the xgcm ‘grid’ object#

Calculate LHS: \(\eta\) time tendency: \(G_{\text{total tendency}}\)#

Plot the time-mean \(\partial \eta / \partial t\), total \(\Delta \eta\), and one example \(\partial \eta / \partial t\) field#

Time-mean \(\partial \eta / \partial t\)#

Total \(\Delta \eta\)#

Example \(\partial \eta / \partial t\) field#

Calculate RHS: \(\eta\) tendency due to surface fluxes, \(G_\text{surface fluxes}\)#

Plot the time-mean, total, and one month average of \(G_{\text{surface fluxes}}\)#

Time-mean \(G_{\text{surface fluxes}}\)#

Total \(\Delta \eta\) due to surface fluxes#

Example \(\partial \eta / \partial t\) impliedy by surface fluxes#

Calculate RHS: \(\eta\) tendency due to volumetric flux divergence, \(G_\text{volumetric fluxes}\)#

Vertical volumetric flux divergence#

Plot the time-mean \(G_{\text{vertical flux divergence}}\)#

Horizontal Volume Flux Divergence#

Plot the time-mean \(G_{\text{horizontal flux divergence}}\)#

Plot one example \(G_{\text{horizontal flux divergence}}\)#

Plot \(G_{\text{horizontal flux divergence}}\) at 100m depth#