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

 * 1) $$\nabla . \Big( \frac{1}{D+H} \nabla \psi \Big) = \zeta$$
 * 2) $$\zeta = - \frac{f}{c_{d}} \nabla . \overline{\boldsymbol{u}} - \frac{1}{\rho_0 c_{d}} \boldsymbol{k} . \nabla \times \overline{\nabla p_{b}}$$
 * 3) $$\overline{u} = - \frac{1}{(D+H)}\frac{\partial \psi}{\partial y} \qquad \overline{v} =  \frac{1}{(D+H)}\frac{\partial \psi}{\partial x}$$
 * 4) $$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] $$
 * 5) $$\bar{\boldsymbol{u}}_{c} = \frac{D+H}{H} \overline{\boldsymbol{u}} - \frac{D}{H} \bar{\boldsymbol{u}}_{p}$$
 * 6) $$\frac{\partial D}{\partial t} + \nabla . \left( D \bar{\boldsymbol{u}_{p}} \right) = \dot{e}$$
 * 7) $$H = h - B - D$$
 * 8) $$\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})$$
 * 9) $$\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})$$
 * 10) $$\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})$$
 * 11) $$\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

 * Need to distinguish between land/sea boundary and open boundary when differentiating - no need for plume velocity as one-sided differencing not used at open boundary only on land/sea boundary. Do this through masks. How about inflow boundary on the grounding line?
 * Investigate the inflow boundary on the grounding line.
 * Play with solvers in PetSc subroutine.
 * Look at open boundary conditions for the barotropic flow.
 * Do we want to use PetSc as a thickness and scalar solver?

Friday 26th February 2010

 * Looked how the streamfunction-vorticity equation is solved. The linear system of equations of the elliptical equation can be expression as $$A \psi =b$$. It turns out that $$A$$ need only be set up once because its elements contain values of $$H(x,y)$$ which do not vary during the simulation. If the inverse of $$A$$ could be calculated on the first timestep then the streamfunction could subsequently be found by a matrix multiplication, $$\psi = A^{-1} b$$.
 * PetSc can calculate inverse matrices directly, but it's not recommended. Steph suggested using the full LU decomposition as a preconditioner for A instead. This means that the inverse is calculated on the first call then the system is solved with one iteration on subsequent calls.
 * The fastest combination so far is pc_type lu and ksp_type preonly which took 4m5s to complete 1 day of simulation with dt=60s. It the default pc and ksp types were used the it took 7m55s.
 * Looked at using PetSc for the thickness solver again. Something is going wrong.

Friday 8th January 2010

 * Implemented three level timestep incorrectly. It removed the oscillation in the simple geometry testcase because the time scheme behaved like the first order implicit scheme with a smaller timestep. With the correct implementation oscillation is accentuated (as noted by Ferziger and Peric) and in fact the PIG set up crashes.
 * Implemented 2nd order Crank-Nichols time stepping scheme for the thickness solver.
 * Reducing the time step (dt=30s) for the simple geometry removes the oscillation along the west wall.
 * First attempt at implementing PetSc to solve the thickness equation. It runs but the plume depth grows too large at times but is advected away for PIG. Also, there is spurious "inflow" at the west boundary along the grounding line.

Friday 18th December 2009

 * Ran PIG simulations upside down to look at the way the north and east boundaries handle the open boundaries. The plume thickens along the southern/northern boundary partly due to the flow hugging the coastline/grounding line not passing out of the domain properly, particularly with the no normal gradient condition, and partly due to a branch of the plume flowing close to the boundary. This thickening is more prominent along the northern boundary. Waves seem to propagate from the branch passing close to the boundary along the northern boundary towards the outflow adjacent to the coastline. These waves may dampen or grow in time depending on the choice of open boundary condition.

