BISMG:VickyL

Flow beneath ice shelves
Aim is to model the flow of ocean water into the cavity beneath an ice shelf and its interaction with the overlying ice. We want to capture the effect of intrusions of warm water from the surrounding oceans on the basal melt rate of the ice shelf.

Conceptual model
Melt water from the base of the ice shelf and from the foot of an ice stream is assumed to form a layer adjacent the underside of the shelf. The layer flows along the underside to the shelf's front. As it flows ocean water in the cavity is entrained into the melt-water layer.

Model Assumptions

 * 1) Pressure is assumed hydrostatic.
 * 2) Pressure is continuous across layer interfaces

Mathematical model
$$\nabla. \Big( \frac{1}{D+H} \nabla \psi \Big) = \zeta$$

$$\zeta = - \frac{f}{c_{d}} \nabla. \overline{\boldsymbol{u}} - \frac{1}{\rho_0 c_{d}} \boldsymbol{k}. \nabla \times \overline{\nabla p_{b}}$$

$$\overline{u} = - \frac{1}{(D+H)}\frac{\partial \psi}{\partial y} \qquad \overline{v} = \frac{1}{(D+H)}\frac{\partial \psi}{\partial x}$$

$$f \boldsymbol{k} \times \bar{\boldsymbol{u}_{p}} + c_{d} \bar{\boldsymbol{u}_{p}} =  f \boldsymbol{k} \times \overline{\boldsymbol{u}}  +  c_{d} \overline{\boldsymbol{u}} +   \frac{g}{ \rho_{0} (D+H)} \left[ \frac{DH}{2} \nabla \rho_{p} - \right. \left. H (\rho_{c} - \rho_{p}) \nabla I + \frac{H^{2}}{2} \nabla \rho_{c} \right] $$

$$\bar{\boldsymbol{u}}_{c} = \frac{D+H}{H} \overline{\boldsymbol{u}} - \frac{D}{H} \bar{\boldsymbol{u}}_{p}$$

$$\frac{\partial D}{\partial t} + \nabla. \left( D \bar{\boldsymbol{u}_{p}} \right) = \dot{e}$$

$$H = h - B - D$$

$$\frac{\partial T_{p}}{\partial t} + \bar{\boldsymbol{u}_{p}} \nabla. T_{p} = K_{h} \nabla^{2} T_{p} + \frac{\gamma_{T}|\bar{\boldsymbol{u}_{p}}|}{D} (T_{B}-T_{p}) + \frac{\dot{e}}{D} (T_{I}-T_{p})$$

$$\frac{\partial S_{p}}{\partial t} + \bar{\boldsymbol{u}_{p}} \nabla. S_{p} = K_{h} \nabla^{2} S_{p} + \frac{\gamma_{S}|\bar{\boldsymbol{u}_{p}}|}{D} (S_{B}-S_{p}) + \frac{\dot{e}}{D} (S_{I}-S_{p})$$

$$\frac{\partial T_{c}}{\partial t} + \bar{\boldsymbol{u}_{c}} \nabla. T_{c} = K_{h} \nabla^{2} T_{c} + \delta_{D} \frac{\gamma_{T}|\bar{\boldsymbol{u}_{c}}|}{H} (T_{B}-T_{c}) - \frac{\dot{e}}{H} (T_{I}-T_{c})$$

$$\frac{\partial S_{c}}{\partial t} + \bar{\boldsymbol{u}_{c}} \nabla. S_{c} = K_{h} \nabla^{2} S_{c} + \delta_{D} \frac{\gamma_{S}|\bar{\boldsymbol{u}_{c}}|}{H} (S_{B}-S_{c}) - \frac{\dot{e}}{H} (S_{I}-S_{c})$$

Indexing
PetSc routines are written in C where vectors and arrays have a lower index of 0. Fortran 90 vectors and arrays start from 1 by default. This how to put a Petsc Vec variable into a 1d allocatable array

real(dp), dimension, allocatable :: f90arr Vec :: petscvector PetscInt :: N PetscInt, dimension, allocatable :: indx

! Create the PetSc vector call VecCreateSeq(PETSC_COMM_SELF,N,petscvector,ierr)

... fill petscvector

allocate(indx(N), f90arr(N))

! Create array of indices for petscvector to extract the data do i=1,N k=i-1 indx(i)=k end do

! Put data in petscvector into the f90 array f90arr call VecGetValues(petscvector,N,indx,f90arr,ierr)

f90arr(1) contains the first element of petscvector f90arr(0) is undefined

Work in progress
Suggestions and comments from JCRP meeting (13/11/09):
 * Why is diffusion of scalars (T,S) important near the inflow along the grounding line? Surely the scalar difference is driving the plume?
 * Use the simplified momentum equation derived under the assumption of stationary cavity for the motion of the plume with the moving cavity. We know this equation works. Continuity of the whole cavity is satisfied through the streamfunction. Velocity in the cavity is simply the difference between the depth weighted velocity in the whole cavity calculated from the streamfunction and the plume velocity. Question: can we really justify using a stationary ocean assumption just to calculate the plume dynamics and then ignore it?
 * Why not ignore the pressure gradient due to the base of the ice shelf and force term due to the slope of the ice shelf in plume momentum? These terms exist even if there is no flow beneath the ice shelf - they do not drive motion so they should be ignored. In the above model the pressure gradient and the force due to the slope of the ice shelf is found from the barotropic momentum equation and then substituted into the plume equation. Crucially, this substitution introduces a coupling to the cavity flow through the reduced gravity term. By ignoring the two 'forces' the cavity only influences the plume through entrainment in the thickness and scalar equations. There is a reduced gravity term in the momentum equation for the cavity. Perhaps we could calculate the cavity flow directly and then the plume velocity would be the difference between the depth weighted barotropic and cavity velocities.
 * What would happen if the cavity flow went to zero? I'm not sure how to switch off the cavity flow in the above model. It was suggested that the drag coefficient for the cavity be increased to slow the flow down. The cavity momentum is not calculated directly in this model. I would need to increase the drag in the barotopic momentum equation.

