A robust Upwind Mixed Hybrid Finite Element method for transport in variably saturated porous media

Younes, Anis; Hoteit, Hussein; Helmig, Rainer; Fahs, Marwan

doi:https://doi.org/10.5194/hess-2022-153

Preprints

https://doi.org/10.5194/hess-2022-153

Preprints

27 Apr 2022

Status: a revised version of this preprint was accepted for the journal HESS.

A robust Upwind Mixed Hybrid Finite Element method for transport in variably saturated porous media

Anis Younes¹, Hussein Hoteit², Rainer Helmig³, and Marwan Fahs¹ Anis Younes et al.

Show author details

¹Institut Terre et Environnement de Strasbourg, Université de Strasbourg, CNRS, ENGEES, UMR 7063, 67084 Strasbourg, France
²Physical Science and Engineering Division, King Abdullah University of Science and Technology (KAUST), Thuwal, Saudi Arabia
³Institute for Modelling Hydraulic and Environmental Systems, University of Stuttgart, Pfaffenwaldring 61, 15 70569 Stuttgart, Germany

Received: 14 Apr 2022 – Discussion started: 27 Apr 2022

Abstract. The Mixed Finite Element (MFE) method is well adapted for the simulation of fluid flow in heterogeneous porous media. However, when employed for the transport equation, it can generate solutions with strong unphysical oscillations because of the hyperbolic nature of advection. In this work, a robust upwind MFE scheme is proposed to avoid such unphysical oscillations. The new scheme is a combination of the upwind edge/face centred Finite Volume (FV) method with the hybrid formulation of the MFE method. The scheme ensures continuity of both advective and dispersive fluxes between adjacent elements and allows to maintain the time derivative continuous, which permits employment of high order time integration methods via the Method of Lines (MOL).

Numerical simulations are performed in both saturated and unsaturated porous media to investigate the robustness of the new upwind-MFE scheme. Results show that, contrarily to the standard scheme, the upwind-MFE method generates stable solutions without under and overshoots. The simulation of contaminant transport into a variably saturated porous medium highlights the robustness of the proposed upwind scheme when combined with the MOL for solving nonlinear problems.

Anis Younes et al.

Status: final response (author comments only)

CC1:
'Comment on hess-2022-153', Thomas Graf, 25 May 2022

This manuscript introduces a new upwind Mixed Finite Element (MFE) method to solve for solute transport under saturated and unsaturated conditions. The new scheme is shown to avoid unphysical oscillations that is otherwise potentially caused by advection.
The manuscript is clearly written, and the mathematical context is sound and complete. The two illustrative results are useful and demonstrate very well the capabilities of the new method.
There are some minor comments that I ask the authors to address:
1. There are a number of filling words (especially "indeed") that need to be removed. It is also uncommong to use "the" in e.g. line 129 when generally addressing water content etc.
2. Description of boundary conditions of the two examples in section 4. is not complete. Both figures 4 and 7 should show all BCs for both flow and transport. As is, this is not the case and must be changed. Also, the text does not give the full description of all BCs for both flow and transport, which also has to be completed. A figure that shows the meshes for both examples is missing and would be very helpful. I understand it is an unstructured mesh which explains the fact that result in Fig. 5 is not symmetric. This should be mentioned. For simplicity, it would have been more efficient to use only the upper half of the domain due to symmetry reasons. The authors should give a reason why they did not do so.

