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