Implementation of a prognostic free surface problem

In prognostic free surface problems, the kinematic boundary condition for the free surface is a partial differential equation (PDE) and hence has to be solved. The domain of the solution, though, is a “(problem dimension - 1)” entity. I.e., in case of a 3D glacier, the domain is a surface, in case of a flow-line model being applied, the domain is a line.

To this end, Elmer has a feature that enables a boundary to be declared as an additional body. The syntax is the following:

Boundary Condition 2
  Name = "surf"
  Target Boundaries = 2
  ! this is the assigned body number of the free surface domain
  Body ID = 2
  ...
End

This declares the free surface to act as the second body in the simulation (if only one 3D body is declared). For this body, separate sections (body force, initial conditions, …) have to be declared in a separate body declaration:

Body 2
  Name= "free surface"
  Equation = 2
  Material = 2
  Body Force = 2
  Initial Condition = 2
End

The different sections attached to this body then may read as:

! NB:  units in meters MPa years
!-------------------------------
Equation 2
 Active Solvers(1) = 5 !(the Solver ID of the FreeSurfaceSolver)
 Flow Solution Name = String "Flow Solution"
 Convection = Computed
End

Initial Condition 2
! this sets the initial position of the free surface variable
 FS = Equals Coordinate 3
! and this the reference to this initial condition
 ReferenceFS = Equals Coordinate 3
! this sets a zero deviation from the input mesh
 Mesh Update 3 = Real 0.0
 Mesh Update 2 = Real 0.0
 Mesh Update 1 = Real 0.0
End

Body Force 2
  ! the accumulation flux in the 3 Cartesian coordinates
  FS Accumulation Flux 1 = Real 0.0e0
  FS Accumulation Flux 2 = Real 0.0e0
  ! 1.0 m/a w.e. converted into ice equivalent)
  FS Accumulation Flux 3 = Real $ 1.0 * 1000.0/910.0
End

Material 2
 Density =  Real MATC "910.0*1.0E-06*(31556926.0)^(-2.0)"
 ! the minimum value of the free surface variable
 ! keeping a minimum flow depth of 10 meters
 Min FS = Variable Coordinate 3, Height
      Real MATC "tx(0) - tx(1) + 10.0"
 ! allow a maximum thickening of 100 meters with
 ! respect to the initial position of the free surface
 Max FS = Variable ReferenceFS
      Real MATC "tx(0) + 100.0"
End

There are at least two solvers needed, namely, the free surface solver itself and the mesh update solver, in order to shift the nodes of the mesh to comply with the altered geometry imposed by the change in the free surface. Alternatively, if you are using a structured mesh, you might consider replacing the MeshUpdate solver by the StructuredMeshMapper solver'', as explained here.

! the free surface solver to be solved on Body 2
!-----------------------------------------------
Solver 5
  ! usually, the solver is executed only after the thermo-mechanical
  ! problem has obtained a solution on the time-level
  Exec Solver = "After TimeStep"

  Equation =  String "Free Surface Evolution"
  ! the name of the variable
  Variable = "FS"
  Variable DOFs = 1

  Procedure = "FreeSurfaceSolver" "FreeSurfaceSolver"

  ! this enables the limitation of the free surface
  ! by upper and/or lower limits (see material section above)
  ! using a variational inequality formulation
  Apply Dirichlet = Logical True

  Linear System Solver = Iterative
  Linear System Iterative Method = BiCGStab
  Linear System Max Iterations  = 1000
  Linear System Preconditioning = ILU1
  Linear System Convergence Tolerance = 1.0e-08

  Nonlinear System Max Iterations = 100 ! variational inequality needs more than one round
  Nonlinear System Min Iterations = 2
  Nonlinear System Convergence Tolerance = 1.0e-06

  Steady State Convergence Tolerance = 1.0e-4
  ! switch that off in parallel runs, as it may introduce
  ! partition dependent relaxation factors:
  ! Maximum Displacement = Real 10.0

  Stabilization Method = Bubbles
   Flow Solution Name = String "Flow Solution"

  ! this is needed if the variational inequality method
  ! (Apply Dirichlet = Logical True) is applied
  Exported Variable 1 =  FS Residual
  Exported Variable 1 DOFS = 1

  ! this variable contains the free surface's initial
  ! values (set in initial condition above and needed
  ! for limiting maximum changes as well as in the boundary condition
  ! at the free surface)
  Exported Variable 2 = ReferenceFS
  Exported Variable 2 DOFS = 1
End

! the mesh update solver to be solved on body 1 (the main ice volume)
!--------------------------------------------------------------------
Solver 6
  Equation = "Mesh Update"
  ! usually, the solver is executed only after the thermo-mechanical
  ! problem has obtained a solution on the time-level
  Exec Solver = "After TimeStep"

  Linear System Solver = Iterative
  Linear System Iterative Method = BiCGStab
  Linear System Max Iterations  = 500
  Linear System Preconditioning = ILU1
  Linear System Convergence Tolerance = 1.0e-06

  Nonlinear System Max Iterations = 1
  Nonlinear System Convergence Tolerance = 1.0e-06
End

The material section entry for the mesh update solver then read as (have to be inserted in the glacial body's material section and values are not so important):

Mesh Elastic Modulus = 1.0
Mesh Poisson Ratio = 0.3

(Note: Youngs Modulus/Poisson Ratio are old keywords, which still work, but the new ones are recommended to be used.)

If the model consists of a bedrock boundary, a free surface and a side boundary, the boundary condition section read as follows:

! the bedrock boundary
!----------------------
Boundary Condition 3
  Name = "bedrock"
  Target Boundaries  = 1

  ! zero flow height at bedrock
  Height = Real 0.0

  ! no-slip condition
   Velocity 1 = 0.0
   Velocity 2 = 0.0
   Velocity 3 = 0.0

  ! keep the nodes at the bedrock, as they are
   ! NB.: iso-static rebound could  be linked in here
   ! by linking the change of the bedrock position
   ! to Mesh Update 3, similar as for the free surface
   ! at the upper surface boundary
   Mesh Update 1 = Real 0.0
   Mesh Update 2 = Real 0.0
   Mesh Update 3 = Real 0.0
End

! the free surface boundary
!--------------------------
Boundary Condition 2
  Name = "surf"
  Target Boundaries = 2

  ! this is the assigned body number of the free surface domain
  Body ID = 2

  ! zero flow-depth at the free surface
  Depth = Real 0.0

  Mesh Update 1 = Real 0.0
  Mesh Update 2 = Real 0.0

  ! links the free surface to the mesh update variable
  Mesh Update 3 = Variable FS, ReferenceFS
      Real MATC "tx(0) - tx(1)"
End

! the side wall
!--------------
Boundary Condition 3
  Name = "side"
  Target Boundaries  = 3
  ! free slip, no penetration, hydrostatic pressure
  ! -----------------------------------------------
  Normal-Tangential Velocity = True
  Velocity 1 = 0
  Velocity 3 = 0
  Flow Force BC = True
  Normal Pressure = Variable Depth
         Real MATC "tx*910.0*1.0E-06" ! in mPa

  ! glacier boundary remains unchanged
  Mesh Update 1 = Real 0.0
  Mesh Update 2 = Real 0.0
  Mesh Update 3 = Real 0.0
  ! freeze the free surface
  FS = Equals ReferenceFS
End

Additional two instances of the FlowDepthSolver might be needed to inquire the flow depth/height of the glacier (a non-trivial problem on an unstructured FEM-mesh).
“Mass Consistent Normals = Logical True” might be needed in the boundary conditions to ensure correct calculation of fluxes (through bedrock and lateral sides).