Citation: https://doi.org/10.5194/hess-2022-153-CC1
- AC1:
  'Reply on CC1', Anis Younes, 06 Aug 2022
  CC1: 'Comment on hess-2022-153', Thomas Graf, 25 May 2022
  This manuscript introduces a new upwind Mixed Finite Element (MFE) method to solve for solute transport under saturated and unsaturated conditions. The new scheme is shown to avoid unphysical oscillations that is otherwise potentially caused by advection.
  The manuscript is clearly written, and the mathematical context is sound and complete. The two illustrative results are useful and demonstrate very well the capabilities of the new method.
  Answer: We highly appreciate the Reviewer’s positive appraisal of our work as well as careful reading of the paper. As detailed below, all comments are accounted for in the new revision.
  There are some minor comments that I ask the authors to address:
  There are a number of filling words (especially "indeed") that need to be removed. It is also uncommong to use "the" in e.g. line 129 when generally addressing water content etc.
  
  Answer: Corrected in the entire document.
  Description of boundary conditions of the two examples in section 4. is not complete. Both figures 4 and 7 should show all BCs for both flow and transport. As is, this is not the case and must be changed. Also, the text does not give the full description of all BCs for both flow and transport, which also has to be completed.
  
  Answer: Figures 4 and 7 have been changed to include all flow and transport boundary conditions. A full description of all BCs is also given in the text of the revised version.
  A figure that shows the meshes for both examples is missing and would be very helpful. I understand it is an unstructured mesh which explains the fact that result in Fig. 5 is not symmetric. This should be mentioned.
  Answer: Although we have used an unstructured mesh, the results of Figure 5 show almost symmetrical results, except the unphysical oscillations obtained with the standard MFE scheme.
  Note that in the revised version, the first test problem is also simulated using different mesh refinements to investigate time and space convergence of the new developed method.
  For simplicity, it would have been more efficient to use only the upper half of the domain due to symmetry reasons. The authors should give a reason why they did not do so.
  Answer: Yes, it is more efficient to simulate only the upper half of the domain because of symmetry, but we prefer to simulate the whole domain to have more representative results and show the capability of the new upwind MFE model to reproduce symmetrical results.
  
  Citation: https://doi.org/10.5194/hess-2022-153-AC1
RC1:
'Comment on hess-2022-153', Anonymous Referee #1, 07 Jul 2022

This paper deals with 2D numerical simulations for a coupled system arising in flow and transport in heterogeneous media. The mathematical model under consideration is the flow and transport in variably saturated porous media using Richard’s equation. A numerical scheme is developed for the discretization of this system by combining a mixed finite element method and a new upwind scheme for the convective term. 2D numerical results are presented to see the performance of the scheme for two tests for numerical simulation of contaminant transport into a variably saturated porous medium. The obtained results are satisfactory.

The subject is of interest and of current events. The authors made an interesting contribution for a difficult problem. The paper is well written and the results are of current interest. I deeply recommend the publication of this article.

The authors should clarify the following points:

- '{ }' notation should be defined to avoid confusion.

- The time discretization, the strategy used for the choice of the time step, the resolution of the nonlinear system and the linear systems should be specified.

- To ensure reproducibility of the results of the two tests presented, all necessary data including discretization and solvers etc. should be specified.

- It would also be interesting to give information about the environment in which the simulations were performed and the CPU times for each simulation.

- Can you comment on the extension of this approach to the 3D problem?

Citation: https://doi.org/10.5194/hess-2022-153-RC1
- AC2: 'Reply on RC1', Anis Younes, 06 Aug 2022
  
  RC1: 'Comment on hess-2022-153', Anonymous Referee #1, 07 Jul 2022 reply
  This paper deals with 2D numerical simulations for a coupled system arising in flow and transport in heterogeneous media. The mathematical model under consideration is the flow and transport in variably saturated porous media using Richard’s equation. A numerical scheme is developed for the discretization of this system by combining a mixed finite element method and a new upwind scheme for the convective term. 2D numerical results are presented to see the performance of the scheme for two tests for numerical simulation of contaminant transport into a variably saturated porous medium. The obtained results are satisfactory.
  The subject is of interest and of current events. The authors made an interesting contribution for a difficult problem. The paper is well written and the results are of current interest. I deeply recommend the publication of this article.
  Answer: We thank the Reviewer for his/her positive appraisal of our work. As detailed below, all comments are accounted for in the new revision.
  The authors should clarify the following points:
  - '{ }' notation should be defined to avoid confusion.
  Answer: The symbol designates the contribution from the adjacent element . This has been specified in the revised version.
  - The time discretization, the strategy used for the choice of the time step, the resolution of the nonlinear system and the linear systems should be specified.
  Answer: Done in the revised manuscript.
  - To ensure reproducibility of the results of the two tests presented, all necessary data including discretization and solvers etc. should be specified.
  Answer: We agree and add all data in the revised version.
  - It would also be interesting to give information about the environment in which the simulations were performed and the CPU times for each simulation.
  Answer: We agree and add information about the environment for the simulations and CPU times.
  - Can you comment on the extension of this approach to the 3D problem?
  Answer: The scheme can be seen as a combination of the upwind edge/face centred Finite Volume (FV) method with the lumped formulation of the hybrid MFE method and, as such, it can be easily extended to 3D problems.
  
  Citation: https://doi.org/10.5194/hess-2022-153-AC2
