A Study on Darcy versus Forchheimer Models for Flow through Heterogeneous Landfills including Macropores

Abstract

Flow through heterogeneous landfills that include macropores may occur under Reynolds numbers higher than those where Darcy’s law is valid. Extensions, such as a Forchheimer approach, may be required to include inertial effects. Our aim is developing predictive models for such landfills that are built from the low-level radioactive waste and debris of dismantled nuclear power plants. It consists of different materials, which after crushing result in a spatially heterogeneous distribution of porous-media properties in the landfills. Rain events or leakage, for example, may wash out radionuclides and transport them with the water flow. We investigate here the water flow and consider an inclusion of macropores. To deal with possibly high velocities, we choose the Forchheimer model and, taking different Forchheimer coefficients into account, compare it to the Darcy model. The focal points of the study are (i) the influence of the macropores on the flow field and (ii) the impact of the choice of the Forchheimer coefficient both on the solution and the computational effort. The results show that dependent on their size, macropores can dominate the flow field. Furthermore, Forchheimer coefficients introducing more inertial effects are associated with considerably higher runtimes.

1. Introduction

When developing predictive models for flow through heterogeneous landfills, flow velocities beyond the validity of Darcy’s law [1,2] may be encountered. To incorporate inertial effects at higher Reynolds-number regimes, there are approaches related to Forchheimer’s extension of Darcy’s law discussed in the literature [3,4,5,6]. Our aim for the present study is to quantify the influence of employing Forchheimer’s extension on the results as well as on the numerical behavior for a particular set of model scenarios.

More precisely, we focus on flow events in heterogeneous landfills from debris of dismantled nuclear power plants after their clearance. Such debris consists mainly of demolished material such as concrete, plaster, bricks, metal, valves, piping, etc. According to [7], clearance is defined as “the removal of radioactive materials or radioactive objects within authorized practices from any further regulatory control by the regulating body”. Due to its nature, this waste can either be reused or disposed in municipal landfills [8,9,10]. The heterogeneity of the materials translates into a porous matrix of irregular characteristics. Matrices such as those are often called “structured” as they can contain large, continuous, structural pores, which can be called macropores [11]. According to [12], voids such as macropores are readily visible, and it is known that they can be continuous for several meters in both vertical and lateral directions. Most soils contain macropores of some sort; in [12], they are divided in categories depending on how they were formed. Since the cause of their formation is not always the same, their sizes can also vary widely [12,13,14,15]. A convenient approximation would be that pore diameters larger than 0.3–0.5 mm are typically classified as macropores. It should be mentioned that size is not alone a sufficient criterion to define a macropore; structure is important as well [12].

For our study, we follow the previous works from [16,17,18] on the modeling of the water pathway. They focused on 1D models under the assumption of homogeneous spatial parameters such as porosity and permeability. Additionally, they only considered regular rain events. We extend the model scenarios to 2D, incorporate spatial heterogeneity, and investigate heavy rain events, which are expected to happen more frequently in the future [17]. With the higher inflow mass from the heavy rain events and areas of high porosity (i.e., the macropores) in the heterogeneous landfills, the flow velocities inside the domain increase. Therefore, it is necessary to study this topic with an extension of Darcy’s law. We will investigate the differences between Darcy’s law and the Forchheimer extension both in terms of the results as well as of the computational performance, with a particular focus on quantifying the influence of macropores on the water flow field.

This paper is structured as follows: We introduce and explain the governing equations that were used for our model in Section 2. We continue with the description of the model scenarios in Section 3 and present results and an analysis in Section 4. We interprete and discuss our results in Section 5 and end with our conclusions in Section 6.

2. Governing Equations

In the landfills that we have in mind in this study, macropores are a relatively small proportion of the soil volume. However, they play an important role in the flow and transport of water and solutes [14], due to their irregular geometry and characteristics. Their effects on flow and how they can be described in state-of-the-art models have intrigued the interest of many researchers. According to [13], macropore flow can be categorized as a flow phenomenon that is initiated from the soil (matrix) surface and terminated at the deeper profile or groundwater, bypassing the intermediate soil profile. Macropore flow is also denoted in the literature as non-equilibrium and preferential [12,13,15,19,20].

Flow in domains containing macropores is often described with the Richards equation.

$$\frac{\partial \left(\varphi {\rho }_w{S}_w\right)}{\partial t}+\nabla ·\left({\rho }_w{\mathbf{v}}_w\right)-{\rho }_w{q}_w=0,$$

