Positive `H` condition for finite differences scheme in SIA equation

[facusapienza21]

We observed that the total ice mass of a glacier is not preserved during SIA simulations without mass balance component. We further observe that the ice seems to move "uphill" and spread in regions where there should not be ice.

The following simulations are for the period (2014, 2018).

This is the initial shape of the glaciers before running the simulation

Without imposing any boundary condition or condition in SIA inside iceflow.jl, the total mass of the ice when we run the forward model without mass balance terms increases (here the percentage over 100% of mass gain/loss by the glacier): We can observe that for some glaciers this is quite significant and leads to unrealistic solutions. The increase in mass of the glacier is due to the positive truncation happening here:

function iceflow!(dH, H, context, t)
    @timeit to "iceflow! PDE" begin
    # First, enforce values to be positive
    H[H.<0.0] .= H[H.<0.0] .* 0.0

    # Compute the Shallow Ice Approximation in a staggered grid
    SIA!(dH, H, context)
    end
end    

The problem needs to be solved directly in the solver, by making solutions with H < 0 non posible. We can also see how this leads to glacier that seems to expand uphill:

This behavior can be partially fixed by making a cap of the solution for very small values of the ice thickness:

ϵ = 0.0001
 H[H.<ϵ] .= 0.0

and then we obtain

However, the problem still came from the lack of imposing the condition H >= 0 inside the solver. We can solve this by adding the following cap inside SIA():

  η₀ = 1.0
  dSdx_edges .= min.(dSdx_edges,  η₀ * H[1:end-1, 2:end-1]./Δx,  η₀ * H[2:end, 2:end-1]./Δx)
  dSdy_edges .= min.(dSdy_edges,  η₀ * H[2:end-1, 1:end-1]./Δy,  η₀ * H[2:end-1, 2:end]./Δy) 
  dSdx_edges .= max.(dSdx_edges, -η₀ * H[1:end-1, 2:end-1]./Δx, -η₀ * H[2:end, 2:end-1]./Δx)
  dSdy_edges .= max.(dSdy_edges, -η₀ * H[2:end-1, 1:end-1]./Δy, -η₀ * H[2:end-1, 2:end]./Δy)

This correction allows to bound cases where the surface slop is to high and make that the ice thickness in the current grid point goes below zero. After making this fix, we can see that the change in total ice is almost zero: and changes in ice thickness became

This nows looks much better to me that what we had before. I still see the glacier expanding a little bit outside the contours, but I don't think this is unrealistic. @JordiBolibar do you have any thoughts on this?

[facusapienza21]

As a note, I still believe is worth exploring which is the best way of imposing the positive condition in the solver in order not to loss mass. Right now, this condition is strictly sufficient but not necessary and I don't see how this can be use for cases where the ice thickness in a given location needs to strictly reach zero.

I also believe the condition can be mathematically simplified to one single line. Right now, I am doing the cap both as a lower and upper bound, but I am not sure this is necessary, but seems to be working.

[facusapienza21]

After revisiting this in the context of the Halfar solutions, I figured that the right way of implementing the boundary condition is actually

  dSdx_edges .= @views @. min(dSdx_edges,  η₀ * H[2:end, 2:end-1]/Δx)
  dSdx_edges .= @views @. max(dSdx_edges, -η₀ * H[1:end-1, 2:end-1]/Δx)
  dSdy_edges .= @views @. min(dSdy_edges,  η₀ * H[2:end-1, 2:end]/Δy)
  dSdy_edges .= @views @. max(dSdy_edges, -η₀ * H[2:end-1, 1:end-1]/Δy)

This imposes the mathematical condition just where is necessary, by imposing

- \frac{H_{i, j+1}}{\Delta x} \leq [dSdx]_{ij} \leq \frac{H_{i+1, j+1}}{\Delta x} 

However, this boundary condition still falis to resolve the right solution in the borders of the glacier. Here we can see a simulation of the Halfar solution where we can observe the missfit between the theorerical solution and the one obtained by the solver in the border.

I expect this way of implementing the border solution can fail for certain glaciers, specially in the border.