Assessment of simplified momentum equations for free surface flows through rigid porous media

In many applications, free surface flow through rigid porous media has to be modeled. Examples refer to coastal engineering applications as well as geotechnical or biomedical applications. Albeit the frequent applications, slight inconsistencies in the formulation of the governing equations can be found in the literature. The main goal of this paper is to identify these differences and provide a quantitative assessment of different approaches. Following a review of the different formulations, simulation results obtained from three alternative formulations are compared with experimental and numerical data. Results obtained by 2D and 3D test cases indicate that the predictive differences returned by the different formulations remain small for most applications, in particular for small porous Reynolds number ReP < 5000. Thus it seems justified to select a simplified formulation that supports an efficient algorithm and coding structure in a computational fluid dynamics environment. An estimated accuracy depending on the porous Reynolds number or the mean grain diameter is given for the simplified formulation.


Introduction
Geotechnical and coastal engineering applications often require simulation procedures that can model partially and fully saturated granular media. Also in, e.g., biomedical engineering and manufacturing, the flow through porous materials has to be considered. A variety of approaches for unconfined seepage problems can be found in early literature. Applications of Baiocchi's method (Baiocchi et al., 1973;Bruch, 1991;Bardet and Tobita, 2002) and Signorini's method (e.g., Chen et al., 2011) solve the unconfined seepage problem from the viewpoint of saturated zones. For seismic analyses of earth dams and embankments, initial hydromechanical conditions for transient analyses can be introduced in a unique approach with the equations developed for soil dynamics (Zienkiewicz et al., 1999).
In this paper, the governing equations for the simulation of instationary, incompressible, immiscible two-phase flow through rigid porous media are explored. Studying the literature, it is found that different formulations exist for instationary flow, e.g., Van Gent (1992), Liu et al. (1999), Hsu et al. (2002), de Lemos (2006, del Jesus et al. (2012), Uzuoka and Borja (2012), Higuera et al. (2014), Jensen et al. (2014), Higuera (2015), Larese et al. (2015), and Losada et al. (2016). The differences primarily refer to the transport terms of the momentum equation and are particularly influential at the boundary of the porous material or for an inhomogeneous porosity. To assess the predictive discrepancies in formerly investigated 2D and 3D test cases, a unified formulation is sought in this work. Parts of the deviations addressed herein were already outlined by del Jesus et al. (2012), Jensen et al. (2014), Higuera (2015), and Losada et al. (2016). The related discussion of equations will be extended in regards to the Vol. 5, No. 2, 2023, 159-177 Experimental andComputational Multiphase Flow https://doi.org/10.1007/s42757-022-0133-y continuity equation and momentum equation in this work, and a direct comparison of the formulations of Jensen et al. (2014) and Higuera (2015) will be done in the simulation studies using a uniform framework. The approach of this paper is to firstly list variations of the governing equations from the literature. Subsequently, two options are implemented in a cell-centered Finite Volume solver using a volume of fluid approach to solve the two-phase flow. The two options are very close to the approaches described in Jensen et al. (2014) and Higuera et al. (2014) but were developed based on the work of del Jesus et al. (2012). Moreover, results of a third edge-based formulation by the authors (see Larese et al., 2015) are compared. The three formulations are applied to illustrative 2D and 3D cases, which feature different porous Reynolds numbers and results are discussed. Results reveal that certain simplifications can be introduced without a traceable predictive quality change.
The paper aims to contribute to the understanding of dominant influences in the equations for a two-phase flow through rigid porous medium and deliver estimated accuracies for simplified formulations. The paper is structured as follows. In Section 2, the governing equations are discussed in the context of the existing literature. Different notations are compared and two different equation sets to be tested are derived from these considerations. A third simplified set of equations is also introduced in this section. Section 3 and Appendix A outline the numerical treatment for a cell-centered, second-order accurate Finite Volume method. Applications of the different formulations are discussed in Section 4, which also ranks the formulation related influences against approximation and porous resistance force modelling aspects. Since Zienkiewicz et al. (1990), Van Gent (1992), and Higuera et al. (2014) and previous sections lead to the conclusion that the influence of different formulations of the nonlinear terms in the momentum equation on accuracy is small for lower porous Reynolds numbers compared to resistance model influences. This hypothesis is tested on a test case with varying Reynolds numbers in Section 4.4. Not only is the maximum formulation error dependency on the porous Reynolds number given but it is also related to the mean grain diameter of the porous material in order to be able to estimate the simplified model's accuracy without knowing the nature of the porous flow at hand. Final conclusions are drawn in Section 5.
The notation uses standard lower case Latin subscripts to mark Cartesian tensor coordinates and Einsteins summation is used over repeated lower case indices. Symbolic notations of vectors and tensors employ underlining (e.g., v ).

Mathematical model
The simulation of two-phase flow through rigid porous material is simulated solving the Navier-Stokes equations in a Eulerian framework, modified so to account for the presence of the porous media. The present paper studies only immiscible fluids featuring incompressible bulk fluids and a sharp interface. In such cases, Volume of Fluid (VoF) or level set (e.g., Hirt and Nicholls, 1981;Sussmann, 1994) interface capturing approaches, which assume a unique velocity field, are frequently used simpler alternatives. They involve one continuity and one momentum equation for the mixture, which are supplemented by an equation to identify the local fluid phase and simple equations of state for the mixture properties. Hence, VoF and level set results are the basis of displayed computed own results.
Following the literature, different formulations of the governing equations are discussed in this section. We start with the definition of a porosity n, defined as the ratio of the fluid volume F V to the total volume V: The total volume dV consists of the fluid volume d F V and the volume of the rigid porous material

Continuity equation
Considering the mass conservation of the fluid inside the fluid volume for an incompressible fluid: leads to the usual zero-divergence differential form for one-phase flows with constant density ρ : Introducing immiscible two-phase flow (i.e., air and water) through a rigid porous medium and using a VoF approach, the density is obtained from an air mixture fraction air, / with constant bulk densities A ρ and W ρ . Introducing Eq. (4) into Eq. (2), again the zero-divergence differential form Eq. (3) follows when the immiscible condition D / D 0 c t = is applied.
Integrating over the total volume V instead of the fluid volume F V nV = and introducing the Darcy velocity as is also stated in del Jesus et al. (2012). The second term obviously vanishes for spatially constant porosities. If the simulation domain includes a porous material with constant porosity and an area without porous material, the second term only influences the results at the boundary of the porous material. In Section 4.1, it will be shown that the term has a neglectable influence on results for an artificial test case without resistance forces using the discretisation of the second term given in Appendix A.1.2.
In Jensen et al. (2014), it is claimed that the continuity equation stated in del Jesus et al. (2012) is not applicable and a 1D example is given to support this statement. A 1D flow with constant volume flux and a decreasing porosity in flow direction is used to show that the fluid velocity gradient is not zero but the Darcy velocity gradient is. This is true but neglects that the divergence of the fluid velocity has to be zero in the fluid volume as stated in Eq. (3). For a 1D flow, Eq. (3) can be expressed as Therefore if n is decreasing, also F A is decreasing and with that an increasing fluid velocity follows from Eq. (3). Also it can be noted that the second term of Eq. (5) vanishes for 1D problems.
To study the influence of the momentum equation formulations, the continuity equation is applied as in, e.g., Hsu et al. (2002), Jensen et al. (2014), or Larese et al. (2015 where the second term is neglected and therefore the divergence of the Darcy velocity ˆi v vanishes

Momentum equation
As a starting point, the well-known momentum equation of a fluid inside a fluid volume is used to derive the momentum equation in porous media. For each fluid volume d F V : The porous momentum equations are stated with minor modifications in the existing literature. An overview of the respective terms of a generic momentum equation (Eq. (9)) is given in Tables 1 and 2 for a selection of publications. In order to compare the differences exposed by these formulations, all equations are written in a unified form. The selected generic momentum equation refers to Eq. (9) which consists of the inertia term of the Darcy velocity I i a in addition to a convective C i a and a diffusion  Uzuoka and Borja (2012) Hsu et al. (2002) i ng Liu et al. (1999) -Van Gent (1992 term .  Generally, two different versions of the convective term exist, one based on the divergence of the fluid velocity (de Lemos, 2006;del Jesus et al., 2012;Uzuoka and Borja, 2012;Higuera et al., 2014;Jensen et al., 2014;Losada et al., 2016) and one using the divergence of the Darcy velocity (Van Gent, 1992;Liu et al., 1999;Hsu et al., 2002;Larese et al., 2015). Both formulations only agree if the porosity is homogeneous. For a spatially non-constant porosity, the first version follows from the derivation of the porous momentum equations, see above or, e.g., del Jesus et al. (2012). The second version is derived by considering the fluid velocity transporting the Darcy momentum. This imposes a lack of momentum transport at the boundaries of the porous material as highlighted in del Jesus et al. (2012). The amount of error depends on the face interpolation of the porosity in this case. Higuera et al. (2014), Jensen et al. (2014), and Losada et al. (2016) consider a variable density in the convective term. In Jensen et al. (2014), the momentum equation is given in the conservative form, which can also be applied to compressible fluids.
The fluid friction within the porous material is neglected by Uzuoka and Borja (2012). This is related to the focus of their simulations on predicting the deformation of poroelastic solids rather than the fluid flow inside the solid. Uzuoka and Borja (2012) concentrate on materials with a low intrinsic permeability. Van Gent (1992) neglects the viscous term of the fluid flow and instead includes related influences in a nonlinear term of the porous resistance. Except for minor differences as regards the fluid properties, all other formulations besides del Jesus et al. (2012) and Higuera et al. (2014) Larese et al. (2015), and Losada et al. (2016). In del Jesus et al. (2012), it is stated that Hsu et al. (2002) use a different pressure term which gives a smaller pressure drop. This cannot be confirmed by the present results. The formulation of de Lemos (2006) again coincides with the other formulations for a spatially constant porosity.
The three main variants of the porous resistance force are used by the tabulated formulations: Darcy's formula (Uzuoka and Borja, 2012), a Forchheimer formula (de Lemos, 2006;Larese et al., 2015), and a Forchheimer formula with an added mass term (Van Gent, 1992;Liu et al., 1999;Hsu et al., 2002;del Jesus et al., 2012;Higuera et al., 2014;Jensen et al., 2014). The use of the formulation depends on the nature of the investigated flow field which can be described by the pore Reynolds number. A detailed description of resistance laws and their range of application is given in Van Gent (1992) and Losada et al. (2016). Further details will be provided in Section 2.3 and Appendix B.
In order to compare different versions of the equation system, two momentum equations featuring different sets of variables are discretized in this work. Using the unified form of the equations, i.e., ˆ1 1ˆî Apart from the usage of a different velocity formulation, the equations differ in the diffusion term where Eq. (10) includes the divergence of the porosity hidden inside the Darcy velocity. In the diffusion term in Eq. (11), the porosity is excluded and therefore the divergence of it is not considered. This difference is especially relevant for non-constant porosities as well as at the boundaries of porous materials. For temporally inconstant porosities also the time term of Eqs. (10) and (11) leads to a different discretisation. Additionally numerical results are compared against the results from Larese et al. (2015), which provides a third version of the governing equations. Equation (10)  Earlier models like Van Gent (1992), Liu et al. (1999), and Hsu et al. (2002) gave already very good results for free surface flows through porous media without complying to the latest VARANS formulation (e.g., found in Losada et al. (2016)). Losada et al. (2016) also emphasize the high influence of the porous resistance model on the results and this is confirmed in Section 4.2. This suggests that a correct implementation of the nonlinear terms in the momentum equation might be of minor importance for certain engineering problems. Results from a simplified approach of the momentum equation where only the porous terms are added to the momentum equation are compared against results from a momentum equation where the porosity is taken into account in all terms in Section 4.4. A method that relates these deviations purely to the porous material properties in order to estimate the correctness of the simplified method beforehand (to calculate the porous Reynolds number the velocity magnitude already has to be known) is suggested. Other simplified approaches are also mentioned in Higuera et al. (2014), it is thought that the level of deviation given in Section 4.4 for simplified methods can be transferred to other simplified momentum equation approaches as well and therefore give a guidance on choosing an appropriate model for engineering problems. It is important to note that the correct treatment of the free surface concentration transport equation has to be secured in order to use the accuracies mentioned above.

Porous forces
The analysis of the flow through a porous media is well established in an engineering context. The initial model relates to Darcy's linear law, i.e., ˆR which is known since 1856. Its range of applicability is discussed in length by Polubarinova-Kochina (1952), Vant Gent (1995, and Larese (2012). For small porous Reynolds numbers: where 50 D refers to the mean grain size of the porous material. A linear law is sufficient to represent the resistance forces acting on the fluid. Nonlinear behaviour starts at [1 10] P Re Î -, see Gu and Wang (1991), and turbulence appears around [60 150] P Re Î as stated in Larese (2012). In 1901, Darcy's law was extended by Forchheimer (1901) using a nonlinear term to represent the contribution of turbulence and parts of large-scale convective transport. An additional transient term was added by Polubarinova-Kochina (1952) which takes added mass phenomena into account. The latter leads to an extended Forchheimer equation, viz.
Expression (13) still represents a state-of-the-art approach for flows through porous material. The resistance force factors A  , B  , and C  are best obtained by experiments.
Many authors tried to derive general expressions which depend on the porosity n, the physical properties of the fluid, the mean grain diameter 50 D , and additional constant, cf. Forchheimer (1901), Kozeny (1927), Ergun (1952), Engelund (1953), and Bear (1972). A table referring to the models used by an excerpt of the previously mentioned authors is given in Appendix C. Further models can be found in, e.g., Losada et al. (2016).

Two-phase model
The majority of authors used an interface capturing approach (Hirt and Nichols, 1981;Sussman et al., 1994;Anderson et al., 1998) which identifies the position of the interface from a scalar indicator or concentration or phase field function, cf. also related reviews in Mirjalili et al. (2017) and Tryggvason et al. (2011). Larese et al. (2015) rely on the level set method to detect the phases and only solve the momentum and pressure equation for the water phase. The approach offers the advantage of an inherently sharp interface for immiscible flows.
In this work, the fluid is assumed to consist of two immiscible phases, i.e., air and water. Using a VoF approach, the fluid properties are obtained from an air mixture fraction where ,

A W ρ ρ and ,
A W μ μ refer to the bulk properties of the two fluids, which are deemed to be constant in the present study. The employed equation governing the mixture fraction c follows from the continuity equation (3) and the immiscibility condition, i.e., D D 0 c t / = : Mixture fraction values c = 0.5 are used to identify the interface.

Numerical method
The present paper refers to a pressure-based Finite Volume (FV) formulation using a cell-centered, co-located variable arrangement on unstructured polyhedral grids, cf. Ferziger and Peric (2008). Such frameworks are frequently employed in fluid engineering applications. Integrals are approximated using a second-order accurate mid-point integration rule. Implicit first-order time discretisations are used and convective fluxes are approximated by first-order upwind (UDS) baseline formulae and subjected to explicit deferred corrections for higher-order approximations. In the present study, a simple flux blending scheme featuring 70% second-order central differencing are used. Diffusion fluxes are obtained from central differences, which also employ a deferred correction approach to account for non-orthogonality and face interpolation related issues. The procedure largely follows Ferziger and Peric (2008). Further details are given in Yakubov et al. (2015) and Völkner et al. (2017). We firstly outline the governing system based upon ˆi v and .   9), which is in line with e.g., del Jesus et al. (2012), for a rigid inhomogenous porous material with exception of the diffusion term (see Eq. (10)). The difference in the diffusion term is attributed to our aim at obtaining a conservative form that serves a FV implementation. This conservative form follows from dividing Eq. (8) by n which is valid for rigid porous media, i.e., Since in this work the porosity is constant in time, this equation is completely in line with Jensen et al. (2014). Note that to let the second bracket on the left hand side vanish, the conservative air mixture fraction transport must employ the fluid velocity and not the Darcy velocity, cf. (15), viz., for a fluid volume d .
and a conservative air mixture fraction equation is used to close the system, viz., This formulation is in line with Higuera et al. (2014) and has the advantage of only adding two resistance terms to the momentum equation in a standard fluid solver. It uses a different diffusion term as in Jensen et al. (2014) or Eq. (16). The difficulty of this approach lies in the correct discretisation of the continuity equation, as is also stated in Higuera et al. (2014).

Simplified momentum equation for the ˆi v , p F formulation
The simplified formulation is the set of governing equations in Section 3.1 with a simplification of the aforementioned momentum equation (16). The reduced version of the momentum equation neglects the division of the left hand side terms of the momentum equation by the porosity: Similar methods are mentioned in Higuera et al. (2014), albeit a discussion on the exactness of the method is not given.
for demonstrative purposes. The pressure-correction scheme requires to derive a reduced momentum equation prior to modifying the continuity equation. The parts of the convective and diffusive terms which are expressed over the neighbouring Darcy velocities are usually neglected. Also the gravity term is usually neglected. This leads to the semi-discrete reduced momentum equation: One argument is, that for n ≈ 0.5 (which is a typical porosity) and a domination of the porous resistance forces, gravity forces, and pressure term, the simplification of the time term, convective term, and diffusive term has a low influence on the result of the momentum equation: Furthermore for lower Reynolds numbers, creeping flow situations, featuring small velocities, and small velocity derivatives can be expected. In the limit of small coefficients *, n v P A and * nb A , one immediately recognizes that the pressure-correction equation (30) reduces to a subset of an equilibrium of forces between the right hand side terms of Eq. (16). In this case, the truncated momentum equation (23) essentially converges to the "pure" Darcy flow equation. The previous paragraph can also be summarized as: due to small velocities and small velocity derivatives inside the porous medium, the left hand side terms of Eqs. (16) and (21) both approach zero and therefore the simplification is acceptable. The advantage of using Eq. (21) is the simplicity of the implementation into a standard fluid solver. No modification of the momentum equation is needed (the variable F v is changed to v ) except the addition of the porous forces terms as well as no modification of the pressure correction scheme is needed at all. A simple modification of the free surface equation as explained above (see Eq. (17)) has to be performed. Both alterations may be implemented in a user coding environment in commercial software. A further advantage of the simplified method is its higher robustness towards different treatments of the porosity value at faces since no face values of the porosity are needed in the discretisation of the simplified method.

Simulation studies-2D/3D test cases
In this section different 2D and 3D test cases, i.e., variants of unsteady dam break flows and a steady embankment flow, are investigated using different porous materials. The cell-centered porosities (n) and the constants of a reduced Engelund law ( A  , B  , 0 C º  ) are defined in line with data used by previous computational and experimental work. Face values of the porosity follow from a linear interpolation. To estimate nonlinear influences on the porous resistance, average and peak porous Reynolds number P Re are calculated from the present simulations. The unconstraint external flows feature a fairly low Reynolds number. Thus the fluid simulation model follows a laminar representation-i.e., no turbulent eddy viscosity, etc., is considered-and potential turbulent effects inside the porous material are deemed to be incorporated in the porous resistance force model. Results of the present simulations utilize a compressive HRIC approximation of convective fluxes (Muzaferija and Peric, 1998) in test cases B.1 and B.3 in Appendix B as well as the CICSAM scheme (Ubbink, 1997) in test case B.2 in Appendix B. An implicit first-order accurate scheme is used to approximate time derivatives.

Modelling influences
Modelling influences are outlined by Fig. 2, which compares experimental data for the free surface evolution published in Liu et al. (1999) with numerical results obtained from different computational models. Simulations refer to an edge-based level set solver reported by Larese et al. (2015) and the two VoF formulations as discussed in Section 3 on an isotropic homogeneous grid that features h/Δx = 42 which approximately refers to Δx = 0.332 cm.  Note that Larese et al. (2015) increased the experimentally reported initial column width from w = 28 cm to w = 29.8 cm. The increased width was also used by the present VoF simulations, to support the verification and assess formulation related differences. Moreover, we supplement FV results obtained from the ˆi v -based formulation for the experimental column width w = 28 cm. No visible differences can be detected for the results obtained from the two VoF approaches, i.e., using the ˆi v -based or F i v -based formulation. Comparing the results for the two different initial column widths predicted by the same FV formulation, it is seen that the correct initialisation provides a significantly better agreement with experimental data for all time instants except t = 4.0 s. A possible reason for this is described by Liu et al. (1999), as the free surface tends to adhere to the lateral glass wall. Therefore, the experimental data might overestimate the free surface elevation at the late time instant. Figure 2 also supports a comparison of different numerical methods using the same initial conditions and similar resolutions. Mind that Larese et al. (2015) did use a level set approach which is restricted to the water phase whereas the present approach resolves both fluid phases using VoF. Albeit the different governing equations, the results of the different solvers are in acceptable agreement for all time instants. Noticeable differences occur upstream the porous dam, and at the entrance of the water into the porous material. The transit of free surface elevation predicted by the FV solver is more continuous at t = 0.8 s, t = 1.6 s, and t = 4.0 s, whereas slightly more damming is observed in the results of Larese et al. (2015). Since the convection term of both formulations differs along the porous dam, we attribute these predictive differences to the different formulations. However, slight damming due to dynamic effects also appears for the FV study when changing to the experimentally reported initialisation, which reveals a small gap of air between the liquid dam and the porous medium. It is noteworthy that the differences between the FV VoF solver and the edge-based level set solver are far more considerable than the differences between the two VoF formulations described in Section 3.
A second test case of water flow through an inclined dam which consists of homogeneous rockfill material (see Appendix B.2) studies the difference of the formulations for a higher porous Reynolds number. Numerical results are compared with experimental data for the steady state in Fig. 3. Moreover, Fig. 4 provides a comparison of the numerical data for the transient buildup of the flow. Both aspects are discussed with regards to shape of the free surface, particularly in the porous regime. Experimental data is available at three locations inside the downstream half of the dam. A supplementary measurement position is located just downstream the dam. When attention is directed to the steady state depicted in Fig. 3, the free surface shapes returned by the investigated two FV formulations display hardly any difference and agree very well with measured data. Therefore, the formulation differences are deemed insignificant for this case. The edge-based solver predicts a free surface which is located slightly below the VoF solutions and-on average-slightly underestimates the experimentally reported elevation. The steady state findings are confirmed by the comparison of transient results. Figure 4 reveals that the free surfaces obtained from the FV and the edge-based solvers travel with a similar forward speed through the porous material, but a vertical/up lift tendency is observed by the VoF simulations when compared to the edge-based results. The latter might be attributed to the different formulation of the convective term and an increase of the porosity index in vertical direction. Also a marginal difference in the vertical velocity of the two VoF free surfaces can be observed in Fig. 4. This is thought to stem from the difference in the diffusion term of the two formulations, where the ˆ, i v F p formulation leads to a slightly higher diffusion term since the multiplication of the porosity and the fluid velocity is done inside the diffusion operator. It has to be noted that for a higher pore Reynolds number, the differences in formulations are increasing.
To demonstrate the influence of different formulations for a more dynamic flow, a 3D dam break case reported by del Jesus et al. (2012) is studied. The case is illustrated in Fig. 17 and refers to a three-dimensional extension of the crushed rocks case discussed in Section 4.3. Additional to the HRIC scheme, an explicit interface sharpening algorithm described in Manzke (2019) Larese et al. (2015). Physical features such as a sink behind the corner of the porous material at t = 0.4 s and the backwash reflected behind the dam at t = 1.33 s are represented in all three models. The shape of the advancing water front at t = 0.4 s and the free surface shape at t = 1.33 s differ due to differences in discretisation methods (FV-VoF/edgebased level set) and the applied model equations (with and without turbulence models). Due to similar discretisation methods, the present results generally agree with the results from del   Fig. 6 for two lateral planes at y = 10 cm and y = 50 cm. Both VoF methods discussed in Section 3 deliver very close results. The difference in the formulation of the diffusion term in Eqs. (10) and (11) seems to cause only minor deviations at the boundaries of the porous dam.
A reasonable overall agreement of the free surface development inside the porous material is found for all time steps and all displayed results. The main deviations between the present results and del Jesus et al. (2012) are an upstream shift in the free surface elevation of the present predictions at t = 0.5 s in the plane which does not cut through the porous dam. This indicates an increased damming experienced in the present predictions, which also explains why more water is found at t = 0.5 s in the results of del Jesus et al. (2012) in the y = 50 cm plane the dam. This indicates that a different mass transport through the porous boundaries is caused by the dissimilar treatment of the continuity equation. The same effect is thought to be responsible for the substantial differences of the observed free surface elevation at t = 1.38 s next to the dam.
The initial situation follows from Fig. 17. Since for all three cases above, the differences for the two VoF formulations are small, an artificial case is used to show the maximum difference obtainable. The neglect of porous resistance forces provides an illustrative example to emphasize the expected differences of the two formulations described in Section 3. A snapshot of the liquid body obtained from simulations without porous resistance, i.e., 0, 0, Fig. 7. The simulated liquid body is coloured by the Darcy velocity for the ˆi v -based formulation and fluid velocity for the F i v -based formulation. As expected by definition, the fluid velocity strongly accelerates in the porous medium without porous resistance due to the reduction of the wetted regime and decelerates when leaving the porous media. On the contrary, the Darcy velocity displays a more continuous evolution over the fluid domain. Therefore, the ˆi v -based formulation yields a smoother velocity field in the transition region, whereas the F i v -based approach returns large gradients that even yield slight oscillatory velocities in the transition region (see Fig. 8). This might relate to the assumption F A nA = in the continuity equation. Adding the neglected term of the full continuity equation (5), discretised as given in Appendix A.1.2 leads to no visible differences in results for this extreme case which underlines that the assumption made in Section 2.1 is valid. In the absence of the usually dominating porous resistances, the formulation differences related to the transport terms lead to a visible difference of the free surface elevation. For the simplified momentum equation, a high deviation in free surface shapes can be noted leading to the conclusion that the accuracy of the simplification depends on the relation between porous resistances and the nonlinear terms of the momentum equation.

Porous force model influences
Porous force model influences are displayed in Fig. 9 for FV simulations using the ˆi v -based VoF formulation in conjunction with the experimentally reported initial conditions described in Appendix B.1. The references to the porous models can be found in Appendix C. As indicated by the two Van Gent results, the influence of the added mass coefficient is negligible, which explains why it was hardly considered by other authors. In contrast to the formulation and modelling differences assessed above, changes to the resistance law have a large effect on the development of the free surface. This highlights the dominance of the resistance law to balance pressure and gravity forces inside the momentum equation. The resistance law fitting best to the experimental data is the one used by Larese et al. (2015) which was also used in the first test case of the modelling study above.

Porous Reynolds number influences
Porous Reynolds number influences were studied by a 2D dam break case which involves a higher porosity index and a different initial water column. Modifications were extracted from the crushed rocks case described in Liu et al. (1999). Results aim to assess predictive differences revealed by the two formulations discussed in Section 3 for a higher porous Reynolds number. ) was employed. The added mass coefficient C  was neglected. Results are again assessed by the temporal evolution of the free surface and compared with experimental data reported by Liu et al. (1999). Figure 10 reveals that the results returned by the two VoF formulations discussed in Section 3 are only marginally different and a fair predictive agreement can be observed. Particularly, the elevation dip observed downstream the porous layer from t = 1.0-1.6 s is represented well. Small deviations between the two simulations and between simulated and measured data can be observed inside the porous material for 0.8 s. t £ The numerical modelling influences are thus slightly more pronounced for the higher porous Reynolds number but still remain small.

Deviation of simplified momentum equation results
The deviation of simplified momentum equation results as a function of the porous Reynolds number provides an estimation of the errors due to the neglect of variable porosity influences inside the transport terms. Results obtained in the previous subsections suggest that the influence of transport terms is limited, whilst the pressure, volume, and drag forces on the right hand side of the momentum equation (16) balance each other towards an equilibrium. Therefore, the errors due to modifications of the transport terms might be limited. The goal is thus to obtain an approximate correlation between the related error and the porous Reynolds number, which can be used to justify a simpler approach  for a particular porous Reynolds number regime. To this extent, simulations were conducted with the Darcy velocity using a simplified momentum equation: For continuity reasons, the term inside the second square bracket of Eq. (25) vanishes and the remaining expression features the usual conservative fluid dynamic transport term on the left hand side together with the right hand side force terms that occur in the Darcy velocity based momentum equation (16). Equation (25) denotes a traditional conservative fluid dynamic formulation which is only augmented by a porous force term on the right hand side. The approach is therefore of interest, when either parts of the . Plotting the error measure l D / 0 l over the maximum porous Reynolds number in Fig. 11, a clear dependency of the deviation on the porous Reynolds number is observed. This indicates that, for small enough pore Reynolds numbers, the left hand side terms are obviously much smaller than the right hand side terms which in turn seem to equilibrate. Using a nonlinear least-squares (NLLS) Marquardt-Levenberg algorithm for acquiring a fit leads to = .
For porous Reynolds numbers lower than 10 4 , the maximum deviation of results obtained from the simplified momentum equation (25) falls below 5%. This motivates a possibility to simplify the momentum equation for all porous materials with mean grain diameters below 1 cm.

A priori estimation of deviation for simplified model
To estimate the deviation of the simplified momentum equation method from the correct implementation prior to the simulation, i.e., before having access to the porous Reynolds number, the computed deviations are plotted against the porosity for the different mean grain diameters D in Fig. 13. It is interesting that, for very small and very large mean grain diameters, the deviations seem to be almost   . constant over the porosities whereas for medium mean grain diameters, an increase of deviation with rising porosity is noticed. The latter is explained by the resulting higher porous Reynolds numbers for higher porosities. Clearly for some coastal engineering problems, mean diameters in the range of meters are possible and the simplified approach is not recommendable. For many other engineering problems dealing with sand or finer grains, a nonlinear simplified approach only yields very small deviations.

Conclusions
The aim of this work was to assess different formulations of the governing momentum equations for free-surface flows through rigid porous material. To this end, two exemplary formulations based on the Darcy and the fluid velocity were implemented in a VoF-based finite-volume procedure, and simulation results for 2D and 3D cases were compared against results from other literature reported formulations. It is found that for flows at porous Reynolds numbers up to 5000 P Re £ , formulation-based changes of momentum transport terms do only yield minor if not negligible deviations of the results, which are usually concealed by even subtle differences of the employed resistance law and the related parameters. Therefore, the particular choice of the employed momentum equations follows different motivations and further simplifications might be defensible from an engineering perspective. For a simplified formulation, the deviations from the solution can be estimated depending on grain sizes. In this regard it is found that for materials with mean grain sizes below 1 cm simplifications of the momentum equation lead to deviations in the free surface representation under 5%.
In conclusion, the optimal formulation might depend on the structure of the employed simulation procedure and the attainable efficiency and robustness. In this regard, reduced run-time and enhanced robustness point to a preference of the Darcy velocity ˆi v -based non-simplified formulation.

Appendix A Discretised equations
Here ( ) refers to porous mass flux across a face, the subscript f marks the cell faces and related variable values (see Eq. (14)), A denotes the face area and i A the outward pointing face vector. The scalar distance between the cell center P and adjacent neighbouring centers NB is labeled d. For each time step, the governing equations are iterated to convergence in a segregated manner. The superscripts n and m denote to the time step and the iteration index of the iterative procedure, respectively. Porous resistance terms are implemented in an implicit manner which improves the stability. The source term ˆi v S includes explicit terms which arise from different deferred correction contributions, e.g., higher-order convection, non-orthogonality, and interpolation corrections.

A.1.2 Discretised continuity equation using ˆi v and F p
A SIMPLE-type pressure correction scheme is used to fulfill the continuity equation, cf. Ferziger and Peric (2008) and Yakubov et al. (2015). The procedure follows a fairly standardised approach and is only briefly outlined here. It starts from the values ( 1m Note that the second term of Eq. (5) is neglected at this point. Subsequently Eq. (29) is used to determine ˆ. i v ¢ To avoid pressure oscillations occurring for a co-located variable arrangement, a special interpolation of the fluxes in ( ) is used (Rhie and Chow, 1983). It is important to note that the classical SIMPLE pressure correction scheme can be adopted without any change. Porous media influences are implicitly considered through the main diagonal coefficients ˆ, n v P A cf. Eq. (28).
To study the influence of the second term of Eq. (5) in Section 4.1, the full equation is discretised by adding two explicit terms in the equation including the scalar product of the gradient of the porosity and the Darcy velocity divided by the porosity. The velocity correction is discretised by linking it to the pressure correction from the previous iteration with Eq. (29): A.1.3 Discretised air mixture fraction equation using ˆi v and F p Using a first-order implicit time integration scheme and an implicit upwind-difference scheme part for the approximation of the convective term, the discrete mixture fraction transport reads , 1 , , where c S hosts the deferred correction terms of the compressive approximation. Regarding this, both the HRIC (e.g., Muzaferija and Peric, 1998) scheme and the CICSAM (e.g., Ubbink, 1997;Ubbink and Issa, 1999) scheme are used in the present study. With the air mixture fraction field, the values of the cell/face-centered density and viscosity are updated at each iteration by the cell/face-centered values of c using the simple linear equation of state (14).

A.2 Formulation using
As can be seen the time terms, the convection terms, the diffusion terms are reduced. The flux blending scheme is consistently simplified.

Appendix B Definition of test cases
B.1 Dam break flow through rigid porous material- Lin (1998) The first case refers to a dam break through porous material as proposed in Lin (1998) and Liu et al. (1999). Experimental data from Liu et al. (1999) is used to validate the results. The experiment features a water column with an initial height of h = 14.0 cm and an initial width of w = 28.0 cm, which passes through a matrix of glass beads with an average diameter of D 50 = 0.3 cm that forms a porous material with a porosity value of n = 0.39. The geometry is depicted by Fig. 15. Using the present results, a low average Reynolds number of Re P = 55 follows from the space/ time-mean Darcy velocity. The corresponding maximum velocity yields a Reynolds number of Re P = 125. These values indicate, that nonlinear or turbulent effects should be of relevance for the resistance modelling. Following Table 3, three different representations of the porous forces ˆR i f are found in the literature. Larese et al. (2015) use the Engelund equation, i.e., to determine the porous resistance from α = 1000 and β = 0.25, which leads to A  = 64.7 s −1 and Moreover, Liu et al. (1999) employ the Van Gent model with α = 200 and β = 1.1 according to Table 3. Using a Keulegan-Carpenter number, an estimated characteristic time of T = 1 s, as well as α = 200 and β = 1.1 yields A  = 67.3 s −1 , B  = 1470.5 m −1 , and an added mass parameter C  = 0.21.

B.2 Flow through an inclined porous dam
The second case refers to 2D experiments of water flow through an inclined dam which consists of homogeneous rockfill material as described in Larese et al. (2011). The selected case is illustrated in Fig. 16. It refers a constant volume flux of Q = 25.46 s -1 entering a channel through an orifice of 5 cm (height) × 246 cm (width), in which a dam of rockfill material is present. The porosity index of the dam reads n = 0.4052 and the average particle diameter refers to D 50 = 35.04 mm. This leads to a (computed) Reynolds number of Re P = 1335 using a space/time-averaged Darcy velocity and Re P = 4939 using the maximum Darcy velocity. The Ergun theory is used for the resistance law by means of ( ) In accordance with the experiments, a Dirichlet condition using a homogeneous velocity 1 F v = 0.207 m/s is used to simulate the inflow and zero-gradient conditions were used at the outlet. The size of the considered 2D domain is indicated by Fig. 16. The employed water density refers to 1000 kg/m 3 and the dynamic viscosity is assigned to μ = 0.001 Pa·s. Pressure data was experimentally recorded for the steady state flow at four different locations to detect the free surface and is used for comparison in this study.

B.3 3D dam break
A porous dam consisting of the crushed rocks material described in Section 4.3 is positioned as given in Fig. 17 in a domain which is 0.6 m in height, 0.6 m in width, and 1.2 m in length. The dam features a height of 0.6 m, a width of 0.3 m, and a length of 0.3 m. The initial water column is 0.6 m in width, 0.3 m in length, and 0.4 m in height. As in the 2D cases, a free surface height of 2.5 cm is initiated in the rest of the domain. For the space/time-averaged Darcy velocity, this yields a pore Reynolds number Re P = 1761. Using the maximum Darcy velocity observed in space and time, the pore Reynolds number reads Re P = 15099. To obtain a similar resolution as del Jesus et al. (2012), the homogeneous horizontal grid spacing was assigned to Δx = Δy = 1 cm and the vertical to Δz = 0.5 cm.