(1)

which only describes the flow of the water phase [14,20,21,22]. The limitations of the Richards equation with respect to the full interaction of water and gas phase can be overcome by using a two-phase flow approach [23]. Then, the mass conservation equations can be written as

$$\frac{\partial \left(\varphi {\rho }_{\alpha }{S}_{\alpha }\right)}{\partial t}+\nabla ·\left({\rho }_{\alpha }{\mathbf{v}}_{\alpha }\right)-{\rho }_{\alpha }{q}_{\alpha }=0;\alpha \in \{w,n\},$$

(2)

with $\varphi$ denoting the porosity, ${\rho }_{\alpha }$ denoting the density of phase $\alpha$, $S_{\alpha }$ denoting the phase saturation, ${\mathbf{v}}_{\alpha }$ denoting the Darcy velocity of phase $\alpha$, and $q_{\alpha }$ denoting the phase source term.

The phase velocity from Equation (2) can be calculated using the multi-phase extension of Darcy’s law as

$${\mathbf{v}}_{\alpha }=-\frac{k_{r\alpha }}{{\mu }_{\alpha }}\mathbf{K}(\nabla p_{\alpha }-{\rho }_{\alpha }\mathbf{g})\left[\mathrm{m}{\mathrm{s}}^{-1}\right];\alpha \in \{w,n\},$$

(3)

with $k_{r\alpha }$ indicating the phase relative permeability, ${\mu }_{\alpha }$ indicating the phase mobility, and $\mathbf{K}$ indicating the intrinsic permeability. We remark that we assume a rigid porous medium, namely, that the porosity and permeability distribution is fixed over time.

The preferential flow introduced in the domain because of the macropores (acceleration and deceleration in the pore throats) might also include inertial effects, which are characterized by a Reynolds number larger than unity ($Re>1$). For such cases, the multi-phase extension for Forchheimer’s law

$$-(\nabla p_{\alpha }-{\rho }_{\alpha }\mathbf{g})=\frac{{\mu }_{\alpha }}{\mathbf{K}{k}_{r\alpha }}{\mathbf{v}}_{\alpha }+{\beta }_{\alpha }{\rho }_{\alpha }{\mathbf{v}}_{\alpha }\left|{\mathbf{v}}_{\alpha }\right|$$

(4)

is an approach that has been used in the literature [3,23,24,25,26,27] to describe higher Re-number flow regimes in porous media, fractured, or macroporous media and also for flow near well regions [5,6]. It is a nonlinear relation with second-order corrections and, thus, can be viewed as an extension of the linear Darcy law for higher velocity regimes [26]. The nonlinear term of the equation accounts for inertial effects and contains a coefficient $\beta$ [${\mathrm{m}}^{-1}$], which is described as the effective non-Darcy-flow or Forchheimer coefficient. This coefficient is an empirical value specific to the porous medium [28,29] and has been a source of many controversies and ambiguities [30].

The inertial coefficient varies with geometrical characteristics of the pore structure. Methods for computing the $\beta$ coefficient can be roughly split into two categories: (i) the first one is concerned with fitting measurement data directly to the assumed model, while (ii) another approach assumes a direct correlation of $\beta$ to such fundamental properties of the porous medium as permeability ($\mathbf{K}$), porosity ($\varphi$), and in some cases tortuosity ($\tau$) [30].

The latter category can be further divided into theoretical and empirical correlations [24]. The empirical correlations of the inertia coefficient were initially determined experimentally for single-phase flow [25]. In order to account also for multi-phase flow, they were extended in terms of saturation and effective permeability correction. This way, the corresponding change in the Forchheimer flow coefficient due to saturation difference can be described [6]. The $\beta$ parameter values for typically studied porous media are reported to be greater than 5000 ${\mathrm{m}}^{-1}$ [26].

In this study, we investigate flow in a heterogeneous macroporous domain. The problem is formulated for building rubble and debris of a nuclear power plant on top of a landfill. Our aim is to assess if the flow in such a porous-media domain can be studied using the standard multi-phase Darcy approach or if an inertia correction should be added. The flow models are implemented in DuMu${}^{\mathrm{x}}$ [31]. DuMu${}^{\mathrm{x}}$ is an open-source simulator and research code in C++; its main intention is to provide a sustainable and consistent framework for the implementation and application of porous-media model concepts and constitutive relations. For more information, we refer to the official website dumux.org (accessed on 28 January 2022). The flow velocities are calculated by the extended multi-phase Darcy model, Equation (3), or Forchheimer’s model, Equation (4), respectively. For Forchheimer’s law, the Forchheimer coefficient is calculated with two approaches; thus, effectively, we consider two different Forchheimer models and evaluate their differences.

