COMSOL activities

From Interactive System for Ice sheet Simulation
Revision as of 10:09, 6 August 2009 by Sprice (Talk | contribs)

Jump to: navigation, search

Contents

Overview

Now, let's see if COMSOL can be used to solve problems of glaciological relevance. We'll look at shallow ice approximation flow, and shallow shelf approximation flows.

Isothermal Shallow Ice Approximation

Begin with the often used shallow ice form for ice thickness evolution, which casts evolution as a non-linear diffusion problem

\frac{\partial H}{\partial t} = - \nabla D \nabla H + M

where

D = \frac{2A(\rho g)^n}{n+2} H^{n+2} \left[\nabla H \cdot \nabla H \right]^{(n-1)/2}

with boundary condition H=0 on the edge of the computational domain.

Comsol Modeling

We will use the PDE, General Form transient mode to solve this equation. For convenience, make the dependent variable H.

Geometry

You should not find it difficult to create a unit square. Once it's made, you can double click it to change it's size and do other transformations. Read below to find the appropriate domain.

Field equations

This equation mode solves equations of the form

e_a\frac{\partial^2 H}{\partial t^2} + d_a \frac{\partial H}{\partial t} + \nabla \cdot \gamma = F

Which is just what we want if we recognize that in our system e_a=0, d_a=1, F=M, and

\Gamma_x = -\frac{2A(\rho g)^n}{n+2} H^{n+2} \left[\nabla H \cdot \nabla H \right]^{(n-1)/2} \frac{\partial H}{\partial x}
\Gamma_y = -\frac{2A(\rho g)^n}{n+2} H^{n+2} \left[\nabla H \cdot \nabla H \right]^{(n-1)/2} \frac{\partial H}{\partial y}

Now the problem has been reduced to one of typing. It will make the COMSOL model easier to read if you create a scalar expression for D . Then your \Gamma_x = -D \frac{\partial H}{\partial x} and \Gamma_x =- D\frac{\partial H}{\partial y} are very clear.

Boundary conditions

This type of problem requires a Dirchlet boundary condition. Set H = 0 on all four sides.

Other

You'll also be needing to know how to tell COMSOL to use a derivative. That is Hx, Hy, and Ht for \frac{\partial H}{\partial x}, \frac{\partial H}{\partial y}, and \frac{\partial H}{\partial t} respectively.

Exercises

  1. Complete the model, and do the isothermal fixed margin experiment Huybrechts (1996)[1]. You'll find all values of constants there as well. Verify that your model is providing results consistent with those reported in the paper.
  2. Now alter your model (the accumulation field) to do the isothermal moving margin, again verify that it's at least just as wrong as the other models. You're going to have to come up with something to deal with the negative values of thickness that you'll get...

Shallow shelf approximation

Now, consider the equations describing a flow that is vertically integrated. The equations are

\frac{\partial}{\partial x}\left ( 2 \eta H 
\left(2\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)\right)
+\frac{\partial}{\partial y}\left(\eta H\left(
\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\right)
=\rho gH \frac{\partial s}{\partial x}

\frac{\partial}{\partial y}\left ( 2 \eta H 
\left(2\frac{\partial v}{\partial y}+\frac{\partial u}{\partial x}\right)\right)
+\frac{\partial}{\partial x}\left(\eta H\left(
\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\right)
=\rho gH \frac{\partial s}{\partial y}

Exercises

As a first exercise in solving these equations, try the experiments described in the EISMINT ice shelf models, but never published [1]. Get the self-descr.pdf, or the first hyperlink on the page. Let's do experiments 3-4 on page 6 of the document (Note that we will see and work with the solution to these experiments again when we do some exercises with the higher-order dynamics routines in Glimmer/CISM).

References

  1. Huybrechts et al. The EISMINT Benchmarks for Testing Ice--Sheet Models. Ann. Glaciol. (1996) vol. 23 pp. 1-12 pdf