Issues with settling fluxes

[fdufois]

Hi all,

I recently noticed a weird behavior in my ROMS/sediment config.
I have several sand to gravel sediment classes, and all the bigger sizes (with settling velocities ranging from 2cm/s to 25 cm/s) behave similarly in the water column (i.e. the relative sediment profiles are the same, meaning that all classes are falling at the same rate). This is obviously wrong (cf. Rousse profile), and I'm trying to figure out what's going on.

To make it clearer and reproducible by others, I describe below what can be observed when running the 1DV sedbed_toy test case. I ran the model for 20m depth, for 4 different sediment sizes (with settling velocities of 1,5, 10 and 40 cm/s) and for time steps ranging from 5 to 160s. I show in the figures the relative concentration profile for a specific time (when there is a bit of resuspension going on).

There are two main behaviors that I'd like to highlight:
(i) the sediment relative profile (C/C(k=1)) tend to be similar for larger settling velocities, even more as dt increases (fig 1). In figure 1 each panel is for a specific DT, and the different lines represent the different settling velocities.

(ii) the sediment relative profile (C/C(k=1)) is time-step dependent (fig2), even more as the settling velocity increases. In figure 2 each panel is for a specific settling velocity, and the different lines represent the different DT.

For those larger sediment classes there shouldn't be much suspended load most of the time, and the advection/diffusion scheme should be able to deal with that easily but does not seem to handle it.
I assumed that the issue comes from the sed_settling.F routine imposing the settling fluxes. The scheme used for the settling fluxes doesn't have a CFL, but seems to be higly time-step and settling velocity dependent, which is an issue.

Is there another scheme that can be used for the sediment settling fluxes?

Hope this is clear enough, otherwise I can add some more details.

Cheers,

Francois

[jcwarner]

I haven't looked at this in a while, but i can dig back into it.
I suppose I am not clear what you plotted in fig 1. You said that each line is a different wsed, but then what are the different panels? One thing to check is that the bottom stress gets to be the same, especially as you change dt. You need to run these types of tests until they reach steady state. Are you running these all at the same time (ie with 6 sed classes at the same time) or are you running this with 6 different runs, one sed class at a time?
-john

[fdufois]

Hi John,
In the results I showed the model reached a steady state. I actually showed the results at the end of the 1 N/m2 wind stress step, and the concentrations weren't changing anymore at this stage for all the 4 sed classes.

I run the 4 sed classes at the same time, and repeat this run for 6 different time steps. I double-checked and the bottom stresses are identical in all the runs, reaching 1N/m2.

In figure 1, each line is a different wsed (4 lines), and each panel is for a different DT (6 panels for 6 different dt indicated above each panel). It's a bit confusing, and sorry for that, but I just noticed that fig 1 is the one title fig1_sedtoy.png and is actually the 2nde figure attached to the post.

I did a few checks and my gut feeling is that there is nothing wrong with the settling flux scheme itself. I printed out the concentration after that step, and for the fastest settling sediments everything has been removed from the water column as expected. The issue probably comes from the fact that the scheme along the vertical is explicit, and erosion and diffusion come by steps after the settling step. I'm wondering if in case of an explicit scheme, it wouldn't make sense to have settling last, otherwise the diffusion keeps on adding up stuff that should have been removed by the massive settling fluxes? Does that make sense or am I messing up everything?

Cheers,

Francois

[jcwarner]

so the sedbed_toy might have several variations over the years. the one i have ramps the forcing up and down to look at several cycles of resuspension and deposition. So to be sure that we are using the same cases, can you send the .h ocean.in and any ana files to me? jcwarner@usgs.gov
thanks
john

[jcwarner]

sorry for taking a long time to respond. I received your email and was able to run the test case. as you said it was basically a shorter version of the distributed case. you can get different solutions depending on how many sed classes you choose, but the effect you are talking about is still there. I think this all due to the choices made, and it identifies issues that users should be aware of. Here are some comments.

