open TELEMAC-MASCARET The mathematically superior suite of solvers

Dear All,

the issue is connected to two previous posts: Bed celerity -against- the directions of flow and sediment transport and Sisyphe in transcritical flows.

We have come across an ugly situation in Sisyphe coupled with Telemac-2D. The situation is easy to explain: Single nodes gets dry in the middle of the domain mostly due to bed load accumulation and not due to the water level changes. There is no velocity on the node and the water depth is zero. However, the bed load transport does not stop and the bottom -- now seemingly a dry land -- grows steadily up and up... Looks ugly and shows that we have some principal numerical problem here.

The pattern is easy to recognise -- these single nodes are surrounded by nodes with non-zero velocity in a patch, and this pattern remains as it is during many time steps of the simulation because of the topography and resolution. These nodes can also appear as a frontal nodes (against the current carrying the bed load) of a line of dried-out nodes forming a "dry edges" line.

A naive explanation could be that the sediment still gets into a volume associated with the "dried-out" node, but never gets out, so that the accumulation takes place.

I try to attach a few illustrations explaining the situation: A strongly exaggerated view, a plot showing that although on a node the water depth is zero, the bed load is not, and two pictures with the characteristical element patches around the nodes in quetion (the color represents elevation).









Well, of course the first thought would be to refine the area, so that maybe the drying out can take place not only node-wise, but also element-wise, compatible with Telemac-2D wetting/dryinge approach. The second possibility would be to refine in such way, that the situation mentioned above does not appear in terms of the sediment transport. Both wetting/drying and positive depths algorithms are extremely mesh-dependent, so maybe taking care that edges of elements in the problematic places are not parallel to the global flow direction could help. Also: Taking care that the hydrodynamics is stabilised (there are some badly treated local transcritical flows there).

But, anyway, this seems to me not to be a small bug, but a principle fault in the transport routine of Sisyphe, which seems to have trouble by modelling of drying out in the middle of the domain due to the accumulation of the material brought by the bed load -- and not due to the water level changes.

What are your comments, practical recommendations and solutions? Have you come across such situations? How the algorithm could be repaired / improved?

All comments welcome!

Best regards,
jaj

Dear Jacek,
Many thanks again for your comments and feedback and for point us this. We have been discussing about the possibility of expand the number of numerical schemes available to solve sediment transport problems. Very soon we have a consortium meeting and this will be certainly a topic to be discussed.

Pablo

Dear Pablo,

come on, do not give up so quickly! And, by all respect, such problems will not be solved by any "consortium meetings", but by clever people like you!

And maybe it is better to improve an existing algorithm than to "expand the number of algorithms"... I prefer only one option working robustly than a choice of ten fancy things which do not at all deliver what the people really want.

It seems to me that this (relatively) new positive depths algorithm relies on some upwinding method which remains blindly one-directionally upwinded even when the situation does not actually fit to the upwind, isn't it so?

And, when it comes to the "number of available algorithms": Where is the old FEM algorithm gone? Deleted? Just in such situation one would like to fall back into a previously existing option which -- maybe -- would have been more compatible with Telemac-2D "tidal flats" algorithms...

So, please think about it -- maybe there are some recommendations how to deal with the problem without necessity to re-develop the algorithms?

Looking forward for your answer,
jaj

Dear jaj,

I have no solution but an answer to your second question, giving an example where I have a similar observation.

A long-term sediment transport (only susp.load, no bed load) simulation (1-year) gives deposition heights of some meters in boundary nodes which should be dry during the whole simulation and hence they also should be charged with zero deposition.

Testing the case by means of a relatively simple and short hydrodynamic calculation to a steady state (4 hours), Telemac-2D calculates water depths in boundary nodes which should by dry since they the bottom elevations there are higher than the expected water surface elevations. Enclosed a picture to visualize the situation. The node numbering (time series of water depths) within the polygon is from left to right. Interestingly the second node has exactly zero water depth which should be the case for all the boundary nodes. Why this node? Mesh dependency?

Now if you simulate sediment transport for one year based on this configuration these errors accumulates drastically..

Below also my standard settings.

Best regards,
Clemens





