Effective Forchheimer Coefficient for Layered Porous Media

Inertial flow in porous media, governed by the Forchheimer equation, is affected by domain heterogeneity at the field scale. We propose a method to derive formulae of the effective Forchheimer coefficient with application to a perfectly stratified medium. Consider uniform flow under a constant pressure gradient ΔP/L\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta P/L$$\end{document} in a layered permeability field with a given probability distribution. The local Forchheimer coefficient β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} is related to the local permeability k via the relation β=a/kc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta =a/k^c$$\end{document}, where a>0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a>0$$\end{document} being a constant and c∈[0,2]\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c\in [0,2]$$\end{document}. Under ergodicity, an effective value of β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} is derived for flow (i) perpendicular and (ii) parallel to layers. Expressions for effective Forchheimer coefficient, βe\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _e$$\end{document}, generalize previous formulations for discrete permeability variations. Closed-form βe\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _e$$\end{document} expressions are derived for flow perpendicular to layers and under two limit cases, F≪1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F\ll 1$$\end{document} and F≫1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F\gg 1$$\end{document}, for flow parallel to layering, with F a Forchheimer number depending on the pressure gradient. For F of order unity, βe\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _e$$\end{document} is obtained numerically: when realistic values of ΔP/L\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta P/L$$\end{document} and a are adopted, βe\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _e$$\end{document} approaches the results valid for the high Forchheimer approximation. Further, βe\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _{e}$$\end{document} increases with heterogeneity, with values always larger than those it would take if the β-k\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta - k$$\end{document} relationship was applied to the mean permeability; it increases (decreases) with increasing (decreasing) exponent c for flow perpendicular (parallel) to layers. βe\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _{e}$$\end{document} is also moderately sensitive to the permeability distribution, and is larger for the gamma than for the lognormal distribution.


Introduction
An understanding of the interplay between nonlinear effects in porous media flow and domain heterogeneity is of great importance in several engineering and geological disciplines. Darcy's law, usually adopted to model fluid flow in porous media, implies a linear relationship between flow rate and pressure (or head) gradient. However, for high flow rates, the pressure gradient is higher than that predicted by Darcy's law, the deviation increasing with the flow rate. This phenomenon, known as non-Darcy or Forchheimer flow, occurs in several civil, environmental, and industrial engineering applications, such as flow in rockfill dams; flow in coarse-grained, fractured or karstic porous media; flow in the vicinity of pumping or injection wells; gas flow in natural or artificial porous media; industrial filtration processes; reservoir exploitation. The importance of the aforementioned applications generated a large body of scientific and technical literature in the past decades aimed at understanding the mechanisms governing non-Darcy flow in porous media. The nonlinear relationship between flow rate and pressure drop is usually attributed to the insurgence of inertial effects within the laminar flow regime; these inertial effects are commonly represented by adding to Darcy's equation an additional term proportional to the fluid density and to the second power of the flow rate; the coefficient of proportionality is known as Forchheimer, or inertial, or non-Darcy coefficient. For a review of the different existing formulations of Forchheimer's law see, among others, Trussell and Chang (1999), Sidiropoulou et al. (2007), and Huang and Ayoub (2008). An alternative to Forchheimer's law is the Izbash equation, which states that the hydraulic gradient is a power function of the specific discharge (Bordier and Zimmer 2000;Moutsopoulos et al. 2009); more complex polynomial equations also exist (Balhoff et al. 2010;Lofrano et al. 2020), based on the fact that inertial flow in porous media may be classified into different subregimes depending on the form of the inertial correction (Agnaou et al. 2017).
In field-scale modeling of porous media flow at scales ranging from local to regional, heterogeneity in model parameters plays a crucial role; as a consequence, a relevant effort was devoted in the scientific literature at characterizing heterogeneity in hydraulic and transport properties, and at deriving representative properties as functions of statistical parameters describing heterogeneity. While most of the effort was directed at hydraulic conductivity (Sanchez-Vila et al. 2006), heterogeneity in the Forchheimer coefficient was also recognized in the field (Jones 1987;Narayanaswamy et al.. 1999), prompting researchers to investigate the concept of a representative Forchheimer coefficient at a given upscaled scale, as a function of its local value (Fourar et al. 2005;Auriault et al. 2007;Garibotti and Peszynska 2009;Aulisa et al. 2014). The nonlinearity of the flow complicates the upscaling problem; related work in the context of geologic media includes the determination of the effective conductivity for non-Newtonian fluid flow (Di Federico et al. 2010;Airiau and Bottaro 2020).
The value of the Forchheimer coefficient has long been recognized as being empirically correlated to other properties of the porous medium, namely porosity, tortuosity, and permeability; a recent summary of proposed correlations was provided by Arabjamaloei and Ruth (2017). In particular, the correlation of the Forchheimer coefficient with permeability k was found to be a power-law inverse one since the work of Ergun (1952) and Geertsma (1974). The present paper incorporates the − k correlation into deriving an upscaled Forchheimer coefficient for flow under a uniform pressure gradient in a stratified porous medium.
The Forchheimer seepage law, related experimental findings, and the nature of the − k correlation are illustrated in Sect. 2, while effective parameters for Forchheimer flow in one-dimensional stratified porous media are derived in Sect. 3 for flow perpendicular and parallel to layers; in the latter case, two analytical approximations are presented for lowand high-Forchheimer number case, together with a general option for numerical evaluation. These results are discussed as a function of problem parameters in Sect. 4. A set of conclusions closes the paper (Sect. 5), while details on the permeability distribution functions adopted and special values of the quantities of interest are presented in Appendices A through D.