1.

A generalized equation for all media [32], which also is viable for multi-phase flow, formulated here according to [33]:

$${\beta }_{\alpha }=\frac{c_F}{\sqrt{\mathbf{K}}{k}_{r\alpha }}\left[{\mathrm{m}}^{-1}\right],$$

(5)

where $c_F$ is the so-called Forchheimer constant. A typical value of the constant for a range of investigated porous media is 0.55 [32].

If we insert Equation (5) in the general form of the Forchheimer Equation (4), we obtain the following equation for the first case:

$$\to {\mathbf{v}}_{\alpha }+c_F\sqrt{\mathbf{K}}\frac{{\rho }_{\alpha }}{{\mu }_{\alpha }}{\mathbf{v}}_{\alpha }\left|{\mathbf{v}}_{\alpha }\right|+\frac{\mathbf{K}{k}_{r\alpha }}{{\mu }_{\alpha }}(\nabla p_{\alpha }-{\rho }_{\alpha }\mathbf{g})=0.$$

(6)

2.

An equation that accounts for multi-phase flow more explicitly [6]:

$${\beta }_{\alpha }=\frac{c_{\beta }\tau }{\left(\mathbf{K}{k}_{r\alpha }\right)\left(\varphi S_{\alpha }\right)}\left[{\mathrm{m}}^{-1}\right];\alpha \in \{w,n\},$$

(7)

where $c_{\beta }$ is a numerical constant which is found to be approximately equal to $1.5\times {10}^{-4}\mathrm{m}$ for a specific range of investigated porous media [6].

By inserting this coefficient of Equation (7) in the Forchheimer Law, we can write the velocity of the fluid phase $\alpha$ for the second case:

$$\to {\mathbf{v}}_{\alpha }+\frac{c_{\beta }\tau }{\left(\varphi S_{\alpha }\right)}\frac{{\rho }_{\alpha }}{{\mu }_{\alpha }}{\mathbf{v}}_{\alpha }\left|{\mathbf{v}}_{\alpha }\right|+\frac{\mathbf{K}{k}_{r\alpha }}{{\mu }_{\alpha }}(\nabla p_{\alpha }-{\rho }_{\alpha }\mathbf{g})=0.$$

(8)

For each of the different models (based on Equations (3), (6), and (8), respectively), different macropore scenarios are simulated, and then, the results are compared with respect to their accuracy and also the efficiency of the simulations. Furthermore, we study the relationship between the parameters of each model, and the inertia flow is identified using the dimensionless flow criteria, Reynolds (Re) and Forchheimer (Fo) numbers (see Section 4.1).

The Forchheimer number equals the ratio of pressure drop caused by liquid–solid interactions to that by viscous resistance [34] and indicates how microscopic effects in the porous-media structure lead to significant macroscopic nonlinear effects [35]. Due to the consistent definition and physical meaning of the variables involved, it might be a better approach for the characterization of the flow regime [6]. In porous media, there is no clear threshold number that defines the transition between the flow regimes, and the nonlinearity in non-Darcy flow is introduced from inertia instead of turbulent effects. Therefore, in porous media, non-Darcy flow can occur for small Reynolds numbers [6]. Roughly, creeping Darcy flow is observed for $Re<1$ [36]. Non-Darcy inertial flow (or Forchheimer flow) occurs for a range of $1<Re<10$ [36]. Furthermore, the critical Forchheimer number is between 0.005 and 0.2, as stated in [6].

3. Model Scenarios

A major challenge in modeling landfills and their different materials as described above is the variation in their properties. Since it is impossible to predict the exact spatial distribution of the materials, we choose to create a spatial distribution with the R-project package gstat [37,38]. It requires specification of a variogram model; therefore, we chose the Matérn covariance function

$$F\left(h\right)=\frac{1}{2^{\nu -1}\mathsf{\Gamma }\left(\nu \right)}{\left(\frac{h}{r}\right)}^{\nu }{K}_{\nu }\left(\frac{h}{r}\right)$$

