Preliminary Model¶

The preliminary model solves a system of coupled, linear equations. It was initially set up on a 1D domain with a zero background velocity field imposed, and was compared with an analytical solution for consistency.

Next, on a 2D domain. Three different divergenceless background velocity fields were imposed: the zero field, a single cell overturning field, and a two-cell overturning field with upwelling at the domain edges and convergent downwelling in the centre. The single-cell and two-cell velocity field calculations are initialized on the zero-field steady state. Two different advection schemes were implemented: an upstream scheme, and a flux based scheme. The flux based scheme has so far demonstrated more robust stability.

Points of interest are the time elapsed to reach steady state, the grid resolution required for 90% similarity, and of course the distribution of [Th]/[Pa].

The objective of the model is to provide insight into the overturning rate of the Arctic ocean by exploring how the ratio of [Th]/[Pa] responds to forcing from overturning currents.

Coupled 1D Model¶

The steady state is reached after 50 years.
The grid is 100 x 1, but 20 would be sufficient.
The solution propogates from surface to seafloor, and linearly increases in magnitude with depth.

The curve corresponding to 10 years shows how the steady state solution, a linear positive slope, is reached first at the surface, where the initial distribution falls out and is not replaced by anything falling from above. The initial distribution has fallen to the 2000m depth level after 10 years, 4000m after 30 years, and completely after 50 years.

The time for a particle to sink from the surface to the seafloor is given by t = zmax / S, however, the time to reach steady state is greater than this, because sinking particles can desorb into the dissolved phase, which do not sink.

Coupled 2D Model: Underlying Physics¶

The steady state for Th is reached after 100 years, while the steady state for Pa is reached after 1000 years. The order of magnitude difference between the two is caused by the order of magnitude difference in the production, adsorption, and desorption constants that govern the two elements.
Initializing the single-cell and two-cell calculations on the zero-field steady state does not change the time to reach steady state.
Both solutions propogate from surface to seafloor.
Upwelling currents are reflected in the concentration profile as regions of high concentration. This is caused by the upward force that suspends particulate matter in the water column.
The domain is 5000m by 1000 km, and is discretized into a 20 x 30 grid, which is quite course, but sufficient to resolve 90% of the solution.

Results from upstream calculation:

Results from flux based calculation:

Results from upstream calculation:

Note the unstability that develops after 600 years in the upstream calculation:

Results from flux based calculation:

Results from upstream calculation:

Coupled 2D Model: Steady State¶

[Th]

zero velocity: 100 years
single cell velocity: 100 years
two cell velocity: 100 years

[Pa]

zero velocity: 1000 years
single cell velocity: 1000 years
two cell velocity: 1000 years

[Th] / [Pa]

zero velocity:
single cell velocity:
two cell velocity:

Coupled 2D Model: Th2D.py¶

Th2D.adflow(T, V, u, nz, nx, k_ad, k_de, Q, flowfig)¶

Compute and store the dissolved and particulate [Th] profiles, write them to a file, plot the results.

Parameters:

T (int) – scale for tmax such that tmax = T*(g.zmax - g.zmin)/S
V (int) – scale for ux, uz, which are originally order 1.
u (float) – 3D tensor of shape [nz, nx, 2]. Stores z component of velocity in [:, :, 1], x component of velocity in [:, :, 2]
nz (int) – number of grid points in z dimension
nx (int) – number of grid points in x dimension
k_ad (float) – nz x nx adsorption rate matrix
k_de (float) – nz x nx adsorption rate matrix
adscheme (function) – function to implement the desired advection scheme

Th2D.u_simple(xmin, xmax, zmin, zmax, nx, nz)¶

Compute a simple rotational, divergenceless flow field on a specified grid.

Parameters:	xmin – minimum x on the grid xmax – maximum x on the grid zmin – minimum z on the grid zmax – maximum z on the grid nx – number of points in x dimension nz – number of points in z dimension

Th2D.u_complex(xmin, xmax, zmin, zmax, nx, nz)¶

Compute a rotational, downwelling velocity field.

Parameters:	xmin – minimum x on the grid xmax – maximum x on the grid zmin – minimum z on the grid zmax – maximum z on the grid nx – number of points in x dimension nz – number of points in z dimension

Th2D.k_sorp(string, xmin, xmax, zmin, zmax, nx, nz)¶

Compute adsorption,desorption, & production constants for Th or Pa.

Parameters:	string – a string, either ‘Th’ or ‘Pa’ xmin – minimum x on the grid xmax – maximum x on the grid zmin – minimum z on the grid zmax – maximum z on the grid nx – number of points in x dimension nz – number of points in z dimension

Th2D.plotratio(DTh, DPa, PTh, PPa, xmin, xmax, zmin, zmax, nx, nz, T)¶

Plot the ratio T/P and output to notebook.

Parameters:

DTh – 2D profile of dissolved Th
PTh – 2D profile of particulate Th
DPa – 2D profile of dissolved Pa
PPa – 2D profile of particulate Pa
xmin – minimum x on the grid
xmax – maximum x on the grid
zmin – minimum z on the grid
zmax – maximum z on the grid
nx – number of points in x dimension
nz – number of points in z dimension
T (int) – scale for tmax such that tmax = T*(g.zmax - g.zmin)/S