RC2:
'Comment on hess-2022-153', Anonymous Referee #2, 11 Jul 2022
The paper under review presents a new method to obtain an upwind discretization of the advection--diffusion equation, using mixed hybrid finite element. This is a non--trivial problem, as one would like to obtain a stable ans monotone method, and because the requirement for a hybrid discretization means that the balance equation needs to be local to each element. This is difficult to reconcile with the need for upwinding so as to obtain a stable method. The author review an existing method, and propose a new variation. The proposed method is validated on two numerical examples.

In this reviewer's opinion, the main advantage of the proposed method would that it avoids the need to discretize in time to derive the hybrid discretization. As argued by the authors, this enables the use of a higher order, or variable time step, discretization.

However, the paper as it is currently written suffers from several deficiencies, some major, some minor, that I would like to see addressed before publication.

Specific comments

I do not understand the "derivation" of the method that is given is Section 3 (starting on page 12). Why is it allowed, or useful, or even correct, to use a steady state form of the basic mass conservation equation (12) ? This approximation is used to obtain equation (24), and then again in equation (29).

As a consequence of the previous point, why would the proposed method be consistent ? And why do the authors expect that the method is stable (or satisfies a maximum principle) ? I understand that the paper, or even the journal, is not a numerical analysis paper, but any discussion shedding some heuristic light in the two issues of consistency and stability would be welcome.

If no analysis is possible, then the authors should at least include a numerical convergence study, say for example 1, to quantify the accuracy of the proposed method.

For the two examples discussed in Section 4, the proposed method is compared to the basic hybrid MFE method (that is equation (18) of the paper). But this method is already known to be non-monotone. A more significant comparison would be with the previous upwind variant from the Radu et al. (2011) reference quoted (that is the method described on page 11). This would allow the authors to discuss the possible pros and cons of their proposed method, in terms of accuracy, robustness and versatility. Even though I do not ask lightly that the authors do more numerical experiments, I feel this is really needed before the paper can be published.

The paper under review has a significant overlap with the recently published paper (see below) by the same authors. The authors should clarify the relationship between both papers.

A robust fully mixed finite element model for flow and transport in unsaturated fractured porous media,

Anis Younes, Hussein Hoteit, Rainer Helmig, Marwan Fahs, Advances in Water Resources, Volume 166, 2022, https://doi.org/10.1016/j.advwatres.2022.104259.