- The vertical balance is essentially from 3 terms: vertical mixing, settling, and erosion flux.
-The vertical mixing is implicit and not time step limited.
- The settling is performed by a WENO scheme that reconstructs the vertical profile and advects that profile downward a distance of ws*dt. The amount that is advected out of the bottom is placed on the seafloor and is available for resuspension.
- The erosive flux is basically Erate*dt*a_bunch_of_other_stuff (poros, density, sed frac, etc).
- For some of the tests that you have, ws*dt is 40cm/s*160s=6400cm=64.00m. That distance of settling is greater than the total water depth and clears the water column each time step. Then a hope that the erosive flux will be diffused upwards in the same step to give a nice profile. Even though the model is stable, i think this is pushing the practical limits of the physics. There needs to be a lower practical limit for the time step.
- For the lower settling velocity, we have ws*dt is 1cm/s*160s=160cm=1.6m. That is probably more than the bottom cell thickness, probably into the second layer. As you can see those profiles (different dts for ws=1cm/s) do not vary much for all those time steps. This is because the settling is only removing a smaller total fraction.
-I think the Figure 1 (actually the second figure) is just the same thing but looking at it based on time step, not ws.
- For the coarser material, i recommend a smaller dt. I am not sure what that really should be, i was thinking that ws*dt<0.5*Hz(1). but that may be too restrictive.
- I would be happy to run a few more tests, but i think it is working as it should. We just need to be careful to allow the physics to evolve for the time step.

-john

[fdufois]

Hi John,

Thanks for looking at this!
I agree with your comments, but I really think this is an issue when doing sediment transport modelling at the regional scale. I guess that the way I plotted the outputs didn't give justice to the issue. That's probably because I was mostly showing the profiles (C/Cref) and not the actual absolute values.

In the figure below I show this time the concentration (each pannel is a different Ws, and each line a different dt) for the same configuration as before (20m depth, settling velocities of 1,5, 10 and 40 cm/s). And we see massive differences in the concentrations we get, and that whatever the settling velocities (after a certain threshold) the results are the same

NB: the units are in mg/l everywhere not kg/m3 as stated in the legend...

And to emphasis the issue I used smaller particles (settling velocities), and also smaller water depth (2m). Each sigma level is roughly 0.1m in this config. That range of settling velocities is for floculated mud to medium sand (~500 um), so nothing to big.
Roughly we see that if you exceed the CFL along the vertical it starts diverging.

That's typically what's gonna happen in my regional config, and I can't really go down to 10s time step to be able to correctly simulates the fine to medium sand dynamics. That would really be too restrictive for shelf scale applications. And also with the sigma level I use in my regional config (refined towards the bottom), I expect to have big biases in a lot of shallow regions, and those biases are gonna be depth dependent because of the sigma discretisation (time steps, settling velocities, or layer thickness all impact the results the same way).

I tried to do a few different changes within the code, in order to change the order of the erosion, settling and vertical advection/diffusion processes.
None of my tests have been really satisfying, but I ended up with something which I find a bit better.

Now I do: 1-erosion, 2-vertical adv/diff and 3- settling. Because the WENO scheme is not positively defined, I used a simple upstream scheme with a sub time-stepping (dtt=0.1s).
The concentration profiles are still time-step dependent but at least for the larger time steps I don't end up with floating gravels. And as well for a given dt the bigger the settling velocity, the smaller the concentration, which was a critical issue with the former scheme.

Below are the results with this new scheme for the 20m depth case, with coarse particles...

Below the plot is for the 2m depth case, with smaller settling velocities...

I'm now convinced that if we don't respect the CFL there is no way we can have the right profiles (C/Cref) doing erosion, diffusion and settling by steps. Indeed putting the settling fluxes last doesn't solve this issue at all, but at least when it can't solve properly the suspended transport for bigger particles it tends to underestimate it. And previously my only option was to remove suspended load for all the sands, because I can't risk having floating heavy gravels in the shallow regions. The way it is now, I don't have to switch off suspension for sands, and where it is far from respecting the CFL it just gonna underestimate it.

I'm left wondering if it is not possible to just add the selling velocity to W in the diffusion/advection implicit step? Would that spoil everything?

Cheers,

Francois

[fdufois]

I updated the last 2 figures of the post above.
I indeed had an issue because the scheme wasn't conserving mass.
I double-checked and everything should be OK now, although what I did probably only works for that specific case with its specific CPP keys.

I had an issue because I got tricked by step3d_t. The concentration unit changes from m.Tunit to Tunit within the step3d_t routine. And also for some reasons I couldn't plug the settling after the step3d_t routine(mass was not conserved). I couldn't figure out why, but I ended up putting the settling inside step3d_t just after the diffusion step.

I attached the modified routines. Please let me know if you feel that there is anything dodgy in what I did.

Cheers,

Francois