Forchheimer Flow Law
For high seepage velocities, Darcy's law is inadequate to represent fluid flow in porous media due to inertial effects, which are no longer negligible when the pore-scale Reynolds number exceeds a threshold value between 1 and 10 (Bear 1979), while turbulence usually occurs at much higher values of the pore Reynolds number (Huang and Ayoub 2008). In turn, inertial effects derive from the nonlinear terms in the general momentum balance equation.
The inertial effects are commonly represented by adding to Darcy's law an additional term proportional to the fluid density and the second power of the flow rate; the resulting equation is termed Forchheimer's flow law; its isotropic form reads (Auriault et al. 2007) where P = p + gz is the generalized pressure including gravity effects, p the pressure, and the fluid density and dynamic viscosity, the specific discharge vector [LT −1 ] , k F the Forchheimer permeability coefficient [L 2 ] , larger than, but close to, the Darcy permeability coefficient k D appearing in Darcy's law (El-Zehairy et al. 2019). Finally, is the Forchheimer coefficient [L −1 ] , also termed velocity coefficient, inertial coefficient, or non-Darcian coefficient; when = 0 , Eq. (1) reduces to Darcy's law with k D in place of k F . In the following, we will assume k D = k F = k.
All experimental values cited above pertain to water flow, while values of found experimentally for gas flow tend to be higher by a couple of orders of magnitude: values of the Forchheimer coefficient measured by Jones (1987) on a total of 364 core plugs were in the range 10 5 to 10 13 m −1 ; Zeng and Grigg (2006) performed laboratory experiments with nitrogen gas, finding = 2.88 × 10 8 − 1.57 × 10 10 m −1 ; Wells et al. (2008) obtained for airflow values of in the range 1.41 × 10 6 − 5 × 10 8 m −1 with a field permeameter. For fracture proppant packs in hydrofracturing, lies in the range 7.2 × 10 4 − 3.8 × 10 7 m −1 (Friedel and Voigt 2006).