See Technical corrections in the attached file.
Citation: https://doi.org/10.5194/hess-2022-153-RC2
- AC3:
  'Reply on RC2', Anis Younes, 06 Aug 2022
  RC2: 'Comment on hess-2022-153', Anonymous Referee #2, 11 Jul 2022 reply
  The paper under review presents a new method to obtain an upwind discretization of the advection--diffusion equation, using mixed hybrid finite element. This is a non-trivial problem, as one would like to obtain a stable and monotone method, and because the requirement for a hybrid discretization means that the balance equation needs to be local to each element. This is difficult to reconcile with the need for upwinding so as to obtain a stable method. The author review an existing method, and propose a new variation. The proposed method is validated on two numerical examples.
  In this reviewer's opinion, the main advantage of the proposed method would that it avoids the need to discretize in time to derive the hybrid discretization. As argued by the authors, this enables the use of a higher order, or variable time step, discretization.
  However, the paper as it is currently written suffers from several deficiencies, some major, some minor, that I would like to see addressed before publication.
  Answer: We thank the Reviewer for the constructive comments which helped us to improve the quality of our manuscript. As detailed below, all suggestions are accounted for in the revised version.
  Specific comments
  I do not understand the "derivation" of the method that is given is Section 3 (starting on page 12). Why is it allowed, or useful, or even correct, to use a steady state form of the basic mass conservation equation (12) ? This approximation is used to obtain equation (24), and then again in equation (29).
  
  Answer: We agree and improve the presentation of the method in the revised version. More details are provided to explain how and why it is useful to use the steady state form of the mass conservation equation.
  The main idea of the new method is to combine the upwind edge centered finite volume method with the lumped hybrid MFE method. A new section is added to better explain the principles of the new method. Indeed, the standard hybrid MFE method is based on two stages: (i) the discretization of the transient mass conservation equation over the element E and (ii) the continuity of fluxes across the edge i sharing the elements E and E’. With the lumped formulation, the previous two stages are transformed to new ones as following: (i) the transient mass conservation over E is transformed to a steady state one and (ii) the continuity equation, which is seen as a steady state mass conservation equation at the edge level, is transformed to a transient equation. The new scheme keeps the time derivative continuous which allows the use of efficient high order temporal discretization methods via the MOL.
  These points are better specified in the revised version
  As a consequence of the previous point, why would the proposed method be consistent ? And why do the authors expect that the method is stable (or satisfies a maximum principle) ? I understand that the paper, or even the journal, is not a numerical analysis paper, but any discussion shedding some heuristic light in the two issues of consistency and stability would be welcome.
  
  Answer: In the new method, the term of mass (with time derivative) has a contribution only on the diagonal terms of the final system matrix. This improves the monotonous character of the hybrid-MFE solution. For dispersive transport, in the case of an acute triangulation, the maximum principle is respected whatever the heterogeneity of the porous medium (Younes et al., 2006). In the case of transport by only advection, the system matrix with the new upwind hybrid MFE scheme is always an M-matrix (a non singular matrix with ). The M-matrix property insures the stability of the scheme since it guaranties the respect of the discrete maximum principle i.e. local maxima or minima will not appear in the solution in a domain without local sources or sinks.
  Note that in the case of 2D triangular elements, the lumped formulation is algebraically equivalent to the nonconforming Crouzeix-Raviart (Crouzeix and Raviart, 1973) finite element method (see Younes et al, 2008). The nonconforming Crouzeix-Raviart method uses the chapeau functions as basis functions to approximate the concentration, like the standard finite element method, but seed nodes are the midpoints of the edges.
  These points will be specified in the revised version
  
  If no analysis is possible, then the authors should at least include a numerical convergence study, say for example 1, to quantify the accuracy of the proposed method.
  
  Answer: We agree and will present results of simulations for the first test case (with the analytical solution) using different mesh refinements in order to investigate the order of convergence in space and time of the new upwind-hybrid MFE method.
  
  For the two examples discussed in Section 4, the proposed method is compared to the basic hybrid MFE method (that is equation (18) of the paper). But this method is already known to be non-monotone. A more significant comparison would be with the previous upwind variant from the Radu et al. (2011) reference quoted (that is the method described on page 11). This would allow the authors to discuss the possible pros and cons of their proposed method, in terms of accuracy, robustness and versatility. Even though I do not ask lightly that the authors do more numerical experiments, I feel this is really needed before the paper can be published.
  
  Answer: In the revised version, we add numerical experiments which show that the new scheme has optimal first order convergence in time and space as the upwind scheme of Radu et al. (2011) and the non-hybrid upwind mixed method of Dawson (1978).
  Note that the upwind scheme of Radu et al. (2011) has good convergence properties but uses the approximation that the concentration at the edge i is the average of the mean concentrations of the two adjacent elements sharing i. This can be a rough approximation since with the hybrid MFE method, the trace of concentration can be significantly different from the average of mean concentrations of adjacent cells, especially in the case of a heterogeneous domain where dispersion can vary with several orders of magnitude from element to element. Further, with this scheme one cannot ensure continuity of both advective and dispersive fluxes at the interface since the continuity is prescribed only for the total flux. These two drawbacks are avoided with the new developed scheme.
  These points will be specified in the revised version
  
  The paper under review has a significant overlap with the recently published paper (see below) by the same authors. The authors should clarify the relationship between both papers.
  
  Answer: The aim of the present paper is (i) to develop a new upwind hybrid MFE scheme for advection-dispersion transport, (ii) to detail the main principles of the new scheme and show its benefits as compared to other upwind MFE schemes, (iii) to validate the method by comparison against an analytical solution, (iv) to investigate its order of convergence in time and space (in the revised version) and (v) to show its robustness when combined with the MOL for transport in unsaturated zones.
  The developed upwind scheme for transport has been recently successfully combined with the MFE method for fluid flow to simulate coupled flow and transport in unsaturated fractured porous media (Younes et al., 2022). The main objective of the paper of Younes et al. (2022) was to investigate the efficiency and robustness of a 1D-2D discrete fracture matrix (DFM) approach to model nonlinear flow and transport problems in fractured porous media with matrix-fracture and fracture-fracture fluid and mass exchanges.
  These points will be specified in the revised version
  A robust fully mixed finite element model for flow and transport in unsaturated fractured porous media, Anis Younes, Hussein Hoteit, Rainer Helmig, Marwan Fahs, Advances in Water Resources, Volume 166, 2022, https://doi.org/10.1016/j.advwatres.2022.104259.
  
  Technical corrections
  The used first in equation (27) and several times after that is very confusing. Does it mean “additional terms that the authors do not want to detail”, or rather “the same as before, but for element E0 instead of E"? I know the answer is the latter, but I wasn't sure the first time. This really should be clarified, preferably by writing the equation in full, at least by giving an explanation.
  
  Answer: We agree and change the text accordingly
  The exposition of the model in Section 2 should follow a more logical order: I would suggest starting with equation (3) (with (4) and (5)), then give equation (1) to say where comes from (and explain what is, then conclude with equation (6). The discretized equations such as (2) and (7) should only appear later.
  
  Answer: We agree and change the order of equations in the revised version
  I think there is a confusion between and on line 162. And the correct notation
  
  should be
  Answer: Corrected in the revised manuscript.
  In principle, the flux terms in equations (16) and later should have a superscript n + 1 attached.
  
  Answer: We agree and specify in the text that the subscript n+1 is omitted to alleviate the notation
  The physical units in Table 1 should be written in an upright font (as in the text before).
  
  Answer: Corrected in the revised manuscript.
  Example 1 has been used as a test case in the same context by (at least) two papers: [2]
  
  and [3], that should be cited.
  Answer: The two papers are referenced in the revised version.
  
  Citation: https://doi.org/10.5194/hess-2022-153-AC3

Anis Younes et al.

Viewed

Total article views: 422 (including HTML, PDF, and XML)

HTML	PDF	XML	Total	BibTeX	EndNote
301	104	17	422	3	2

HTML: 301
PDF: 104
XML: 17
Total: 422
BibTeX: 3
EndNote: 2

Views and downloads (calculated since 27 Apr 2022)

Month	HTML	PDF	XML	Total
Apr 2022	84	12	2	98
May 2022	91	28	5	124
Jun 2022	21	11	0	32
Jul 2022	53	17	4	74
Aug 2022	41	28	6	75
Sep 2022	11	8	0	19

Cumulative views and downloads (calculated since 27 Apr 2022)

Month	HTML	PDF	XML	Total
Apr 2022	84	12	2	98
May 2022	91	28	5	124
Jun 2022	21	11	0	32
Jul 2022	53	17	4	74
Aug 2022	41	28	6	75
Sep 2022	11	8	0	19

Viewed (geographical distribution)

Total article views: 394 (including HTML, PDF, and XML) Thereof 394 with geography defined and 0 with unknown origin.

Country	#	Views	%

Latest update: 28 Sep 2022

Short summary

Despite its advantages for the simulation of flow in heterogeneous and fractured porous media, the mixed hybrid finite element method has been rarely used for transport as it suffers from strong unphysical oscillations. We develop here a new upwind scheme for the mixed hybrid finite element that can avoid oscillations. Numerical examples confirm the robustness of this new scheme for the simulation of contaminant transport in both saturated and unsaturated conditions.


Total:	0
HTML:	0
PDF:	0
XML:	0