Open boundary conditions \begin{cases} \theta_{B}^{k} & if \, \, c \le 0 \\ \frac{\Delta n \theta_{B}^{k} + c \Delta t \theta_{B-1}^{k+1}}{ \Delta n + c \Delta t} & if \, \, 0 < c \le \frac{\Delta n}{\Delta t} \\ \frac{1}{2} (\theta_{B}^{k}+\theta_{B-1}^{k+1}) & if \, \, c > \frac{\Delta n}{\Delta t} \end{cases}$$. For the simple geometry test case the open boundary behaves like a no normal gradient condition. For the PIG set up with explicit thickness solver the behaviour was similar to the no normal gradient condition except the waves along the northern boundary took longer to dampen. The implicit thickness solver extrapolates the plume depth to the open boundary when setting up the matrix for the problem. This meant that the behaviour along the open boundary was relatively good with flow leaving the domain along the coastline properly and the initial wave along the boundary from the branch of the plume moving close to the boundary was damped quickly. Flow was similar to the backward no normal gradient condition (see below). Phase velocities along the northern boundary for PIG were typically about 10, note that $$\frac{\Delta n}{\Delta t} = 16.66$$, with a few negative values particularly near the branch. \begin{cases} \theta_{ocean} & if \, \, c \le 0 \\ \frac{\Delta n \theta_{B}^{k} + c \Delta t \theta_{B-1}^{k+1}}{ \Delta n + c \Delta t} & if \, \, 0 < c \le \frac{\Delta n}{\Delta t} \\ \frac{1}{2} (\theta_{B}^{k}+\theta_{B-1}^{k+1}) & if \, \, c > \frac{\Delta n}{\Delta t} \end{cases}$$ where $$D_{ocean}=0$$, $$\bar{\boldsymbol{u}}_{p \, ocean}=0$$, $$T_{ocean} = T_{amb}(-B)$$ and $$S_{ocean} = S_{amb}(-B)$$. Test cases ran but the results show that outflow along the coastline pinched off towards the land and the northern boundary was noisy as the sign of the phase velocity alternated.
 * No normal gradients. This is used in the original code probably because it's relatively robust for complex geometries. It's limitations are obvious in the simple geometry test case (bathymetry type 12) were the flux at the open boundary is too small and there is some flow along the northern boundary rather passing out of it. For the (upside down) PIG geometry with the explicit thickness solver this condition works reasonably well. The waves along the northern boundary damp moderately quickly.
 * Extrapolating u, v, D, T and S to the open boundary. Works very well for simple geometry test case where the flow is well behaved but crashes PIG. Problems arise close to the coastline outflow where velocity fluctuate wildly and at the boundary where the branch of plume passes close by because extrapolation accentuates errors, e.g. extrapolating depth can lead to negative values. Negatives depths can be removed by setting the grid value to zero, large values are harder to restricit.
 * Radiative boundary conditions using the Sommerfeld radiation condition $$\frac{\partial \theta}{\partial t} + c \frac{\partial \theta}{\partial n} = 0$$ where $$\theta$$ is any variable, $$n$$ then normal direction and $$c$$ is the phase velocity. Based on Numerical Modeling of Ocean Dynamics, Kowalik and Murty, p311, the phase velocity is calculated using $$c= - \frac{\partial \theta / \partial t}{\partial \theta / \partial n} = - \frac{(3 \theta_{B-1}^{k+1} - 4 \theta_{B-1}^{k} + \theta_{B-1}^{k-1})/2 \Delta t}{( \frac{1}{2} (\theta_{B-1}^{k}+\theta_{B-1}^{k-1})-\theta_{B-1}^{k+1})/ \Delta n} $$. The boundary value is set as follows $$\theta_{B}^{k+1} =
 * Radiative condition with open ocean value for negative phase velocity such that $$\theta_{B}^{k+1} =
 * Backward no normal gradient uses $$\theta_{B}^{k+1}=\theta_{B-1}^{k}$$, where $$k$$ denotes the time step and $$B$$ the boundary node. This is based on the assumption that the flow can be described by a 1st order wave equation. The wave equation has a solution of the form $$\theta(n-ct)$$, i.e. the distribution of $$\theta$$ does not change in time only it's position. If we assume our variables are being advected out of the domain normal to the open boundary, then the boundary value is approximately equal to the value of the variable on the adjacent node at the previous time step. Works well for the simple geometry test case. For PIG, used with the explicit thickness solver the result is similar to the no normal gradient condition, but waves along the northern boundary aren't damped completely. The result with the implicit solver (using extrapolation condition to set the matrix) is good. There is a large wave along the northern boundary initially but is dampened.
 * Recommendations:
 * 1) None of the above conditions worked perfectly, but they did not affect the flow in the main part of the domain. I suggest that applying an open boundary at the ice shelf front will be a mistake. It will be better to keep the open boundary away from complex flow patterns. Make the computational domain larger than is required by extending the ice shelf artificially but only analyse results below the 'real' ice shelf as for the PIG set up.
 * 2) Extrapolating the variables worked best for simple flows so try first.
 * 3) Use extrapolation to set up the matrix for the implicit thickness solver even if another type of condition is used to set the boundary values.
 * 4) If the flow is too complex for extrapolation, use the backward no normal gradient condition with the implicit thickness solver.
 * 5) If using the explicit thickness solver stick with the no normal gradient condition.
 * 6) The radiative condition may be "well-posed" but it's a lot of work for little return.
 * Implemented 2nd order time derivative into the implicit thickness solver routine. Removed the spurious oscillation along the western boundary in the simple geometry test case. Smoothed, but did not remove, the boundary oscillation for PIG set up.

