diffusion for multiple layers

[IgnaceP]

Hello,

I needed multiple sediment layers (a compacted and a freshly deposited sediment layer with varying critical shear stress and density). Therefore, I needed to update the diffuse_bottom method in telemac.py.

I decided to conserve the elevation of the layer boundaries. If the diffusion result in local erosion, the top layer is eroded untill the layer beneath is reached which will further erode, etc. In case of local deposition, the newly deposited sediments will have the density of the top layer. See image below:

I am not sure whether this is a good approach, I believe I am violating conservation of mass (erosion of compacted sediments are deposited nearby as freshly deposited sediment with a lower density). What are your thoughts on this?

The updated code:

def diffuse_bottom(self, nu, dt, t, nproc = 1, launcher = 'mpiexec',
                    step = -1):
     """Smooth the bottom surface by diffusion.

     Args:
         nu (float): Diffusion coefficient (m^2/yr).
         dt (float): Time step (yr).
         t (float): Duration of diffusion (yr).
         nproc (int, optional): Number of MPI processes. Default to 1 (no
             MPI).
         step (int, optional): Telemac time step on which the diffusion is
             applied. Default to -1 for last time step.

     Todo:
         Update for more than one mud layer/class (e.g., different rho
         values).

     Bug:
         Use of serial mode (nproc = 1) makes next Telemac run (whether in
         serial or parallel mode) to crash (observed on the SCC at Boston
         University, never tested elsewhere so far).

     """
     # Class attributes.
     x = self.x
     y = self.y
     tri = self.tri
     rho = self.rho
     u0 = self.u[step, :]
     v0 = self.v[step, :]
     s0 = self.s[step, :]
     b0 = self.b[step, :]
     r0 = self.r[step, :]
     t0 = self.t[0, step, :]
     nl = self.m.shape[0]

     # preserve elevation of layer boundary layers
     rb = np.zeros((nl, 1, 1, self.npoin))
     for i in range(nl):
         if i == 0:
             rb[i, 0, step, :] = b0 - (self.m[i, 0, step, :] / rho[i])
         elif i == nl-1:
             rb[i, 0, step, :] = r0
         else:
             rb[i, 0, step, :] = rb[i - 1, 0, step, :] - (self.m[i, 0, step, :] / rho[i])

     # Diffuse bottom surface.
     # Attention: t is time, not the Telemac variable for cohesive sediments.
     if nproc <= 1:

         ################
         # Serial mode. #
         ################

         # Call diffusion function.
         from demeter import telemac_diffusion
         bi = telemac_diffusion.diffusion(x, y, b0, tri, nu, dt, t)

     else:

         ##################
         # Parallel mode. #
         ##################

         # Create directory to store intermediate input files.
         if not os.path.isdir('./tmp_diffusion'):
             os.mkdir('./tmp_diffusion')

         # Intermediate file names.
         x_global_fn = './tmp_diffusion/x_global.txt'
         y_global_fn = './tmp_diffusion/y_global.txt'
         f_global_fn = './tmp_diffusion/f_global.txt'
         tri_global_fn = './tmp_diffusion/tri_global.txt'
         f1_global_fn = './tmp_diffusion/f1_global.txt'

         # Save intermediate files.
         if not os.path.isfile(x_global_fn):
             np.savetxt(x_global_fn, x)
         if not os.path.isfile(y_global_fn):
             np.savetxt(y_global_fn, y)
         np.savetxt(f_global_fn, b0)
         if not os.path.isfile(tri_global_fn):
             np.savetxt(tri_global_fn, tri, fmt = '%d')

         # Run parallel diffusion module.
         os.system(launcher + ' -n %d python ' % nproc +
                   '$DEMPATH/telemac_diffusion.py %f %f %f' % (nu, dt, t))

         # Load intermediate file.
         bi = np.loadtxt(f1_global_fn)

     # Treat the rigid bed.
     b1 = np.maximum(bi, r0)

     # Update other variables.
     s1 = np.maximum(s0, b1)
     u1 = u0 * (s0 - b0) / np.maximum(s1 - b1, 1e-3)
     v1 = v0 * (s0 - b0) / np.maximum(s1 - b1, 1e-3)
     t1 = t0 * (s0 - b0) / np.maximum(s1 - b1, 1e-3)
     u1[s1 - b1 <= 1e-3] = 0
     v1[s1 - b1 <= 1e-3] = 0
     t1[s1 - b1 <= 1e-3] = 0

     # update mudd mass variables
     for i in range(nl):
         rb[i,0,0,:] = np.minimum(rb[i,0,0,:], b1)
         if i == 0:
             m = (b1 - rb[i,0,0,:]) * rho[i]
         else:
             m = (rb[i-1,0,0,:] - rb[i,0,0,:]) * rho[i]
         m[m < 0] = 0
         self.m[i, 0, step, :] = m

     # Update class attributes.
     self.u[step, :] = u1
     self.v[step, :] = v1
     self.s[step, :] = s1
     self.b[step, :] = b1
     self.t[0, step, :] = t1

[ogourgue]

That's a tricky one!

Ideally, you should not have transfer of mass from one layer to the other in the process. If I understand well, your top surface is diffused, while the interface between the two layers is fixed.

In Demeter v1 (https://github.com/ogourgue/demeter-v1/blob/master/python/demeter.py - line 1684), it seems that I was diffusing "delta z" (that is, the difference in top surface elevation before and after the Telemac run). That is the old HPP version (https://doi.org/10.5194/esurf-10-531-2022). That is not ideal if we want to link this process to soil creeping (https://doi.org/10.1016/j.ecss.2018.05.016).

While developing Demeter v2, I remember considering diffusing the layers themselves separately (variable "m", I would guess) and then recalculate the bathymetry. But I'm not sure it is better.

So, in my humble opinion, I don't think there is any "good" solution :-) So it might be good to consider the possibility of having different options (even if only one is implemented first).

I would suggest to add an option in the function (e.g., parameter keyword), so we can choose between options. Then we just implement the one you propose here and we document that there might be other options needed int he future. Then we also need a new example with two layers.

If you are OK with that, we can create an issue and a new branch where you implement your solution. Then I add some context in comments for future work. Can you also provide the example, similarly to those already available? I can help with that if need as well.

What do you think?