Empirical ˇ− k Correlations
Several researchers found a correlation to exist between local values of the Forchheimer coefficient and other properties of the porous medium; the different formulations of this relationship, either theoretically or empirically based, are summarized by Li and Engler (2001); according to them, a good representation of most of the previous work, either of theoretical or experimental nature, is provided by the formulation in which is porosity, is tortuosity, and c 1 , c 2 , and c 3 are three experimental constants for an assigned porous medium. Similar formulations are reported in Skjetne et al. (2001) and Saboorian-Jooybari and Pourafshary (2015), albeit with = 1 . The former authors note that as [k] = [L 2 ] and [ ] = [L −1 ] , media that are scaled copies of each other have c 2 = 0.5 . According to the analysis of Geertsma (1974), conducted both via dimensional and field data analysis, c 2 = 0.5 ; according to Narayanaswamy et al. (1999), who also cites previous experimental work, c 2 = 1.25 . Jones (1987) analyzed a large number of field samples with differing lithologies, obtaining c 2 = 1.55 ; the reticular model of Thauvin and Mohanty (1998) and the experimental data of Cooper et al. (1998) suggest c 2 = 1 . The values of c 2 in Eq. (2), according to the review by Li and Engler (2001), are generally in the range of 0.5 to 1.88. The same range is reported in the more recent review of Saboorian-Jooybari and Pourafshary (2015). In the sequel of this note, aimed at investigating the impact of heterogeneity in the permeability distribution on the effective Forchheimer coefficient, porosity , and tortuosity will be considered as constants, under the hypothesis that permeability varies across several order of magnitudes while porosity and tortuosity do not (Bear 1979); this allows writing Eq. (2) as where c = c 2 [−] is a non-negative empirical exponent, which will be assumed to lie in the range 0 − 2 , and the empirical parameter a = a(c) has dimensions [L (2c−1) ] ; the case c = 0 is equivalent to assuming a value of the local-scale Forchheimer coefficient independent from permeability.
The value of the parameter a is highly sensitive to the exponent c and the medium properties, porosity and tortuosity . For c = 0.5 , Ergun (1952) for sand packs and Macdonald et al. (1979) and consolidated sandstone. For c = 1 , the experimental data of Liu et al. (1995) result in a = 71274432 m considering as mean values of porosity and tortuosity = 0.25 and = 2 . For c = 1.55 , Jones (1987) obtained a = 187724450 m 2.1 for limestone, and fine-grained sandstone considering = 1 . In what follows, we assume the latter value of a to be valid also for the case c = 1.5 as a = 187724450 m 2 . Additional, different experimental values of a can be found in Saboorian-Jooybari and Pourafshary (2015) for different correlations of and k.

Effective Parameters for Forchheimer Flow in Stratified Porous Media
Consider an heterogeneous, perfectly stratified porous domain of length L, subject to an external, generalized pressure difference ΔP = P(0) − P(L) = P 1 − P 2 in the x-direction; the domain permeability k is taken to vary as a weakly stationary random field characterized by its probability density function (PDF) f(k). The ergodic assumption is assumed to hold; hence spatial averages and ensemble averages are interchangeable, and a single realization can be examined. To derive an expression for the effective permeability k e and effective Forchheimer coefficient e in perfectly layered media, two limiting cases are examined, coincident with those providing the lower and upper bound for the effective permeability (Matheron 1967): (i) flow perpendicular to the layering (serial-type layers), when the pressure gradient is parallel to the permeability variation; (ii) flow parallel to the layering (parallel-type layers) when the pressure gradient is transverse to the permeability variation, as shown schematically in Fig. 1a and b, respectively. Such perfect layering is an idealization of real-world conditions frequently employed in the stochastic hydrology literature (Dagan 2017).

Flow Perpendicular to Layers
For flow perpendicular to layers, the domain is made of N homogeneous layers of equal thickness Δx , each having an area A perpendicular to the flow direction; the i-th layer has a permeability of k i and a Forchheimer coefficient i = i (k i ) given by Eq. (3a-b). By virtue of mass conservation, volumetric flux Q = qA through each layer is the same; assuming that the one-dimensional version of the flow law given by Eq.
(1) holds locally, the pressure difference in each layer is provided by Since Δx = L∕N and ΔP = ∑ n i ΔP i , summing over the layers and then taking the limit as N → ∞ , the length of each layer tends to zero and the discrete permeability variation to a continuous one; then under ergodicity where ⟨⋅⟩ is the expected value, and the relationship = (k) given by Eq. (3a-b) was employed. Eq. (5) extends to a continuous permeability distribution the discrete equation of Narayanaswamy et al. (1999). The result on the r.h.s. is valid unless f(k) is a constant, e.g. for a uniform distribution: see Appendix 1 for this case.
The effective parameters k es and es for flow perpendicular to layers are defined by Comparison between Eqs. (6) and (5) gives for the effective parameters where k H is the harmonic mean.