(9)

as available in gstat. According to [39], it works well for soil variograms. Its smoothness parameter $\nu$, appearing in Equation (9), allows for creating areas with less variation, which we assume to be the case in the landfill. For more information about the Matérn covariance function, we refer to [40]. The parameters used are shown in Appendix A in

Table A2. Parameters for the Matern model in gstat.

Parameter	Value
$\nu$	2.5
r	2
sill	0.02
anisotropy angle	90
anisotropy factor	0.4

.

The most relevant spatial parameters for us are the porosity and the permeability. Since we are interested in macropores, we decided to create a porosity distribution with gstat and then relate the permeability with the Kozeny–Carman relationship [41]

$$\mathbf{K}={\mathbf{K}}_{ref}\frac{{\varphi }^3{(1-{\varphi }_{ref})}^2}{{\varphi }_{ref}^3{(1-\varphi )}^2}.$$

(10)

Here, we use a reference porosity ${\varphi }_{ref}$ and a reference permeability ${\mathbf{K}}_{ref}$, which are calculated based on the grain-size distribution WT1 from the nuclear power plant Würgassen, as published in [42] and shown in

Table 1. Calculation of weighted average (average particle diameter) based on the grain-size distribution.

	Fraction [mm]	Average Fraction ${\mathit{d}}_{\mathit{i}}$ [mm]	Percentage ${\mathit{p}}_{\mathit{i}}$ [%]	Weighted Average ${\mathit{p}}_{\mathit{i}}{\mathit{d}}_{\mathit{i}}$ [mm]
1	<2	1	3.8	0.038
2	2–4	3	3.7	0.111
3	4–8	6	9.6	0.576
4	8–16	12	15.8	1.896
5	16–32	24	41.3	9.912
6	32–64	48	25.8	12.384
l				≈25 mm

. There, we also show the result for the mean particle diameter (l), which will be used in Section 4.1. Since the grain-size distribution is known, we are calculating the average particle diameter l as the weighted average of the range of diameters, namely,

$$l=\sum _{i=1}^6{p}_i{d}_i,$$

(11)

where $d_i$ denotes an average fraction and the percentages $p_i$ of the fractions are used as the weights of the average.

The reference porosity is calculated as in [43,44,45,46,47] with the coefficient of grain uniformity

$$\eta =\frac{d_{60}}{d_{10}}$$

(12)

and the empirical relationship between $\varphi$ and $\eta$, namely,

$$\varphi =0.255(1+0.{83}^{\eta }).$$

(13)

For this grain-size distribution, it results in the reference porosity ${\varphi }_{ref}=0.35$. The reference hydraulic conductivity is calculated after the Terzaghi equation [43]

$$k_f=C_t{\left(\frac{\varphi -0.13}{\sqrt[3]{1-\varphi }}\right)}^2{d}_{10}^2,$$

(14)

with $C_t$ as the sorting coefficient, ranging between $6.1\times {10}^{-3}$ and $10.7\times {10}^{-3}$. Here, we choose the average value of $8.4\times {10}^{-3}$. Based on the the grain-size distribution from Table 1 and the porosity calculated with Equation (13), the resulting reference permeability is ${\mathbf{K}}_{ref}=1.01\times {10}^{-9}{\mathrm{m}}^2$.

The entry pressure $p_d$ [Pa] of the Brooks–Corey capillary pressure–saturation relationship,

$$p_c\left(S_e\right)=p_d{S}_e^{MM-1/\lambda },$$

(15)

is adjusted with the reference entry pressure $p_{d,ref}=1000\mathrm{Pa}$ plus the spatially varying porosity and permeability using the Leverett J-function [48] as

$$p_d=p_{d,ref}\sqrt{\frac{\varphi {\mathbf{K}}_{ref}}{{\varphi }_{ref}\mathbf{K}}}.$$

(16)

From Equations (9), (10), and (16), we have calculated spatial distributions for the porosity, permeability, and capillary entry pressure. As macropores, we insert additional areas with very high porosity values (dimensions shown in

Table A1. Macropore parameters for the three implemented scenarios.

		Porosity [-]	x1 [m]	x2 [m]	y1 [m]	y2 [m]
Scenario (a)	Lens1	-	-	-	-	-
	Lens2	-	-	-	-	-
Scenario (b)	Lens1	0.94	3.0	3.2	0.0	10.0
	Lens2	0.94	7.2	7.6	0.0	10.0
Scenario (c)	Lens1	0.90	3.6	4.0	4.0	10.0
	Lens2	0.94	7.4	7.6	0.0	10.0

), and finally, we end up with the porosity fields, as shown in Figure 1.