The week ahead

 * Modify derivatives of scalars onto scalar grid points, i.e. gradients used in Jacobians, to take account of known boundary conditions. T, S and $$\rho$$ are insulated along the land/sea boundary so the normal gradient is zero.
 * Need to distinguish between land/sea boundary and open boundary when differentiating. Do this through masks. How about inflow boundary on the grounding line?
 * Try solving the plume momentum equations in the form of linear equations by introducing $$\frac{ \partial \boldsymbol{u_{p}}}{ \partial t}$$ and marching the solution forward rather than treating the momentum equations as two simultaneous equations. This has the advantage of not requiring cross derivatives.
 * Investigate the inflow boundary on the grounding line.
 * Play with solvers in PetSc subroutine.

Friday 27th November 2009

 * Played with a PetSc example of a general elliptical problem, the 2d diffusion equation. Wrote a code to solve the example on p261 of Numerical Solution of PDEs, Morton & Mayers. Main difficulty was extracting the solution from a PetSc variable to a fortran array. Stephen helped.
 * Included a petsc subroutine into the plume model to solve the streamfunction-vorticity equation. Initial runs with PIG and simple geometry (bathymetry_type = 12) set ups suggest petsc works with Dirichlet and Neumann boundary conditions.

Friday 20th November 2009

 * Slapsolv cannot solve the streamfunction-vorticity equation with a Dirichlet condition on the land/sea boundary AND a Neumann condition on the open boundary. Peric's linear solves can handle the problem but are painfully slow. Need to try something else.
 * Gethin installed PetSc on Dartagnan.
 * Read through PetSc manual and selected a suitable example.

Friday 13th November 2009

 * Adjusted the cross derivatives of scalar quantities onto velocity grid points. Previously, I used second order one-sided finite difference using 3 scalar points adjacent to the boundary or first order one-sided difference if only 2 scalar point were available. Now the derivatives use known boundary values if available. For the streamfunction the value on the land boundary is known so I use second order difference formulation for a nonuniform grid (the spacing between the boundary value and the adjacent cell centre is $$\Delta x/2$$. For T,S and $$\rho$$ the normal gradient to the land is zero (i.e. insulating boundary). I calculate the gradient half a grid cell above the required velocity point using central difference and then average it with the boundary value (i.e. divide by 2 since the boundary value is zero). The results were not that different from the previous version. It did reduce the variations the divergence of barotropic velocity used in the streamfunction-vorticity equation around the grounding line in the PIG simulations.
 * Introduced two-way entrainment.

$$ \dot{e} = \begin{cases} \frac{c_{L}^{2}}{S_{m}} \left( |\bar{\boldsymbol{u}_{p}}-\bar{\boldsymbol{u}_{c} }|^{2} + \frac{g ( \rho_{c}-\rho_{p} ) D}{\rho_{0} S_{m}}\right)^{\frac{1}{2}} & if \,\, D > D_{\min} \,\, and \,\, |\bar{\boldsymbol{u}_{p}}| > |\bar{\boldsymbol {u}_{c}}| \\ -\frac{c_{L}^{2}}{S_{m}} \left( |\bar{\boldsymbol{u}_{p}}-\bar{\boldsymbol{u}_{c} }|^{2} + \frac{g ( \rho_{c}-\rho_{p} ) H}{\rho_{0} S_{m}}\right)^{\frac{1}{2}} & if \,\, D > D_{\min} \,\, and \,\, |\bar{\boldsymbol{u}_{p}}| \le |\bar{\boldsymbol {u}_{c}}| \\ 0 & otherwise \end{cases} $$

where the Schmidt number is a function of the Richardson number which is now defined as

$$ Ri = \begin{cases} \max \left( 0.05, \frac{g ( \rho_{c}-\rho_{p} ) D}{\rho_{0} |\bar{\boldsymbol{u}_{p}}-\boldsymbol{u}_{c}|^{2}} \right) & if \,\, |\bar{\boldsymbol{u}_{p}}| > |\bar{\boldsymbol {u}_{c}}| \\ \max \left( 0.05, \frac{g ( \rho_{c}-\rho_{p} ) H}{\rho_{0} |\bar{\boldsymbol{u}_{p}}-\boldsymbol{u}_{c}|^{2}} \right) & if \,\, |\bar{\boldsymbol{u}_{p}}| \le |\bar{\boldsymbol{u}_{c}}| \end{cases} $$
 * Values of T and S at the interface are set to the depth weighted average of the plume and cavity values such as $$T_{I}= \frac{D T_{p} + H T_{c}}{D+H}$$. Previously use the mean value between the plume and cavity values. Result for PIG is that temperature in the cavity near the inflow on the grounding line does not drift above its maximum initial value (i.e. 1 degree C).