Flow Parallel to Layers
In this case the porous domain is still constituted by N homogeneous layers, with the i-th layer having an area ΔA = A∕N perpendicular to the direction of the external pressure gradient, a permeability k i and a Forchheimer coefficient i = i (k i ) according to Eq. (3a-b). The volumetric flux ΔQ i in the i-th layer is, solving the quadratic Eq. (1) in terms of flow rate Summing over the layers and switching to a continuous parameter variation gives where the quantity is the dimensionless external pressure gradient, essentially a Forchheimer number, the ratio of liquid-solid interaction to viscous resistance (Ruth and Ma 1992). Substituting the expression of (k) given by Eq. (3a-b) gives F = (4 ak 2−c ∕ 2 )(ΔP∕L) , and Eq. (9) is not amenable to closed-form integration for the most commonly adopted formulations of the probability density function f(k), except for the uniform distribution, see Appendix C. (4) For c = 2 , Eq. (9) with Eq. (3a-b) gives while the flowrate in terms of the effective parameters is Comparing Eqs. (11) and (12) leads to where k A is the arithmetic mean.
For c ≠ 2 , we note that for F both small and large with respect to unity, an approximate closed-form solution for the flow rate and the effective Forchheimer coefficient can be found explicitly. In particular, for F ≪ 1 (case 1), where the quantity F is small, but not so small that inertial effects on the flow could be disregarded, the quantity √ 1 + F in Eq. (9) can be approximated by its second-order Taylor's expansion . When, on the other hand, F ≫ 1 (case 2), we assume in Eq. (9) at leading order Fourar et al. 2005), i.e. the Forchheimer term in Eq. (9) is dominant with respect to the Darcy one.
Both approximations examined may be appropriate for laboratory or field applications, depending on the fluid involved, the values of the porous medium properties, and the pressure gradient. Typically, laboratory and field measurements involve water or gas flow. In the first case, taking for the fluid properties those of water at 10 • C ( = 999.7 kg∕m 3 , = 1.337 × 10 −3 kg∕m∕s ), and further assuming k = 10 −8 m 2 , = 50 m −1 , ΔP∕L = 10 4 Pa ⋅ m (corresponding approximately to a unit hydraulic gradient), yields F = 0.117 , while larger values of permeability, Forchheimer coefficient, and/or pressure gradient bring about larger values of F. The data reported in Moutsopoulos et al. (2009) lead to F = 1.82 − 2025 (depending on the actual temperature during the experiments), those of Bordier and Zimmer (2000) to F = 0.80 − 44.23 , those by Wahyudi et al. (2002) to F = 0.16 − 5.32 . In applications involving gas flow, typically the values of the ratio ∕ 2 , the Forchheimer coefficient , and pressure gradient ΔP∕L are decidedly higher when compared to water flow, while the permeability values are lower. Taking = 1.25 kg∕m 3 , = 1.79 × 10 −5 kg∕m∕s for air at 10 • C , a porous medium with k = 10 −12 − 10 −10 m 2 and = 2 ⋅ 10 8 − 8 ⋅ 10 8 m −1 , yields for a pressure gradient in the range ΔP∕L = 10 6 − 10 7 Pa∕m a wide interval F = 75 − 1.87 × 10 6 . Values of F = 280 − 4.2 × 10 6 , consistently much higher than unity, are associated with the gas condensate field results of Narayanaswamy et al. (1999). The field air permeameter data of Wells et al. (2008) lead to values in the interval F = 0.03 − 27.97 , with five out of sixteen samples providing values of F less than one. For the laboratory measurements of Zeng and Grigg (2006), F = 4.5 × 10 −5 − 2.6 × 10 −3 , while for those of Cooper et al. (1998), F = 1.68 × 10 −3 − 1.28 × 10 −1 .

Low Forchheimer case ( F ≪ 1)
Inserting Eq. (3a-b) in Eq. (14) gives Now taking an upscaled equation like (6) written in terms of the effective parameters for parallel-type layers in this case, namely k ep1 and ep1 , solving it for Q = qA and adopting the same second-order Taylor's approximation employed in Eq. (14), gives Comparison between Eqs. (15) and (16) yields where k A is the arithmetic mean (Matheron 1967).