We assign Neumann no-flow boundary conditions at the right and left side of the domain. At the top boundary, a Neumann boundary condition represents the water in-flowing rain. Its mass is calculated from the official threshold for heavy-rain events $q=25\mathrm{l}{\mathrm{m}}^{-2}{\mathrm{h}}^{-1}=0.007\mathrm{kg}{\mathrm{m}}^{-2}{\mathrm{s}}^{-1}$ provided by the German Weather Service [49]. At the lower boundary, we set Dirichlet conditions with a fixed pressure and saturation.

The residual water saturation for the domain is set to $S_{wr,d}=0.1$, except for the macropores. Several studies [23,50,51] assume the residual saturation in the macropores to zero, namely, $S_{wr,m}=0.0$.

In Figure 2, we show the porosity distribution histograms for the respective scenarios (a), (b), and (c). The spikes at the tails of the latter two represent the high porosities assigned to the macropores. Obviously, this confirms the bimodal porosity behavior of macroporous domains.

4. Results and Analysis

4.1. Flow-Regime Characterization

In flow problems, it is important to assess which flow regime dominates. As already mentioned, in this study, we use the dimensionless Re and Fo numbers as indicators for the prevailing flow regimes. Given their values, the required flow model can be chosen accordingly. This is why for the phase-velocity value, we use the cell data results from the Darcy models.

The Reynolds number Re is based on the average grain diameter (l) of Table 1 and the phase velocity,

$$Re=\frac{\left|{\mathbf{v}}_{\alpha }\right|{\rho }_{\alpha }l}{{\mu }_{\alpha }}[-],$$

(17)

which represents the ratio of inertial force to viscosity force. The Forchheimer number Fo has an expression that includes multi-phase and inertia effects [6],

$$Fo=\frac{{\rho }_{\alpha }k{k}_{r\alpha }{\beta }_{\alpha }{\mathbf{v}}_{\alpha }}{{\mu }_{\alpha }}[-];\alpha \in \{w,n\},$$

(18)

where ${\beta }_{\alpha }$ is calculated from Equation (7).

The non-Darcy effect (E) is the error caused by ignoring the non-Darcy behavior. According to [34,35], it is defined as “the ratio of the pressure gradient consumed in overcoming the liquid-to-solid interactions to the total pressure gradient” and can be calculated from [34] as

$$E=\frac{Fo}{1+Fo}.$$

(19)

From this equation, we understand that the Forchheimer number is directly dependent on E. If we denote $E_c$ as the critical value for the non-Darcy effect, then the critical Forchheimer number would be

$$F{o}_c=\frac{E_c}{1-E_c}.$$

(20)

For $E_c=10\%$ [4,34], the critical Forchheimer number, or the lower limit for inertial flow to be considered, is equal to $F{o}_c=0.11$. This value is larger than the one given from literature and will be used as our critical Fo.

In Figure 3, we show the Reynolds and Forchheimer numbers for the three scenarios at the end of the simulation time. In Scenario (a), the Reynolds numbers barely reach up to $Re=1$, with the majority of the domain in a smaller range of values. In Scenarios (b) and (c), where macropores are added, it reaches up to $Re=2.5$. Again, the majority of the domain behaves similar to Scenario (a), but we know that this variation of values at the tails of the histograms is caused by the macropores.

According to the typical range of the Forchheimer number for transition to inertial flow from the literature, all scenarios depicted here belong in the inertial regime. Even if we take into account the “stricter” $F{o}_c=0.11$, we reach this limit already in Scenario (a). In Scenarios (b) and (c), the Forchheimer number reaches up to $Fo=0.35$. We can see that with both dimensionless criteria, we have a large spread of values on the right tail for the macropore scenarios, along with some small spikes. We assume that these belong to the macroporous regions of the domain.

The results for the Forchheimer number confirm that the inertial effects are strong enough and they cannot be ignored in all three scenarios. The Reynolds criterion $Re<1$ is only valid for Scenarios (b) and (c) in regions where flow is not affected by the macropores.

4.2. Velocity Streamlines