TREATMENT OF NEGATIVE DEPTHS =2
MASS-LUMPING ON H =1.0
FREE SURFACE GRADIENT COMPATIBILITY =0.9
TREATMENT OF THE LINEAR SYSTEM =2
SUPG OPTION =0;0;2;2
TYPE OF ADVECTION =1;5;1;1

Dear Consonaut,

you make me really worry.

It seems that we have a problem with the Sisyphe algorithm which cannot deal with dry nodes properly. Maybe it is important: In both cases, ours and yours, we have dry nodes which are not a part of a dried-out element. In our case, they are single nodes or a row of nodes in the middle of the domain, in yours, they are on the boundary. And it seems it does not matter, if we have suspended sediment transport or pure bed-load transport.

It seems that when these zero-velocity nodes are not a part of an element which would be marked by Telemac-2D as dried, they are a part of the area where the sediment transport takes place. And, due to some "feature" of the algorithm, the sediment gets in, but never gets out, so it deposites. And deposites. And deposites...

As mentioned in a post above, I suspect some upwinding blind to the real situation, but who knows. Of course there could be other explanations.

Have you solved your problem somehow?

Best regards,
Jacek

Dear jaj,

to clarify my post:
the picture is from a pure hydrodynamic calculation with Telemac-2D, where the bottom elevation at the boundary nodes is 0.5m higher then then expected water surface elevation.

Therefore, connecting this discussion and coming back to the thread - Sisyphe in transcritical flows - , in the first step the problem comes from Telemac-2D.

So, playing with thresholds in Sisyphe would be an option.

I didn't solve the problem.
I'm completely with you:
"And maybe it is better to improve an existing algorithm than to "expand the number of algorithms"... I prefer only one option working robustly than a choice of ten fancy things which do not at all deliver what the people really want."


Best regards,
Clemens

Dear Clemens and Jacek

If you don't mind I will invite myself to your discussion :)

The problem is, as well checked up by both of you, is linked to the wetting drying phenomenon. More precisely, it is due to the over-estimation of the velocity magnitude in very shallow regions. These velocities are used to compute the shear stress (and then the Shields number) which is the main parameter that triggers the sediment transport for both bedload and suspension.
That's why, in a disagreement with the Jacek conclusion, I think that this is not a problem of how the algorithm of wetting-drying performs (the recommanded algorithm for FEM, corrects "manually" the computation of the gradient of free surface in order to recover the lake at rest (the well-balaceness or C-property)). but, as I said earlier, it is a problem of unphysical velocities.

The overestimation of the velocity in drying area is due to the dominating effects of friction terms who are horribly non linear.
For Jacek case, the bedload capacity is overestimated between a wet node and its dry neighbor, that's why a continuous sediment flux bring material that settles at the dry node. For Clemens case, the same effects is observed since the equilibrium near-bed oncentrations formulas are based on the critical Shields parameter.

I think that this is still an open problem, and this is also the case with friction treatment. However, a short term solution is to have a special care on the wetting and drying zones in your models. This can be through a better mesh refinement of these areas (Jacek, you mesh looks a bit coarse ) and in the vicinity of high bathymetry gradients. Moreover, more diffused hydrodynamics can help seriously to stabilize the decoupled resolution algorithm of Telemac-Sisyphe (this brings us to the endless dilemma of the smart amount of diffusion), and then velocities are less overestimated.

Finally, as mentioned by Clemens, you can try to use the finite volume kernel which uses a different wetting and drying algorithms which seems to be more efficient and that gave nice results when used for sediment applications. This might cause a non negligeable increase of CPU time, it's up to you !

I hope that this can help to advance to a clearer view on these philosophical items.
With my best regards

Riadh

Dear Riadh,

thank you very much for your careful and plausible explanation of the possible cause of the problem. Fell free to disagree with me, I am happy with it and I would be also vary satisfied to discover that the problem is not connected with upwind/downwind treatment in a principle algorithm. (And anyway, I am definitely not a fan of morpho-logo-dynamics!)

Yes, some devilish manipulation of the friction terms in the shallow areas depending on the water depth could be also a very plausible explanation for some upwind/downwind debalance in the bed load transport in this case! Or, any other case which results in a huge and unnatural differences in stresses driving the bed load.