High Forchheimer Case ( F ≫ 1)
Performing the same steps as in the previous case and adopting in Eq. (9) the zero order approximation

Effective Parameters ˇe s , ˇe p1 , ˇe p2 for given Permeability Distribution
The effective parameters are evaluated for two given PDFs of permeability, the lognormal and the gamma distributions, both suitable to represent non-negative permeability distributions. Their expressions are illustrated in Appendices A and B, respectively.

Numerical Evaluation of Effective Forchheimer Coefficient ˇe pN for Flow Parallel to Layers and any F
Considering the specific flow rate q = Q∕A directly, Eq. (9) can be rearranged as For c = 1 , Eq. (28) is amenable to a simpler solution for the lognormal distribution and a closed-form solution for the gamma distribution; these are reported in Appendix D. In general, once q is evaluated numerically from Eq. (28) for a given permeability distribution, the numerical effective Forchheimer coefficient epN for the parallel-type layers will be Note from Eq. (29) that, in variance with Eq. (7a-b) for flow perpendicular to layers, the effective Forchheimer coefficient for flow parallel to layers depends on the boundary conditions.

Discussion of Results
The effective properties, permeability and Forchheimer coefficient, are functions of the parameters describing heterogeneity. While the results for effective permeability recover well-known results (Matheron 1967;Sanchez-Vila et al. 2006), those for the effective Forchheimer coefficient are novel. Some general tendencies of its behavior are evident from Eqs. (22a-b)- (27): (i) for flow perpendicular to layers, the effective Forchheimer coefficient increases with increasing heterogeneity; (ii) the same is true for flow parallel to layers and low Forchheimer number; (iii) the reverse is true for flow parallel to layers and high Forchheimer number; (iv) the shape of the distribution does not influence the aforementioned tendencies; (v) the parameter c influences the effective Forchheimer coefficient for serial and parallel arrangements in a non-trivial way.

Analytical Results for Flow Perpendicular or Parallel to Layers
We express Eqs. (22a-b)-(27) in the dimensionless form adopting as a scale for e the value of the effective Forchheimer coefficient would take if the correlation given by Eq. (3a-b) was applied to a representative permeability equal to the mean value ⟨k⟩ of the random permeability field. One then obtains for the lognormal distribution and gamma distribution the following expressions: Specific expressions of the effective Forchheimer coefficient for special values of c are listed in Table 1.
In the following, the dimensionless effective Forchheimer coefficient ̂e for flow perpendicular and parallel to layers is illustrated as a function of the permeability coefficient of variation C vk = k ∕⟨k⟩ for the lognormal and gamma distributions (see Appendix A and Appendix B). Note that for a lognormal distribution the classical upper limit of validity of the first-order approximation in stochastic hydrology, 2 y = 2 Y = 1 , corresponds to C vk = 1.31 , and conversely C vk = 1 is equivalent to 2 y = 0.693 . The deviation of ̂e from unity represents the relative error made upon adopting the correlation Eq. (3a-b) and the mean permeability to evaluate e . Figures 2, 3 and 4 show ̂e as a function of C vk for different values of exponent c and for the lognormal and gamma distributions, respectively for flow perpendicular to layers, parallel to layers with low Forchheimer number (case 1), and parallel to layers with high Forchheimer number (case 2). Figure 2a and b depict the effective Forchheimer coefficient for lognormal and gamma distributions for flow perpendicular to layers. The effective Forchheimer coefficient ̂e s increases as a function of increasing heterogeneity (higher C vk ); an increase in the value of c has a similar effect. For c = 0 (no correlation between and k), the effective Forchheimer coefficient coincides with the value calculated inserting the mean permeability in the empirical relationship Eq. (3a-b).
Even for flow parallel to layering (Figs. 3 and 4), as the heterogeneity grows, the value of ̂e p1 and ̂e p2 increases. However, unlike in serial geometry, the increase in the coefficient c corresponds to a smaller rise of ̂e p1 and ̂e p2 as C vk grows; in this case, a value c = 0 corresponds to the maximum difference between the value of the actual Forchheimer coefficient and its local value. The values of ̂e p2 are smaller than that of ̂e p1 , showing that the high Forchheimer case implies smaller upscaled e . We remark 1 1 that the estimation of high Forchheimer values for c = 0.5 and c = 1.5 , shown by a solid line and black dots in Fig. 4a, are identical, as is evident from Eq. (33). Finally, upon comparing the results relating to the two adopted distributions in the two limit cases of low and high Forchheimer numbers, it is noted that with the same C vk and the value of the exponent c, for flow perpendicular to layering the values of e for the lognormal distribution are lower than those for the gamma distribution; the reverse is true for flow parallel to layering. The differences between distributions are modest for low-Forchheimer numbers and more significant for high-Forchheimer numbers. We note that in the low Forchheimer number regime (see Fig. 3), the calculated value of e for c = 1 is identical for the lognormal and gamma distributions.