In Figure 4, we show the streamlines representing the water-phase velocity over the porosity fields for the two macropore scenarios after the last timestep. We present only the results for the Forchheimer–Case 1 model, as no significant differences between the others are visible. The strongest contribution of the macropores to the flow pattern is the preferential flow path they obviously offer. In both scenarios, we notice that the velocity streamlines around the macropores lead into them. This is caused by

• The different geometrical characteristics–continuous length and large width,
• The larger porosity and permeability compared to the rest of the domain, and
• The residual saturation, which is $S_{wr,m}=0.0$.

Moreover, the water-phase velocity in the macropores is considerably higher than in the rest of the domain. As we can see, in both scenarios, the highest values occur in the thinner macropores, where we observe also a slightly higher pressure. Between the three models used to simulate flow, the Forchheimer–Case 2 is the one that differs the most compared to the other two. We assume that the two cases differ the most in regard to the approach and consideration of the inertial effects. These should have the largest impact when the velocity is changing. To distinguish the two cases, we look into the green highlighted cells in Scenario (c), where the velocity is slowed down at the end of the left (wide) macropore.

According to [52,53], the microscopic inertial effects at high velocities distort the flow lines and, therefore, increase the gradients of velocity and the pressure drop. Additionally, ref. [35] explains this with the dissipation of energy as fluid particles accelerate in smaller pores and decelerate while entering large pore spaces. As we study our cases, there should be a visible faster acceleration or deceleration of water flow for the case with the higher inertial effect.

In

Table 2. Model liquid phase velocity results [$\mathrm{m}{\mathrm{s}}^{-1}$] for scenario (c) of Forchheimer–Case 1.

	Left Cell	Right Cell
Macropore cells	2.47 × ${10}^{-5}$	2.46 × ${10}^{-5}$
Domain cells	1.87 × ${10}^{-5}$	1.82 × ${10}^{-5}$

and

Table 3. Model liquid phase velocity results [$\mathrm{m}{\mathrm{s}}^{-1}$] for scenario (c) of Forchheimer–Case 2.

	Left Cell	Right Cell
Macropore cells	$2.75\times {10}^{-5}$	$2.75\times {10}^{-5}$
Domain cells	$1.85\times {10}^{-5}$	$1.80\times {10}^{-5}$

, the velocities of the green highlighted cells in Figure 4 are shown. The cells are the last in the macropore and the first afterwards. Generally, we can say that Case 2 has higher velocities in the macropores. Additionally, we can see that the velocity in Case 2 drops lower in the domain cells. If we follow the previously referenced literature, it would imply that Case 2 accelerates and decelerates more than Case 1 and therefore has more inertial effects.

4.3. Breakthrough Curves

Figure 5 and Figure 6 show the breakthrough curves for all three scenarios. Breakthrough stands here for the time at which the water saturation in the lowest domain cell reaches a threshold of $ϵ=1\times {10}^{-5}$ above the residual saturation. For the Scenarios (b) and (c), the macropores are dominating, which is best displayed on a logarithmic timescale. In the macropores themselves, there are also no significant differences, especially because the arrival times are so small.

The differences between the cases can be seen in Scenario (a) and the left macropore of Scenario (c), as shown enlarged in Figure 6. There, the Darcy model and the Forchheimer–Case 1 are very similar, while Case 2 clearly differs. In the areas where the water reaches the bottom faster, Case 2 is slower, while in the areas with a slower water front, Case 2 is faster. Since Case 1 is so similar to the Darcy model, this supports the approach of Case 2 having the Forchheimer coefficent with a higher inertia effect. This was also tested with higher influx boundary conditions, i.e., $q=0.014\mathrm{kg}{\mathrm{m}}^{-2}{\mathrm{s}}^{-1}$. There, the differences between the models are more significant, and the results are even more clear with regard to the acceleration and deceleration.

4.4. Numerical Behavior

In

Table 4. Computational effort and numerical behavior for Scenario (a).

	Total Time [s]	No. Global Timesteps	No. Newton Iterations	No. Failed Timesteps
Darcy	57.08	121	1107	2
Forchheimer–Case 1	202.02	211	2009	1
Forchheimer–Case 2	1104.53	1274	12,583	0

,

Table 5. Computational effort and numerical behavior for Scenario (b).

	Total Time [s]	No. Global Timesteps	No. Newton Iterations	No. Failed Timesteps