I have checked in the users parameter file for the friction, this is Nikuradse, in friction_calc.f, but it is rudimentarily secured against bad values (see AUX):

! NOTE: 11.036 IS 30.D0/EXP(1.D0)
DO I = N_START, N_END
AUX=MAX(1.001D0,HC%R(I)*11.036D0/CHESTR%R(I))
CF%R(I) = 2.D0 / (LOG(AUX)/KARMAN)**2
ENDDO

I know implementations in which below some (user-given) small depth the hugely steep Nikuradse is replaced by a linear function and it is secured that even by depth=0 there is some friction. How it works with the bed load physics... no idea. See Delft-3D manual for example, they simply assume the beginning of the logarothmic velocity profile is a few z_0's below the bottom.

Yes, the user (it is not me, I have a whole group hiding shyly behind my back) has been recommended to stabilise the hydrodynamics -- there is quite a lot of Telemac-2D Parkinson-like shivers of the free surface and an advection scheme for tracers (14) is applied for the momentum, etc. And, of course, a massive refinement and lowering of the time step is required.

According to the "finite volume kernel" I would say we deal here with a lowland rivers confluence, but agreed, with a very ambitious setting of the problem, and Froude numbers swiching from subcritical to well over one above the sand bank. Please note this is one of these 2D-models calibrated over very variable discharges/waterlevels with the "small roughness -- huge Elder dispersion technology"... Nothing good for hyperbolic Godunov/Riemann/etc.

What roughness coefficient / roughness treatment would you recommend for this case? My experience for shallow areas says Manning, but then we have a very, very different calibration... Any other hints for your stress-debalance-theory?

Best regards,
jaj

Dear Jacek

I'm also a fan of the Manning/strickler law. The problem is that I do not have a magic way that gives straitforwardly a well calibrated/well-optimized-shear-stress model.

In a technical point of view, there are several keywords that help tuning the computation of shear stress/Shields critical number. Just have a glance on the user manual otherwise I'll send you some slides that sum up all the stuff.

A more interesting idea, in my opinion, is to investigate the way to dynamically vary the friction term (or coefficients) as a function of sediment characteristics and the hydrodynamics at these very shallow areas such as we can reproduce more physical behaving of water and bed.

To be continued ...

with my best regards

Riadh

Dear Riadh,

I have followed the slides you have delivered to me to the people in question here, thank you, I hope it will help them in their pursuits.

I know people modelling with Nikuradse in the deeper areas and with Manning in the shallows, and they are very happy with their approach. However: This is done in a 2D-model with horizontal z-layers. But: It translates easily to 2D.

Yes it is very true that the roughness coefficient (or, speaking broader and more precise, bottom stress driving the bed load) for such indirectly coupled bed load transport and hydrodynamics should be the function of the material at the bottom and -- at least for 2D-hydrodynamics -- the very function of the water depth and "flow conditions". This is the very crucial coupling parameter!

Unfortunately, so far as I know, this approach is not applied here, people compute with constant (albeit zonally varying) roughness coefficients mostly obtained via calibrations without transport. It is then applied blindly and independently from the fact what material travels where and how the water levels / depths change in this process...

Additionally, please remember "Elder turbulence model" is widely applied (with its original huge multiplicators) for momentum dispersion, which is the linear function of the water depth and roughness velocity...

And then a sand/gravel bank in the middle of the river grows up to the water surface... It has the original roughness coefficient of the river bed and the momentum dispersion goes down to zero over its crest.

I hope that your "overestimated bed stress theory" works in our case, because if it is true, we should have not only dried-out nodes growing up and up to the sky, but also one-node scours getting deeper and deeper -- until some slope effects switches on. The bed form is really very nervous in this case, but due to the inadequate resolution it is not clear what is what -- especially with the infamous Telemac-2D Parkinson shivering of the free surface above it. So I cannot say over-scouring in single nodes appears there as well...

But anyway, please have a look at the positive depths algorithm in a patch where a node with velocity equal to zero is surrounded by nodes with some higher velocity vectors pointing roughly in the same direction and analyse the situation. OK, you suggest diminishing these vector values or making them less visible in the bed stress computation and effectively lowering the transport. But: Is the fluxes computation in such patch really OK?

Best regards,
jaj