Friday 11th December 2009
Adding only the time derivative to the reduce physics model does not help stability. Instead it introduces spurious waves. As Adrian Jenkins suggested, it's the nonlinear advection terms that prevent oscillations.
 * Tried 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. In the end I used the full nonlinear subroutine and commented out the relevant terms. For both the simple geometry (with extrapolated bc) and PIG setup (no gradient bc) I found the following
 * 1) Removing viscosity and advection terms: waves propagate from the corners of the inflow and corners of the lateral boundaries. Looks exactly like ripples in water. For the simple geometry the simulation crashes when 'waves' hit the open boundary.
 * 2) Removing advection, but retaining viscosity: waves propagate from the corners as above but are damped. At long time damped waves travel along the southern wall but not along the west wall. Flow leaves the open boundary without a problem.
 * 3) Removing viscosity, but retaining viscosity: no noticeable difference with fully nonlinear momentum.
 * For the reduced physics model I tried removing spurious oscillation along the western land boundary for the simple geometry case by:
 * 1) upwinding rather than central differencing. Didn't stop oscillation. Delayed onset as the velocity and depth were damped slightly.
 * 2) introducing radiative open boundary conditions (Kowalik and Murty, Numerical Modeling of Ocean Dynamics p311). A mixture of the field variable at different times. Didn't stop the oscillation.
 * 3) try solving the momentum equations as above and introduce nonlinear terms. Did stop oscillation
 * Extrapolating open boundary didn't work for PIG. Got negative depths and spurious velocities. Wave developed across the boundary when I tried to correct for these.
 * North and east boundaries handled slightly differently from south and west ones. This is due to the (incorrect) handling of the staggered C-grid. In the original code the variables in cells adjacent to the edge of the computational domain are not calculated directly, they are set using the no normal gradient boundary condition instead. However, there is no difference in the variable arrangement for v(:,1) as v(:,n-1) and u(1,:) as u(m-1,:) but in the code the first u and v are set by the boundary condition where as the the second are calculated.

v-  |        | |  T    u   |        | --

Friday 4th December 2009
$$f \boldsymbol{k} \times \bar{\boldsymbol{u}_{p}} + c_{d} \bar{\boldsymbol{u}_{p}} =  f \boldsymbol{k} \times \bar{\boldsymbol{u}_{c}}  +  c_{d} \bar{\boldsymbol{u}_{c}} +   \frac{g}{ \rho_{0}} \left[ \frac{D}{2} \nabla \rho_{p} - \right. \left. (\rho_{c} - \rho_{p}) \nabla I + \frac{H}{2} \nabla \rho_{c} \right] $$
 * Modified 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. The normal derivative of the scalar adjacent to wall is calculated by first calculating the gradient across the cell face opposite the wall using central difference and then taking the mean value (or linear interpolating) with the known value on the boundary.
 * The momentum equation for plume flow above a stationary, homogeneous cavity is recovered by substituting an expression for barotropic velocity in terms of plume and cavity velocities (eqn 5) into the momentum equation (4) above

and setting $$\bar{\boldsymbol{u}_{c}}=0$$ and $$\nabla \rho_{c}=0$$
 * Tested the plume velocity by switching off the cavity coupling as above on the simple geometry set up (bathymetry type 12). Switching off the couple made very little difference to the flow pattern - the depth of the plume was still shallow and did not travel along the western boundary and the velocity field was weak compared to the reduced physics model. The difference between the current and reduced physics model is the drag term. The current model drag uses a linear law $$c_{d} \bar{\boldsymbol{u}_{p}}$$ but the reduced physics model uses a quadratic law and is a function of depth $$ \frac{c_{d} |\bar{\boldsymbol{u}_{p}}|}{D} \bar{\boldsymbol{u}_{p}}$$. For the reduced model, plume speed is about ~0.1 m/s or less and depth is ~100 m along the western boundary and ~10 m elsewhere, so the magnitude of the drag term here is much smaller than that of the current model. This accounts for the difference between the two flow patterns.
 * Reduced the drag coefficient for the current model with the simple geometry set up. Reducing it by a factor 10 to 0.0003 was okay, though a little messy around the inflow. Much smaller and the code crashed. Depth and velocity of the plume did not greatly increase.
 * Went back to the reduced physics model and looked at solving the velocity fields separately rather than using the solution to the simultaneous equations. One can arrange the equations such that either the velocity in the drag term is the subject or in the Coriolis force. Using the drag term means that only the along stream derivatives are required. However, since the drag is quadratic we would have to divide by the linearised speed which rules this approach out (initial velocity is zero). Using the Coriolis force only the cross derivatives are used. It was particularly problematic near the left hand corner of the inflow for the simplified geometry set up (bathymetry type 12). Numerical oscillation appeared to radiate away from the wall. When the solution to the simultaneous equation is used this is not a problem. The use of derivatives in two directions seems to (spatially) stabilise the flow. It should be noted that even with this method the velocity at the inflow alternates values with each timestep (i.e. temporal effect) for the initial 'spin up' and then settles down.
 * Noticed problems with the reduced physics model for the simplified geometry set up (bathymetry type 12) at long time. The open boundary conditions are not very good; the flow is not leaving the computational domain probably. Instead the boundary is behaving like a porous wall. Also, a spurious oscillation appears to start on the top left hand corner where a wall meets the open boundary and travels along the left hand wall. The plume is relatively fast moving and deep here. This oscillation is large but does not destabilise the flow. It is not present in the fully nonlinear model.
 * Replaced the no normal gradient open boundary condition for the northern boundary of the simplified geometry set up (bathymetry type 12) with extrapolation of T, S, $$\boldsymbol{u}_{p}$$ and D. Much improved. Flow is leaving domain properly. Need to alter the implicit solver to include extrapolation.
 * Spurious oscillation still a problem.

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).