How to solve an evolutive-advection-diffusion problem with source term (subglacial hydrology)?

[antoniojm]

Hi there! I am trying to solve a problem in relation with subglacial hydrology (https://gmd.copernicus.org/articles/8/1613/2015/gmd-8-1613-2015.pdf) which reads as follows:

$$\frac{\partial W}{\partial t}=-\nabla\cdot (VW)+\nabla\cdot(D \nabla W) +m/\rho_w$$

with the vector $V=-k\nabla(P+\rho_w gb)$ and $D=\rho_w gkW$. The term $P+\rho_w gb$ is called the hydraulic potential ($\phi$).

My aim is to model a little glacier with water underneath. W is the water thickness in the layer below ice and is the variable to solve. Depending on the pressure (P) and the topography of the bedrock (b) in which the glacier is, water circulates in a way or another, accumulating in some regions and disappearing from others.

For the variables above, P is the pressure defined in a simple case as $P=P_o=\rho_i g H$ (overburden pressure of the ice) with H the ice thickness above the water layer. In turn, H can be modeled as the difference H=s-b with s the surface elevation of the glacier (normally a square-root function) and b is the bedrock elevation at the bottom of the glacier (for this, my plan is to define a function that tries to represent some mountainous terrain, but for the moment it can be flat). Regarding the parameters, $\rho_w$ is the density of water, $k$ is conductivity and $m$ is the source term which is for my case a small constant number. I have defined all this functions in Python in order to make further modifications.

Regarding the equations, as you can see, there is an advection term, together with a diffusion term and in addition a source term. Do you know how to solve it correctly with FiPy?

I have tried the following code to solve it, but I think that the convection term is producing some issues, because the diffusion and source terms work well separately. Also, I don't know if the way I have tried to implement the gradient in V, with functions previously defined (not shown below) is the correct way. As boundary conditions I am trying to impose W=0 at the edges of the domain.
Here there are the functions I use for defining the variables I need with some cases to choose:

k=0.001
rhow=1000
m=0.0001
rhoi=910

def surface(nx,ny,case):
    zs=np.zeros((nx,ny))
    for i in range(nx):
        for j in range(ny):
          if case==0: zs[i,j]=600
          if case==1: zs[i,j]=6*((i*400+1000)**0.5-1000**0.5)
    return zs

def hydraulic_potential(P,b):
    phi=P+rhow*g*b
    return phi

def grad_hyd_potential(phi,nx,ny,dx,dy):
     phidx=np.zeros((nx,ny))
     phidy=np.zeros((nx,ny))
     for i in range(1,nx-1):
        for j in range(1,ny-1):
            phidx[i,j]=(phi[i+1,j]-phi[i,j])/dx
            phidy[i,j]=(phi[i,j+1]-phi[i,j])/dy  
     phidx=phidx.reshape(nx*ny)
     phidy=phidy.reshape(nx*ny)
               
     return phidx,phidy

def bedrock(nx,ny,case):
    b=np.zeros((nx,ny))
    for i in range(nx):
        for j in range(ny):
          if case==0: b[i,j]=0
          if case==1:
              s=4
              m=0
              b[i,j] = 100*np.exp(-( (np.sqrt((i-nx/2)**2+(j-ny/2)**2)-m)**2 / (2.0 * s**2 ))  )
       
    return b

def calc_pressure(nx,ny,n,W,P0,mode):
    P=np.zeros((nx,ny))
    if mode=='P0':
        P=np.copy(P0)
    if mode=='P2a':
        P=min(rhow*g*W,P0)
    if mode=='P2b':
        P=min(P0*(W/0.1)**3.5,P0)
    return P

And here is the code that I have tried to solve it:

Lx, Ly = 2500, 2500
Nx, Ny = 50, 50
dx, dy = Lx/Nx, Ly/Ny


Np = 20    
Nt = 500
dt=100


# Define the grid/mesh
mesh = Grid2D(nx=Nx, ny=Ny, dx=dx, dy=dy)
x, y = mesh.cellCenters[0], mesh.cellCenters[1]



# Define the model variable and set the boundary conditions
W = CellVariable(name="Water height", mesh=mesh, value=0.5)
b = CellVariable(name="Bedrock", mesh=mesh, value=numerix.array(bedrock(Nx,Ny,1)))
s = CellVariable(name="Surface", mesh=mesh, value=numerix.array(surface(Nx,Ny,1)))
H = CellVariable(name="Water height", mesh=mesh, value=numerix.array((s-b)))
P0 = CellVariable(name="Overburden pressure", mesh=mesh, value=numerix.array(rhoi*g*H))
P = CellVariable(name="Water pressure", mesh=mesh, value=numerix.array(calc_pressure(Nx,Ny,0,W,P0,'P0')))
hyd_pot=CellVariable(name="hydraulic potential", mesh=mesh, value=numerix.array(hydraulic_potential(P,b)))


X, Y = mesh.faceCenters

meshBnd = mesh.facesLeft | mesh.facesTop | mesh.facesRight | mesh.facesBottom

# Functions for bedrock, surface elevation, ice thickness, etc. For defining the gradient in the hydraulic potential.
bv=bedrock(Nx,Ny,1)
sv=surface(Nx,Ny,1)
Hv=np.maximum(sv-bv,0)
P0v=rhoi*g*Hv
Pv=calc_pressure(Nx,Ny,0,W,P0v,'P0')

phi=hydraulic_potential(Pv,bv)
phidx,phidy=grad_hyd_potential(phi,Nx,Ny,dx,dy)

# Solving the problem
eq = TransientTerm() == DiffusionTerm(coeff=rhow*g*k*W)\
     + ExponentialConvectionTerm(CellVariable(mesh=mesh,value=[k*phidx,k*phidy])) \
     + m/rhow

W.constrain(0., where=meshBnd)

# To save the solution in each time step and plot after that
my_sol = np.zeros((Np,Nx*Ny))      
my_sol[0,:] = W
k = 1

for step in np.arange(1,Nt):
    #print(step)
    eq.solve(var=W,dt=dt)
    if np.mod(step,Nt/Np)==0:
        my_sol[k,:] = W
        k += 1

[guyer]

Sorry for the delay; I started to compose answer several days ago, but I'm not sure what happened to it.

I have tried the following code to solve it, but I think that the convection term is producing some issues

What issues? What do you see? What do you expect to see? Once I import the necessary classes and modules, your code runs, but I have no idea what you expect it to do.

• W is the only field that evolves?
• ExponentialConvectionTerm(coeff=k * phi.faceGrad) is better
• Writing things like P = CellVariable(name="Water pressure", mesh=mesh, value=numerix.array(calc_pressure(Nx,Ny,0,W,P0,'P0'))) breaks any dynamic connection between W and P as W evolves. Better to write something like

def calc_pressure(n,W,P0,mode):
    if mode=='P0':
        P=P0.copy()
    if mode=='P2a':
        altP = rhow*g*W
        P = altP * (altP <= P0) + P0 * (altP > P0)
    if mode=='P2b':
        altP = P0*(W/0.1)**3.5
        P = altP * (altP <= P0) + P0 * (altP > P0)
    return P

and then just

H = s-b
H.name = "Water height"
P0 = rhoi*g*H
P0.name = "Overburden pressure"
P = calc_pressure(0,W,P0,'P0')
P.name = "Water pressure"

[antoniojm]

Thanks for your answer @guyer, I have tried something similar these days, and changed the coefficient in the convective term and it worked well.

Answering your questions:

• Yes, W is the only term that evolves, although P could also evolve in a next step of my implementation together with others variables like H. In that case, I suppose that I would have to update the variable in the loop I use before solving the defined problem.
• I was expecting to see a movement of the water toards regions with lower pressure and where there are holes of the bedrock elevation, and now with the corrections implemented is what I get.

Thank you again for your comment. I have recently discovered FiPy and it was my first time using it. I will continue developing my code.

[guyer]

• Yes, W is the only term that evolves, although P could also evolve in a next step of my implementation together with others variables like H. In that case, I suppose that I would have to update the variable in the loop I use before solving the defined problem.

If P is a simple function of W, like your "P2a" and "P2b" cases, then you don't need to do anything. P will automatically update whenever W changes if you write it the way I show.

If P will have its own evolution equation, then you will need to either solve the W and P equations sequentially or you will need to couple them and solve the coupled equation. Either way, you will probably need sweeps to get the solutions to converge.

Be aware that the more your equations start to look like Navier-Stokes flow equations, the coupling between pressure and velocity can get tricky. If you find yourself in that situation, see, e.g., examples.flow.stokesCavity for how to treat the pressure correction.

[antoniojm]

Ok. Thank you again. Good to know because there is also the possibility in he paper I worte down in my first message to get an own equation for P evolving with time. Thanks for the links too.
Best regards, Antonio.