Numerical Estimation of Parallel Flow
We assume water as the reference fluid ( = 1000 kg∕m 3 , = 10 −3 Pa ⋅ s ), and a mean permeability ⟨k⟩ = 10 −10 m 2 to perform the numerical integration of Eq. (28). The realistic values of a = 361480 [−], 71274432 m, 187724450 m 2 are considered for c = [0.5, 1.0, 1.5] , respectively, as discussed in Sec. 2.2. In addition, the pressure gradient is taken to be ΔP∕L = 400 Pa∕m . We then introduce a further set of dimensionless parameters for the numerical estimation of the Forchheimer coefficient of the parallel-layer case as where q , ΔP L , and â are normalized flow rate, pressure gradient, and proportionality constant a, respectively. Then, Eq. (28) can be written in terms of the dimensionless parameters of Eqs. (30a-b) and (37a-c) as Eq. (38) can be specialized for the lognormal distribution as and for Gamma distribution by substituting Eq. (B3) and C vk = k ∕⟨k⟩ The dimensionless effective Forchheimer coefficient for the parallel-type layers obtained from the direct numerical analysis is derived by substituting the dimensionless parameters of Eq. (30a-b) and Eqs. (37a-c) into Eq. (29), as: The outcomes of the direct numerical integration are shown in Fig. 5 for ΔP L = 0.0004 , c = 0.5, 1.0, 1.5 , and their corresponding a values, where the left and right panels depict the lognormal and gamma PDFs, respectively. For c = 0.5 and for both distribution functions, the effective Forchheimer coefficient obtained from the direct numerical integration ̂e pN approaches to the high Forchheimer number estimation ̂e p2 . For the cases c = 1.0 and c = 1.5 , the numerical results are almost identical to the high Forchheimer approximation.
We then perform a parametric study on the pressure gradient and proportionality parameter a to analyze their impact on the numerical evaluation of the effective Forchheimer coefficient for flow parallel to layers considering c = 0.5 , the most common value in the literature (Ergun 1952). Figure 6a and b show the effective Forchheimer coefficient for ΔP L = 0.0001, 0.0004, 0.0020 (or equivalently ΔP L = 100, 400, 2000 Pa∕m ) upon assuming a constant value of â = 361480.
For the larger pressure gradients, the numerical estimation approaches the high Forchheimer approximation ( F ≫ 1 ). Figure 6c and d highlight the variation of â = 36148, 361480, 3614800 where ΔP L = 0.0004 (or equivalently ΔP∕L = 400 Pa∕m ) is assumed to be constant. In this case, larger â values produce ̂e pN values very close to the high Forchheimer approximation. The trends are similar for both distribution functions, while the actual values of ̂e pN are higher for the lognormal distribution.