Darcy	42.02	86	787	0
Forchheimer–Case 1	323.86	316	3019	9
Forchheimer–Case 2	905.74	1032	10,159	6

and

Table 6. Computational effort and numerical behavior for Scenario (c).

	Total Time [s]	No. Global Timesteps	No. Newton Iterations	No. Failed Timesteps
Darcy	48.30	98	907	0
Forchheimer–Case 1	290.04	316	2991	11
Forchheimer–Case 2	1003.27	1162	11,442	7

, the total runtime of the model simulations, the timesteps required for the solution, and the total of Newton iterations to solve the timesteps are listed for the implemented scenarios.

A general observation for the Darcy and Forchheimer–Case 2 models is that the total runtime of the simulations decreases when the macropores, which induce faster flow, are added in the domain. More specifically, for both of them, the lowest total runtime occurs for Scenario (b) (Figure 1). This may be explained by macropores being large areas in the domain where porosity, permeability, and entry pressure are constant.

For the Forchheimer–Case 1 model, the total simulation runtime increases when the macropores are added (Table 5 and Table 6). Interestingly, the longest simulation for this model is for Scenario (b), for which both the other models have their shortest runtime. This will be subject to further in-depth investigations beyond the scope of this study.

4.5. Root-Mean-Square Error (RMSE)

For our goal of comparing the two Forchheimer models and understanding their possible differences, we calculate the RMSE (root-mean-square error) for key parameters of the simulations (saturation [-], mobility [${\mathrm{Pa}}^{-1}{\mathrm{s}}^{-1}$], and velocity [$\mathrm{m}{\mathrm{s}}^{-1}$]).

The RMSE is used to measure the error or the deviation of model-predicted results from observed values [54]. In our case, we want to measure the difference between the Forchheimer–Case 1 and Forchheimer–Case 2 results from the Darcy results. We use the Darcy results as our reference values and the results from the Forchheimer cases as the new values that need to be examined. We use the data of the final timestep.

The RMSE is calculated from the following equation [54]:

$$\mathrm{RMSE}=\sqrt{\sum _{i=1}^n\frac{{(\widehat{y_i}-y_i)}^2}{n}},$$

(21)

where $\widehat{y_i}$ are the new values, $y_i$ are the reference values (Darcy model results), and $n=2500$ is the number of cell data pairs.

We compare the Darcy and Forchheimer models, using three different response parameters with different units or no units and different scales [54]. That is why we calculate a normalized root-mean-square error (NRMSE), where we have a relative rather than an absolute squared difference. We calculate our NRMSE by dividing the squared differences with the reference values [54]:

$$\mathrm{NRMSE}=\sqrt{\sum _{i=1}^n\frac{{\left(\frac{\widehat{y_i}-y_i}{y_i}\right)}^2}{n}}.$$

(22)

From the results of the comparison between the Darcy and the Forchheimer–Case 1 models (

Table 7. NRMSE values for the mobile (liquid) phase parameters between the Darcy and Forchheimer–Case 1 models (t = TEnd).

	NRMSE	Saturation [-]	Mobility [-]	Velocity [-]
Case 1
Scenario (a)		2.61 × ${10}^{-5}$	4.63 × ${10}^{-4}$	9.50 × ${10}^{-5}$
Scenario (b)		0.12 × ${10}^{-2}$	1.90 × ${10}^{-2}$	0.49 × ${10}^{-2}$
Scenario (c)		0.16 × ${10}^{-2}$	2.33 × ${10}^{-2}$	1.28 × ${10}^{-2}$

), we can see that the differences of Scenario (a) are a lot smaller compared to the rest and the respective differences of the Darcy and Forchheimer–Case 2 models. This would mean that the Darcy and Forchheimer Case 1 models are more in “agreement” for a heterogeneous domain when no macropores are added. For the rest of the scenarios, in both model differences, we notice a similar tendency with the largest differences for the water-phase velocity and water-phase saturation parameters. The differences in saturation are a result of the propagation of the flow field.

5. Discussion

In this study, we investigated flow in an unsaturated heterogeneous domain with and without macropores. In order to better understand the mechanisms that affect flow, we implemented three different model approaches: an extended multi-phase Darcy model that simulates creeping flow and two multi-phase Forchheimer models that consider the possible inertia effects. The Forchheimer models include a Forchheimer (inertia) coefficient, which is subject of our investigation. Its calculation is mainly empirical, and the theoretical equations that exist vary depending on which parameters of the domain and the phases are included. Here, we use two equations to calculate it according to [6,33].