Conclusions
Nonlinear seepage in heterogeneous, perfectly layered porous media is investigated in this work with the aim of deriving expressions for the effective permeability and effective Forchheimer coefficient. The local Forchheimer coefficient is expressed as a function of the local permeability value with an inverse power relationship of exponent c and constant of proportionality a. The effective permeability for flow perpendicular or parallel to layers coincides respectively with the harmonic and arithmetic mean, and therefore decreases or increases as the heterogeneity of the domain increases. The increase in heterogeneity determines an increase in the effective Forchheimer coefficient for all setups and dimensionless parameters. Its value is always larger than that it would take if the − k correlation was applied to the mean permeability. The effective increases (decreases) with increasing (decreasing) exponent c for flow perpendicular (parallel) to layers. The influence of the permeability Effective dimensionless Forchheimer coefficient ̂e p for flow in heterogeneous layers arranged in parallel as a function of the permeability coefficient of variation C vk , for c = 0.5, 1.0, 1.5 , for lognormal (left panels) and gamma distribution (right panels). The numerical results obtained from the integration of Eq. (41) are shown compared with their corresponding low and high Forchheimer number approximation distribution adopted is most evident for flow perpendicular to layers and the high Forchheimer limit for flow parallel to layers: in both cases, the gamma distribution yields larger values of the effective Forchheimer coefficient than the lognormal. When realistic values of the pressure gradient and of the a parameter are adopted, the effective Forchheimer coefficient for flow parallel to layers is very close or identical to the result of the high Forchheimer approximation ( F ≫ 1 ). For flow parallel to layers, the effective Forchheimer coefficient depends on boundary conditions, in variance with the effective permeability. The outcomes of this work, despite model simplifications, provide additional insight into nonlinear effects on flow in heterogeneous porous media, emphasizing the sensitivity of the effective Forchheimer coefficient to heterogeneity.
A desirable extension is the determination of the effective Forchheimer coefficient for two-and three-dimensional geometry, to be obtained either: (i) via a conjecture, assuming that the effective Forchheimer coefficient for a 2-D isotropic domain is the weighted average (geometric mean) of the results obtained for the two limit geometries (flow Left panels represent the outcomes of the lognormal distribution, and right panels represent gamma distribution results. Low ( ̂e p1 ) and high Forchheimer approximation ( ̂e p2 ) results are over-imposed to the figures perpendicular and parallel to layers) examined here; the conjecture extends to the Forchheimer coefficient the procedure classically adopted for the effective permeability, equal to the geometric mean in 2-D; (ii) by composing the 1-D expressions to model flow geometries at the field scale that are intermediate between the two limit situations of motion parallel and perpendicular to the medium stratification; or (iii) by means of appropriate formal developments in line with the theory of composites. In all cases, results can then be validated via Monte Carlo simulations. A further avenue for validation is the comparison of values obtained at the laboratory scale with larger values obtained in the field.

Appendix A The Lognormal Distribution
The lognormal distribution has the form where k G = ⟨k⟩ exp(− 2 y ∕2) is the geometric mean, ⟨k⟩ the expected value and 2 y the variance of y = ln k . The permeability variance is 2 k = ⟨k⟩ 2 (e 2 y − 1) = k 2 G e 2 y (e 2 y − 1) , and the coefficient of variation is C vk = √ e 2 y − 1.

Appendix B The Gamma Distribution
The gamma distribution with a zero lower bound of permeability takes the form (Loáiciga et al. 2006) where and ( , > 0 ) are the shape and scale parameters, respectively, and Γ(⋅) is the gamma function; its mean and variance are given respectively by ⟨k⟩ = and 2 k = 2 ; the coefficient of variation is C vk = √ 1∕ . For = 1 , Eq. (B2) reduces to the exponential distribution.
An alternative expression of the gamma distribution (B2) as a function of its mean ⟨k⟩ and variance 2 k is

Appendix C Results for the Uniform Distribution
The continuous uniform distribution on support [k 1 , k 2 ] takes the form the mean and variance are given respectively by ⟨k⟩ = (k 1 + k 2 )∕2 and where F is the Forchheimer number, k is a dummy variable, 2 F 1 (⋅, ⋅; ⋅ ;⋅) is the hypergeometric function of parameters −1∕2 , c∕(2 − c) , 2∕(2 − c) and argument F and the transformation formulae (9.130) in Gradshteyn and Rhyzik (1994) can be used for analytic continuation if F > 1 , as 2 F 1 in Eq. (C8a-ba) converges in the entire unit circle. Comparing Eqs. (C8a-ba) and (12) leads to numerically determined values of k ep and ep for any c.

Author Contributions
All authors contributed to the study conception and design. AL and FZ have contributed equally to this work. All authors read and approved the final manuscript.