Regarding the final steady-state flow field, we can see from the breakthrough curves (Section 4.3), as well as from the velocity streamlines (Section 4.2), that the three models do not differ substantially. This shows that the two-phase flow problem is affected by the domain’s characteristics (porosity, permeability), the traversing fluid’s characteristics (saturation, relative permeability, mobility), and their interaction (capillary pressure).

The differences in the final velocity values are a result of the different approaches to the solution of velocity. The Forchheimer coefficient based on Equation (5) depends on permeability and phase saturation via the relative permeability. On the other hand, in Equation (7), the coefficient is directly proportional to tortuosity and inversely proportional to permeability, relative permeability, porosity, and saturation. This is reflected in Figure 7, where the beta-coefficient maps of both Forchheimer cases for Scenario (a) are shown. The distribution of the beta coefficient for Case 1 resembles the distribution of porosity and permeability and shows the same spatially smooth behavior. Case 2, on the other hand, is spatially more disrupted, even though also the areas with very high and very low porosity are visible.

In Section 4.2 and Section 4.3, we supposed that the Forchheimer–Case 2 has the higher inertial effects. This is also visible in Figure 7, where the coefficient can become several orders of magnitude larger. This is probably influenced the most by the difference in how the permeability is considered in the calculation of the coefficient. In Case 1, it is with a square root, while in Case 2, it is not.

Since the context this study is related to water flow and the transport of low-level radioactive substances in landfills [17], the transient flow field is of high importance. This holds in particular since the differences between the breakthrough curves increase with higher velocities. Additionally, we keep in mind that in reality, macropores reaching from top all the way through to the bottom as in Scenario (b) should be unlikely. This puts more emphasis on the behavior of the coefficients in Scenarios (a) and (c).

Choosing the “correct” Forchheimer coefficient for our case is de facto impossible, since there are so many different possibilities with different parameters [6,24,29], for which we have currently no experimental data available yet. However, we can conclude that the coefficent for Case 1 needs less computing time and represents less inertial effects. In particular, the differences to the Darcy model for Scenario (a), shown in

Table 8. NRMSE values for the mobile (liquid) phase parameters between Darcy and Forchheimer–Case 2 models (t = TEnd).

	NRMSE	Saturation [-]	Mobility [-]	Velocity [-]
Case 2
Scenario (a)		1.50 × ${10}^{-2}$	2.27 × ${10}^{-1}$	5.03 × ${10}^{-2}$
Scenario (b)		1.18 × ${10}^{-2}$	2.10 × ${10}^{-1}$	6.07 × ${10}^{-2}$
Scenario (c)		1.53 × ${10}^{-2}$	2.40 × ${10}^{-1}$	5.52 × ${10}^{-2}$

, are marginal. If we follow our thought from above that the most interesting scenario is probably Scenario (a), we may conclude that using the Forchheimer equation with this coefficient is not appropriate and justified. It results in more computing time while providing no significant difference, even though we are above the critical Forchheimer number.

6. Conclusions

Motivated by flow events in heterogeneous landfills from debris of dismantled nuclear power plants, we investigated the differences between Darcy’s law and the Forchheimer extension in the presence of macropores. We reach the following conclusions:

• Based on numerical scenarios, which all featured Reynolds numbers beyond the validity of Darcy’s law and, thus, an expected influence of inertia effects, we found that different approaches to calculating Forchheimer were showing different behavior in representing inertia effects.
• At the same time, the different investigated Forchheimer coefficients are associated with different computational efforts.
• Dependent on the inflow boundary condition, the resulting velocities, and the size of the investigated domain, we could show that heterogeneous setups including macropores require the Forchheimer extension to account for inertial effects for the transient transport solution.
• Macropores dominate the flow-field and the breakthrough times, in particular when they reach over large vertical distances, in this study from top to bottom.

While the present study improves our general understanding of the influence of macropores on the flow behavior by means of synthesized scenarios, further work has to be performed for achieving predictability in practice. In particular, we will extend the spatial dimensions to 3D and to the scale of a realistic landfill. Additionally, the transport of radionuclides and their decay chaines will be considered, where the results of this study can be used to optimize the applied flow models. This further involves concepts for modeling the sorption of radionuclides. Finally, we expect to acquire experimental data to calibrate and validate the models.