Horizontal landfill gas wells: geometry, physics of flow and connection with the atmosphere Damian Halvorsen(a) , Yana Nec(a)† and Greg Huculak(b) (a) Department of Mathematics, Thompson Rivers University, 805 TRU Way, Kamloops, BC, Canada (b) GNH Consulting Ltd., 20 50th st., Delta , BC, Canada † corresponding author: ynec@tru.ca November 15, 2018 Highlights † Landfill gas flow presents unique, but little explored regime and geometry. † This study identifies the significance of medium parameters over fluid properties. † In contrast to many applications gravity is shown to be non-negligible. † Well radius of influence and optimal well spacing are determined. † New accessible method is given to compute pressure, flux and radius of influence. Abstract Landfills form an important link in the chain of waste management. Sound engineering design and effective gas collection can only be attained via a thorough understanding of flow in different parts of the landfill. This study addresses gas flow induced by perforated horizontal wells commonly used for landfill gas collection. Planar flow through a porous medium of varying properties is solved numerically, yielding a twofold insight: quantitative deviations establish the validity and limitations of theoretical models, whilst qualitative comparison illustrates the impact of various flow parameters. It is found that an accurate modelling of geometry and medium resistance is essential in determination of a correct pressure distribution that can then be used to optimise gas collection. Gravity, an often overlooked effect, is shown to be non-negligible. The error source and bounds are given and explained for different landfill configurations. By contrast, temperature is of a lesser importance by an order of magnitude, implying that comprehensive spatial modelling of the exothermal reactions generating the gas as well as temperature dependent fluid properties such as viscosity are of a secondary interest. The well radius of influence, i.e. the maximal collection distance, is a fundamental design parameter impacting the interchange of fluid with the atmosphere. The numerical solution permits to obtain optimal values for effective collection and minimise flux between the landfill and atmosphere. Based upon a comprehensive numerical simulation an analytical approximation is suggested as an accessible design tool. The formula enables the computation of pressure profiles within the landfill and involves a minimal number of required input parameters. The surface flux and radius of influence can also be determined thereby. 1 Introduction Landfills give rise to a unique and little explored flow regime. The system is quite complex for several reasons. One, it combines flow fields of disparate characteristics: within the various landfill layers the flow is hindered 1 y S B ℓ s b A a P̃ a P Ã b B̃ 1 2 ⋯ i ⋯ N −1 N Figure 1: Flow geometry schematic, side view: perforated pipe imbedded within laminae of distinct permeability: gravel pack PA, waste stratum AB and cover layer BS. Dashed lines mark contiguity surfaces. Solid line S is the surface. Subscripts P, A, B, S are used throughout to refer to quantities at these surfaces. Subscript X denotes the outermost surface, i.e. B in the absence of a cover and S in presence thereof. Dotted lines mark planes aligned with perforated pipe cross-sections (cross-section geometry shown in figure 2). Arrows mark general direction of flow. Dimensions not to scale. Reproduced for convenience from [16]. y x S B Figure 2: Flow geometry schematic, cross-section for various mixed domains: pipe (thin solid circle), gravel pack (thin grey circle), waste stratum (thick black curves) and cover layer (thick grey curves). Dashed lines mark surface planes in the presence (S) and absence (B) of a cover layer. Dimensions not to scale. 2 similarly to a porous or fractured medium, whereas within the collecting pipes it is unobstructed. The fluid comprises gas as well as leachate, a liquid. Although early on an attempt was made to model a multi-phase flow [3], leachate flow is more often analysed separately, cf. [1, 19] and references therein. Two, a single landfill cell comprises several layers of distinct permeability, some acting as a filter and not generating. At times multiple cells are involved [5]. Three, within the waste stratum the medium parameters are uncertain and highly variable, depending on the type of waste at each locality. This was confirmed in numerous measurements [22, 14, 10], studies of mechanical properties [18, 15] as well as empirical permeability models [2, 21], all emphasising to a degree that availability of field data on medium properties would pose a limitation on accurate predictions. Amongst the more computationally intensive studies are some using neural network optimisation [13], graph theory [11] and closed commercial solvers [8, 12]. Scarcely few studies have developed analytical models for landfill gas collection as a fluid flow field. Axisymmetric equations for extraction by vertical [23, 5] as well as horizontal [6, 16] wells have been solved for certain combinations of boundary conditions. In all foregoing works gravity is neglected. In fact to date no studies considering gas flow have taken gravity into account. The flow velocities within the landfill are small, potentially rendering gravity non-negligible in the balance of vertical momentum. Substantial evidence thereof was given in [16] by inclusion of a gravel pack – a thin non-generating filter layer around the well – and revelation that most of the imposed pressure gradient dissipated over that layer, when its permeability was lower than that of the waste.1 The purpose of the current contribution is to combine an analytical approach and numerical simulation in order to extend the theoretical solution in [16] and explore more realistic and useful models. In particular the study reveals a local approximation that will be relatively simple to use in the field to estimate flux through boundaries of a landfill of any shape both in design and functional phases. To attain the ends above the partial differential equations governing the landfill gas flow field are solved numerically. Simple cases are verified via exact solutions, including a completely new one in elliptic coordinates. The full numerical solution is used to analyse a gravity inclusive problem with no axial symmetry assumption. Figure 1 depicts the geometry in question: a horizontal pipe with perforated cross-sections is supported by a gravel pack within a waste stratum, where landfill gas is generated. The extraction might be dynamic or passive with sub-atmospheric and atmospheric pressure at the pipe outlet respectively. A simplifying assumption of axisymmetric geometry and flow resulted in a viable leading order solution that recreated several interesting control related phenomena encountered in the field [16]. It was shown that the flow within each perforated cross-section might be deemed planar. This conclusion is called upon here in order to quantify the soundness of the axial symmetry assumption along two distinct lines of inquiry – domain shape and inclusion of gravity – and establish their influence on the flow field. A realistic landfill cell cross-section is not circular. Whilst taking the gravel pack as circular is fairly accurate, the waste stratum’s width often significantly exceeds its depth. The assumption of axial symmetry entails two artefacts. One, the cover, if present, circumscribes the landfill instead of being a surface layer only. And two, the depth of the landfill, being the smaller of the two dimensions, uniquely determines the radius where a boundary condition can be prescribed. Thus the formal solution cannot capture the flow field faithfully in the lateral regions containing the main bulk of waste. To address these problems elliptic and rectangular domains are examined. An elliptic domain rectifies the radial limitation by introducing separate horizontal and vertical axes lengths, whilst retaining a certain affinity to the circle. A rectangle represents real landfills perhaps the most closely, allowing for distinct lateral length and depth parameters, as well as a cover that is a surface layer only. Figure 2 depicts the relation between all aforementioned domains. Hereinafter a radial solution will denote a setting with full axial symmetry, whereas mixed elliptic or rectangular geometry will refer to a circular gravel pack and elliptic / rectangular outer boundaries respectively (tangent ellipses for three laminae). An elliptic solution will appertain solely to the problem in elliptic coordinates. 1 That situation is particularly applicable to demolition waste, highly non-uniform in fragment size and loosely compacted. 3 The steady state equations governing the flow of gas through interconnected pore spaces are (cf. [7, 20]) ∇ ⋅ (ρu) = C, k u = − (∇p − ρg) in Ω, µ (1a) where p and ρ are the gas pressure and density, related through the state equation for ideal gas p = ρRT with R and T being the gas constant and temperature respectively. u and µ are the fluid velocity and viscosity, C is the generation rate per unit volume and k represents an effective medium permeability: k= 2 φ3 r? , 18τ (1 − φ)2 (1b) where 0 < φ < 1, τ and r? denote porosity, tortuosity and the radius of an equivalent spherical particle respectively [4, 17]. g is the gravity vector. The domain Ω has the pipe as its inner boundary (thin black circle in figure 2) and any one of the thick grey or black curves as the outer boundary. Thus there are either two (cover absent) or three (cover present) sub-domains with gas generated in the second one. The boundary condition on the pipe is p(∂ΩP ) = pP , a constant associated with the suction imposed. The boundary condition on the outer circumference is either prescribed pressure p(∂ΩX ) = pX (Dirichlet type), where pX is the hydrostatic pressure distribution in the presence of gravity or a constant in the absence thereof, or no flux u ⋅ n̂ ∣ = 0 (Neumann type), where n̂ is the outward normal. To effect a proper ∂ΩX comparison between elliptic and rectangular domains a condition of only one type is allowed on ∂ΩX . On curves of contiguity between sub-domains continuity of p and u is imposed. The permeability k is deemed an effective quantity. As pointed out in [21], spatially distributed measurements over prolonged periods of time would be essential for accurate predictions, whilst remaining quite impractical since waste degradation, settlement and changes in the landfill configuration would soon nullify the modelling. The waste anisotropy is an inherent property, generally unpredictable and virtually impossible to quantify. Therefore intricate modelling as a porous or fractured medium in one setting does not appear to be applicable elsewhere. By contrast, an effective permeability enables qualitative analysis for design as well as functionality monitoring purposes2 . For instance, modelling leachate flow is outside of this study’s scope, however where leachate is recycled to percolate through the landfill or known to drain poorly, the overall medium permeability can be reduced accordingly. To isolate the impact of distinct physical phenomena, the effects of geometry and gravity were studied with the permeability held constant for each sub-domain. A linear variation of permeability with depth due to consolidation and degradation was also investigated. 2 Numerical solution The FlexPDE solver [24] was used to obtain a finite element solution with an unstructured triangular, dynamically refined mesh to observe a prescribed relative error in p (non-dimensional) of ε = 10−5 . Solutions for ε = 10−3 were found to coincide at least up to 11 significant figures, confirming adequate convergence of the numerical scheme. The numerical simulation was first completed in two cases, where exact solutions were available: radial flow in polar coordinates [16] and separation of variables solution in elliptic coordinates (appendix B). The magnitude of the deviation between numerically obtained and exact solutions (normalised by patm ) was at most O (10−7 ). In addition to confirming script correctness and numerical scheme accuracy, these were useful for qualitative analysis in domains combining circular inner and elliptic outer boundaries. In particular the radial and elliptic solutions feasibly reconstruct the pressure profiles on and in the vicinity of vertical and horizontal rays respectively. It stands to reason the radial solution will be a good approximation within its region of validity, because the axial geometry properly reconstructs the domain centre. The solution in elliptic coordinates requires the focal length as an arbitrary system parameter, somewhat distorting the 2 The permeability k can be estimated via (1b) or from handbook values published for different materials, as obtained with basic measurements by Darcy’s law, the second formula in (1a). 4 pB /patm pP /patm p p pS /patm rA rB r pP /patm 2rB pP /patm rB rS r 2rS p pS /patm p pB /patm rA rA rB r pP /patm 3rB rA rB rS r 3rS Figure 3: Pressure distribution (non-dimensional) for Dirichlet boundary condition: comparison of the numerical solution in mixed elliptic domains (solid) with exact solutions in radial (dashed) and elliptic (dotted) coordinates on horizontal (black) and vertical (grey) rays. Domain: two laminae (left) and three laminae (right). Lateral extent: ℓX = 2rX (top) and ℓX = 3rX (bottom). No gravity. All other parameters listed in appendix A. well’s vicinity. Further mathematical detail is given in appendix B. A typical pressure field is compared in figure 3 for Dirichlet boundary conditions and four geometric settings: lateral extent of double and triple the landfill depth for two (cover absent) and three (cover present) sub-domains. The discrepancies between vertical and horizontal profiles render the significance of domain asymmetry obvious. Hereinafter this is investigated with respect to various physical parameters. The nominal values of all geometric and physical parameters are listed in appendix A and are typical to landfills of medium size, yielding robust results, i.e. deviations in these values entail no qualitative change in the solution. 2.1 Medium and fluid properties The medium permeability k holds the salient responsibility for the resistance to fluid flow. By (1b) it depends on three underlying parameters. Since infinitely many combinations thereof result in the same value of k, the tortuosity τ is chosen to illustrate all impact due to the simplicity of its mathematical function: a proportional modification of all sub-domains’ permeabilities if φ and r? remain fixed. The feasible range of τ is 100 ⪅ τ ⪅ 1000. Figure 4 compares radial pressure profiles with those obtained for a mixed elliptic domain of increasing horizontal axis ℓX . For the boundary condition p(∂ΩX ) = pX juxtaposition of profiles along vertical and horizontal rays shows that all curves corresponding to a fixed value of τ are in good agreement with the radial solution in rP < r < rX . For the lower value of τ the profiles conforming to different ℓX coincide over their common range. The higher value of τ gives a bound of 1-2% to the possible discrepancy of pressure values predicted. However, observe that this deviation corresponds to 1/4 − 1/2 of the pressure 5 1.01 1.14 p p 1 1 pP /patm rP rB 2rB pP /patm rP 3rB rB 2rB 3rB 2rS 3rS r r 1.02 1.21 pP /patm p p 1 1 rP rB rS 2rS pP /patm rP 3rS r rB rS r Figure 4: Pressure distribution (non-dimensional): impact of τ for two (top) and three (bottom) laminae. Boundary conditions: Dirichlet (left) and Neumann (right). Values of tortuosity: τ = 100 (thin curves), τ = 1000 (thick curves). Rays: horizontal (black) and vertical (grey), coinciding for the radial case. Lateral extent: ℓX = rX (dotted), ℓX = 2rX (dashed), ℓX = 3rX (solid). No gravity. All other parameters listed in appendix A. drop pX − pP (compare to the ordinate range), and thus in fact is not small. By contrast, the boundary condition of no flux results in a significant span of pressure values for low and high values of τ even relative to patm . Figure 5 shows the effect of temperature and lends itself to a similar qualitative inference, albeit with respective deviations an order of magnitude smaller. 2.2 Inclusion of gravity By (1a) gravity assists gas collection above the pipe centre and partly counteracts the imposed pressure gradient beneath it. This asymmetry is more significant if the flow is quite slow, occurring when ka < kb , i.e. most of the suction is dissipated over the relatively thin gravel layer. With the introduction of gravity the boundary conditions must be made consistent with the expected asymmetry and type of flow. For that purpose the symmetric Dirichlet condition on the outer boundary was amended to read (y = 0 at the pipe centre) g(rX − y) p(∂ΩX ) = pX exp ( ), (2a) RT gy ), (2b) RT i.e. conform to the hydrostatic distribution of pressure in an ideal gas. Because the pipe radius is small relative to the landfill depth, i.e. rP ≪ rX , the inner condition is for all practical purposes equivalent to p(∂ΩP ) = pP exp (− 6 1.01 1.001 1 pP /patm p p 1 rP rB 2rB pP /patm rP 3rB rB 3rB 2rS 3rS 1.03 p p 1.004 1 pP /patm 2rB r r rP rB rS 2rS 1 pP /patm rP 3rS r rB rS r Figure 5: Impact of temperature: T = 15o C (thin curves), T = 215o C (thick curves). Graphics and flow conditions as in figure 4. p(∂ΩP ) = pP . Figure 6 gives the pressure profiles along horizontal and vertical rays with and without gravity for the mixed elliptic domain with ℓX = rX . With a Dirichlet condition (left panels) using realistic values in (2) gives a measurable difference on the order of magnitude of 0.25kPa between the uppermost and lowermost surfaces of the landfill, or equivalently 5 − 10% if related to the surface to well head loss of 3.75kPa. The difference at the end of a horizontal ray is visible and as expected about 1/2 of the above. With a Neumann condition the surface to well head loss is ≈ 1kPa with a similar gravity induced difference, amounting to 20 − 25%. These observations are independent of the domain’s lateral extent, as confirmed for ℓX = 2rX and ℓX = 3rX (not shown). Figure 6 allows for qualitative inferences as well. Each panel contains six curves, one pair representing the pressure profiles along the vertical ray below the pipe (grey, r < 0), one similar pair above it (grey, r > 0) and one pair depicting the variation along the horizontal ray (black, r > 0). On the top left panel the geometry is radial, therefore the vertical and horizontal profiles without gravity must coincide. Furthermore, the vertical profiles for r > 0 are indistinguishable, since the boundary conditions at the uppermost points of the landfill and pipe are identical with and without gravity. On the bottom left panel there also is no distinction between the two vertical profiles (they do not coincide with the horizontal one due to the presence of a cover layer). In light of the above two conclusions follow. One, the radial solution reproduces the pressure distribution faithfully within r < rX for domains of distorted shapes. And two, for the purpose of surface flux computation without a cover or with a partly permeable cover the gravity is not a material effect, since its impact is mainly felt in the lower part of the landfill. 7 0.97 p p 1.001 1 pP /patm −rB −rA rA pP /patm rB r −rA rA −rS −rB −rA rA rB r 0.97 p p 1.003 1 pP /patm −rB −rS −rB −rA rA pP /patm rB rS r rB rS r Figure 6: Pressure distribution (non-dimensional) with (thick curves) and without (thin curves) gravity on horizontal (black) and vertical (grey) rays. Domain: radial without cover (top) and mixed elliptic with cover (bottom). Boundary conditions: Dirichlet (left) and Neumann (right). Some curves coincide (see text), dominant line style shown. Lateral extent: ℓX = rX . All other parameters listed in appendix A. 2.3 Varying permeability Although it is accepted the permeability of the waste stratum decreases with depth due to consolidation and degradation of the medium, little is known on the diminution rate. Therefore a linear variation of k was studied for a rectangular domain. The permeability at y = rB (top of waste layer) was set to equal kB , the typical value used in all foregoing computations (equation (1b) and appendix A). The value at y = −rB (landfill bottom) was set to f kB with the factor 10−1 ⩽ f ⩽ 10−4 . Thus ⎞ ⎛ y − rB +1 . kAB (y) = kB (1 − f) 2rB ⎠ ⎝ (3) Figure 7 compares the corresponding pressure profiles with those obtained for the mean value kB /2. The asymmetry of the grey curves engendered by (3) is obvious next to the symmetric black ones. The differences are most noticeable near the bottom of the landfill, where the disparity between the very low permeability by (3) and the mean value is greatest. 3 Surface flux One of the quantities of interest in this problem is the flux of gas across the surface plane S or B. For a lateral extent ℓX ℓX ṁs = − ∫ −ℓX k ∂p ρ( + ρg)dx. µ ∂y 8 (4) pS /patm pP /patm p p pB /patm −rB pP /patm −rB rB 0 y pB /patm rB rS rB rS y p pS /patm p pP /patm 0 −rB pP /patm −rB rB 0 y 0 y Figure 7: Pressure (non-dimensional) along vertical sections of a rectangular domain with Dirichlet boundary condition: permeability linearly diminishing with depth in the waste stratum (grey, equation(3) with f = 0.001) and concomitant mean constant permeability (black). Section locations corresponding to pairs of black and grey curves from bottom to top on each panel: x = rP , x = 5rP , x = 10rP , x = 41 (rP + ℓX ), x = 12 (rP + ℓX ), x = 43 (rP + ℓX ). Domain: two laminae (left) and three laminae (right). Lateral extent: ℓX = rX (top) and ℓX = 2rX (bottom). No gravity. All other parameters listed in appendix A. In particular it is desired to determine which gas crosses the surface: landfill gas escapes when the vacuum is insufficient, with an obvious adverse environmental impact, or air is drawn into the landfill if the suction is too strong, decimating the population of anaerobic bacteria responsible for waste degradation. Thus it is important to find the optimal sub-atmospheric pressure to be imposed at the pipe boundary so as to minimise the flux. When the flux function is available, an illative quantity is the radius of influence, defined as the maximal lateral extent where the well collects gas. The first attempt to analyse the surface flux was done in [16], explaining mainly the variation along the well’s length. The radial solution was analytically extended to obtain the flux at the surface, to wit outside of the formal solution domain. Whilst successfully identifying the location of higher pressure gradients and flux, this approach permitted only an approximate prediction of the radius of influence: the flux as a function of horizontal distance from the pipe centre possessed no root, instead descending asymptotically to zero, thereby preventing a definitive estimate. Thence it was possible to define a cut-off point of a very small normal velocity that could be deemed zero for all practical purposes, and thus obtain the radius of influence. A conceptually disappointing feature in that study was the necessity to prescribe the cut-off point, a decision involving judgement and engineering experience. Moreover, the correctness of the threshold choice depends on the characteristic velocity in the landfill that varies with medium permeability. Here the full numerical simulation rectifies the aforementioned drawback by using a mixed rectangular domain, so that the flux surface of interest coincides with the uppermost boundary (figure 2). As a consequence the flux function possesses a simple root, entailing an unambiguous identification of the radius of 9 0.11 ṁs ×104 ṁs ×104 0.056 −0.61 −3rB −2rB −rB 0 x rB 2rB −0.97 3rB 0.09 −3rS −2rS −rS 0 x rS 2rS 3rS −3rS −2rS −rS 0 x rS 2rS 3rS ṁs ×104 ∆ṁs ×104 0 −1.11 −3rS −2rS −rS 0 x rS 2rS −0.16 3rS Figure 8: Surface flux with (dotted) and without (solid) gravity normalised by total generated mass in one longitudinal landfill segment C∆ . Lateral extent: ℓX = rX , 2rX , 3rX . Domain: mixed rectangular with no cover (top left), asymmetric cover (top right) and symmetric cover (bottom left). Bottom right panel shows respective differences between the last two. Dashed lines mark the change of sign, i.e. ṁs = 0, the intercepts with the curves giving the corresponding radii of influence. All other parameters listed in appendix A. influence. Figure 8 depicts typical flux variation as a function of horizontal distance for domains of aspect ratio 1 ⩽ ℓX /rX ⩽ 3. The impact of gravity is small and affects only the column of fluid immediately above the pipe centre. An interesting effect is that the radius of influence does not increase monotonically with the domain’s lateral extent. Figure 9 shows this dependence for a landfill without cover and with an asymmetric cover, the two realistic situations. In both cases for domains of low aspect ratio (ℓX /rX ⪅ 2 here, generally depending on the suction imposed on the pipe boundary) the radius of influence equals the lateral extent of the domain, i.e. the suction suffices to collect gas throughout the landfill. As the aspect ratio increases, the radius of influence reaches a maximum and then slowly diminishes, approaching an asymptotic value. The maximum corresponds to a critical lateral extent, where the side boundaries assist in gas collection by preventing outward flow. As the horizontal dimension increases further, the effect dissipates, giving the true radius of influence as an asymptotic value that depends on the sub-atmospheric pressure pP alone. 4 Local approximation via radial solution Hitherto evidence has been presented of the quantitative inadequacy of the radial solution as a global approximation: when the outer boundary of the domain is no longer circular, the radial solution cannot be expected to capture the pressure field accurately, since the boundary condition is given at a non-uniform distance from the pipe centre. This comprehension leads to the question whether the radial solution is useful in a local approximation. Figure 10 illustrates the difference between the numerical and radial solutions on a set of rays from the pipe centre to the outer boundary. On each ray the radial solution is calculated based 10 (∗) rI rmax (†) rmax rS rB rB rS r† r∗ 5rB 5rS ℓX Figure 9: Radius of influence versus lateral extent for mixed rectangular domains with a top cover (solid, (r∗ , rmax ) marks critical lateral extent and radius of influence) and no cover (dashed, (r† , rmax ) marks respective critical lengths). All other parameters listed in appendix A. (∗) (†) on the actual boundary distance. It is seen that for both mixed elliptic and rectangular domains the error is quite small, well under 1%. The maximal error is obtained on the longest ray, i.e. horizontal for an ellipse and corner ray for a rectangle, deep within the waste layer. In particular near the boundary itself the error is extremely small and thereon identically zero, implying that surface flux can be reasonably approximated by a simple numerical integration, not requiring a finite element solution. Note that (4) includes gravity. Omitting this term might significantly impact the result, since ρg is on ∂p the order of magnitude of 10, whilst ∣ ∣ is not very large. Formula (4) needs to be adjusted to the domain ∂y at hand. The radial solution reads [16] √ µ p= a ln r + b − RT Cr2 , (5a) 2k where a and b are constants depending on boundary conditions and the number of layers in the landfill. Closed form expressions for common situations are given in [16]. The concomitant radial velocity is ur = − ⎞ k dp k ⎛a µ =− − RT Cr , µ dr 2µp ⎝ r k ⎠ (5b) where k and C are the permeability and generation rate of the relevant landfill layer. To compute the flux for a symmetric rectangle, for instance, evaluate ℓX ṁs = −2 ∫ 0 pg ⎞ k p ⎛ ur cos ϑ + dx, µ RT ⎝ RT ⎠ 11 (6) ∆p ×103 ∆p ×103 3 2 1 2 1 0 −1 0 rP rA rB −1 2rB r rP rA rB r 3rB ∆p ×103 ∆p ×103 3 2 1 1 0 −1 2 0 rP rA rB r −1 2rB rP rA rB r 3rB Figure 10: Pressure difference (non-dimensional) between full numerical and ray-specific radial solutions for two laminae mixed elliptic (top) and rectangular (bottom) domains of lateral extent ℓX = 2rX (left) and ℓX = 3rX (right) on 31 rays evenly spaced between horizontal and vertical lines with Dirichlet boundary condition. All other parameters listed in appendix A. where ϑ is the angle measured clockwise from the vertical through the pipe centre. The interval [0, ℓX ] is discretised into √ n segments. The points then are xi = iℓX /n, i = {0, 1, . . . , n}. The radial distance to each 2 + x2 and cos ϑ = r /r . These define p by (5a) and u by (5b). Then the discrete version rX point is ri = i X i i i i of (6) is given by the trapezoidal method: ṁs ≈ − n−1 ℓX k ⎛ pi g ⎞ p0 g pn g p0 (u0 + ) + pn (un cos ϑn + ) + 2 ∑ pi (ui cos ϑi + ) , n µRT ⎝ RT RT RT ⎠ i=1 (7) the accuracy improving with increasing n. It is also possible to estimate the radius of influence by locating pg a point x∗ such that the expression ur cos ϑ + is positive at xi < x∗ and negative at xi+1 > x∗ . RT 5 Conclusion Landfill gas flow as a continuum flow field is by and large an unexplored application of fluid dynamics, combining several types of media and interacting with the atmosphere. Drawing upon a previously established inference of flow planarity within each cross-section aligned with pipe perforations, numerical solution of flow equations in porous sub-domains of distinct permeability and realistic geometry sheds light upon the strengths and limitations of the axisymmetric solution. This solution, with its one and only length dimension, gives a fair approximation in the vicinity of the vertical line between the pipe and surface, intriguingly even in the presence of gravity, but deteriorates with lateral distance. The deviation can amount to a few per cent relative to the pipe pressure with a Dirichlet type condition at the outer boundary, or reach as high 12 as 1/4 with a Neumann type condition, respectively conforming to partly permeable and sealed cover layers. Compared to the surface to pipe head loss, the same deviation ranges up to 1/2. These estimates include the impact of mixed geometry as well as gravity. The error magnitude depends mainly on permeability parameters of the porous media. The effect of temperature was found to be one order of magnitude less. Thence it follows that temperature dependent properties such as fluid viscosity will bear minor impact. This separation between the basic physical landfill parameters and gas characteristics is of the utmost importance for future studies. For instance, it is now evident that modelling heat transfer due to generation within the landfill is unlikely to yield results compelling from the vantage point of fluid dynamics or useful in engineering practices. On the other hand, better modelling of the medium properties might be of interest. In particular it implies a direct relation between the type of waste, its compactness, the thermodynamic work required to collect the generated gas and the landfill dimensions that would minimise the interchange of fluid with the atmosphere. One perhaps counterintuitive conclusion is that gravity is not as immediately dismissible an effect as is often surmised. The common assumption is based upon the relatively small depth of the landfill, further encouraged by the temptation to have one more axis of symmetry and simpler equations. However, the fluid velocity within the porous sub-domains is quite small, rendering the dynamic pressure comparable with potential energy. With a permeable cover the pressure is known on the surface and at the pipe, therefore in between the profile is virtually unaffected by gravity, whilst between the pipe and landfill bottom the hydrostatic pressure difference due to gravity alone is well measurable relative to the head loss due to medium resistance. By contrast, with a sealed cover the pressure is known at the pipe alone, and gravity has a twofold effect: hydrostatic as before, but also an intrinsic modification of the boundary condition: by ∂p = −ρg. The magnitude of the right hand side is what causes the (1a) zero flux on the surface implies ∂y significant disparity between the full and axisymmetric solutions. A crucial landfill design parameter is the radius of influence, the maximal distance where the well would effectively collect gas. Former efforts to estimate this quantity met with only partial success. The numerical solution for the first time attained an unambiguous identification of the radius of influence as well as an interesting feature to be useful in engineering design: there exists a critical lateral extent that maximises the radius of influence for a given suction strength. Lateral placement of adjacent landfill cells in accordance with this threshold should improve the collection efficiency. A vertical placement should be guided by the asymptotic value, beyond which no collection occurs, so that the wells do not work against one another. Then it is possible to find the optimal suction, so as to minimise the surface flux (the function, whose root gives the radius of influence) and thus the exchange of fluid with the atmosphere. Another revelation of import is that the radial flow does not allow for a sound approximation of surface flux when extended beyond the formal solution domain. The error thus begotten stems from artificial head loss over the additional sub-domain the fluid must traverse before reaching the surface, unduly diminishing the pressure gradient at the surface and hence flux. The current analysis suggests that the radial pressure profile can be used as a local approximator based on the actual distance of the desired point from the pipe centre. Hence in the absence of a full numerical solution an estimate of the surface flux might be obtained by a simple discrete integration of the normal velocity predicted on a dense set of rays from the pipe to the landfill surface. The radius of influence is estimated as the root of the integrand. This computation is an accessible tool useful in design as well as functional phases. It was shown that the impact of gravity as well as anisotropy in waste permeability is smallest near the surface, whereby application of this method in practice should entail a reasonable approximation for all quantities relevant to the interaction with the atmosphere. Appendix A. Base parameters Table 1 lists the nominal set of parameters used in computations throughout unless noted specifically in pertinent figure captions. Negative pressure values are relative to the atmosphere. 13 parameter pipe radius tortuosity temperature pressure on P pressure on B pressure on S generation rate segment generation CH4 molar fraction O2 molar fraction CO2 molar fraction parameter symbol rP τ T pP pB pS Cb C∆ xCH4 xO2 xCO2 gravel lamina PA 1m 0.6 0.025m depth φ r? value 0.0762 m (3 in) 100 15o C −3.75 kPa −0.125 kPa pbar 0.004 kg/(m3 hr) 12.6 kg/(m s) 0.5 0.01 0.4 waste lamina AB 8m 0.4 0.05m cover lamina BS 3m 0.7 0.005m Table 1: Parameters common to all examples solved numerically: global (top) and lamina specific (bottom) values (courtesy of GNH Consulting Ltd.) Negative pressure values are relative to pbar = 1atm. Appendix B. B.1 Solution in elliptic coordinates Single sub-domain Confocal ellipses and hyperbolae form an orthogonal coordinate system known as the elliptic coordinates. The relation to the Cartesian system is x = f cosh ξ cos θ, y = f sinh ξ sin θ, (8a) where f is the focal length, common to all ellipses and hyperbolae and deemed fixed. ξ > 0 is a non-dimensional curvilinear coordinate running along the hyperbolae and 0 ⩽ θ < 2π is the elliptic arc length, see figure 11. For the numerical solution in this coordinate system the focal length was taken as f = 0.97rP in order to have ξ > 0 on all relevant ellipses. The scale factors for the basis vectors hξ and hθ by (dx)2 + (dy)2 = h2ξ (dξ)2 + h2θ (dθ)2 are hξ = hθ = f √ cosh2 ξ − cos2 θ , (8b) (8c) and being equal, will be hereinafter renamed as h = hξ = hθ . In any orthogonal system {ê1 , ê2 , ê3 } with respective scale factors h1 , h2 , h3 the gradient of a scalar quantity F is given by ∇F = 1 ∂F 1 ∂F 1 ∂F ê1 + ê2 + ê3 , h1 ∂e1 h2 ∂e2 h3 ∂e3 (9a) and the divergence of a vector F = (F1 , F2 , F3 ) by ∇⋅F= ⎞ ∂ ∂ 1 ⎛ ∂ (h2 h3 F1 ) + (h1 h3 F2 ) + (h1 h2 F3 ) . h1 h2 h3 ⎝ ∂e1 ∂e2 ∂e3 ⎠ 14 (9b) Figure 11: Elliptic coordinate system. With no gravity the flow equation (1a) within a single sub-domain becomes ⎛ ∂2 ∂2 ⎞ + 2 p2 = −4f̃( cosh(2ξ) − cos(2θ)), 2 ∂θ ⎠ ⎝ ∂ξ f̃ = f2 µRT C. 4k (10) Due to the dependence of the right hand side on θ and in contrast to the axisymmetric geometry, it is ∂ ( ⋅ ) = 0, when C ≠ 0. For this reason the boundary impossible to obtain a normal flow, i.e. one where ∂θ conditions used in the solution are modified slightly to accommodate proper matching between sub-domains. Without loss of generality p2 = −f̃( cosh(2ξ) + P (ξ) cos(2θ) + Q(ξ, θ)), (11a) and thence P ′′ − 4P = −4, ⎛ ∂2 ∂2 ⎞ + Q = 0. ⎝ ∂ξ 2 ∂θ2 ⎠ (11b) Solving for P (ξ) is immediate. Seeking a solution by the method of separation of variables for Q(ξ, θ) under the periodicity p(ξ, θ) = p(ξ, θ + 2π) and reflexion conditions implied by the two symmetry axes of the ellipses p(θ) = p(−θ), p( π2 + θ) = p( π2 − θ), yields ∞ p2 = a0 ξ + b0 − f̃( cosh(2ξ) + cos(2θ)) + ∑ (an sinh(2nξ) + bn cosh(2nξ)) cos(2nθ), (12a) n=1 where an , bn , n ⩾ 0, are constants. Any boundary conditions that do not involve the higher harmonics cos(2nθ), n ⩾ 2, will result in p2 = a0 ξ + b0 − f̃ cosh(2ξ) + (a1 sinh(2ξ) + b1 cosh(2ξ) − f̃ ) cos(2θ). (12b) Below the parameters k and f̃ will bear the subscripts ( ⋅ )a , ( ⋅ )b and ( ⋅ )s , and the coefficients an , bn the superscripts ( ⋅ )(a) , ( ⋅ )(b) and ( ⋅ )(s) respectively for sub-domains ξP ⩽ ξ ⩽ ξA , ξA ⩽ ξ ⩽ ξB and ξB ⩽ ξ ⩽ ξS . In particular f̃a = f̃s = 0. 15 B.2 Boundary conditions The domain is [ξP , ξX ] × [0, 2π), where ξX = ξB (two sub-domains, cover absent) or ξX = ξS (three subdomains, cover present). Some conditions are common to all cases: pressure at the pipe √ p(ξP , θ) = p2P + θP cos(2θ) (13a) with pP sub-atmospheric and θP constant; and continuity on the contiguity curve ξ = ξA − + p(ξA , θ) = p(ξA , θ), − + u(ξA , θ) = u(ξA , θ). (13b) On the outermost circumference of the domain either pressure or zero normal flux conditions are prescribed: √ p(ξX , θ) = p2X + θX cos(2θ) or u ⋅ n̂ = 0, (13c) where pX ∈ {pB , pS } and similarly θX ∈ {θB , θS }. If a third sub-domain is present, continuity is imposed on ξ = ξB identically to (13b). The equations for the coefficients an , bn , n = 0, 1 as well as θP and θX , are linear. Solution existence is determined by the rank of this linear system as compared to the number of unknown coefficients – four for each sub-domain – and the degrees of freedom θP and θX . The solutions are given below. B.3 Cover absent: two sub-domains Case p2 (ξB , θ) = p2B + θB cos(2θ) ∶ (14a) (a) −1 a ( 0(b) ) = (kb (ξA − ξP ) + ka (ξB − ξA )) × a0 {(p2B − p2P + f̃b ( cosh(2ξB ) − cosh(2ξA ))) ( (a) b0 ξ −ξ kb ) + 2kb f̃b sinh(2ξA ) ( A B )} , ka ξA − ξP = p2P − a0 ξP , (a) b0 = p2B − a0 ξB + f̃b cosh(2ξB ), (b) (b) −1 cosh(2ξB ) −kb /ka cosh(2ξP ) θP a1 )( )+ ) = ( sinh(2(ξP − ξB ))) {( ( (b) −k /k sinh(2ξ ) sinh(2ξP ) θB a b B b1 (a) f̃b ⎫ ⎛ kb /ka cosh(2ξP )( cosh(2(ξB − ξA )) − 1) ⎞⎪ ⎪ ⎬, ⎝ sinh(2ξP ) − sinh(2ξB ) cosh(2(ξA − ξP ))⎠⎪ ⎪ ⎭ −1 θ k /k cosh(2ξB ) − cosh(2ξP ) a1 ) ( P) + ) = ( sinh(2(ξP − ξB ))) {( a b ( (a) − sinh(2ξ ) k /k sinh(2ξ ) θB B b a P b1 (b) ⎫ ⎛cosh(2ξB ) cosh(2(ξA − ξP )) − cosh(2ξP )⎞⎪ ⎪ ⎬, ⎝ kb /ka sinh(2ξP )(1 − cosh(2(ξB − ξA ))) ⎠⎪ ⎪ ⎭ ⎧ ⎪ ka kb ka sinh(2(ξB − ξA )) ⎪ + θB = −f̃b ⎨1 + cosh(2(ξB − ξA )) + × θP ⎪ kb sinh(2(ξA − ξP )) k − k k a b a − kb ⎪ ⎩ ⎫ ⎪ sinh(2(ξP − ξB )) + sinh(2(ξB − ξA )) cosh(2(ξA − ξP )) ⎪ ⎬. ⎪ sinh(2(ξA − ξP )) ⎪ ⎭ f̃b 16 Case u ⋅ n̂ = 0 ∶ (14b) (a) a0 = 2f̃b kb ( sinh(2ξB ) − sinh(2ξA )), ka a0 = 2f̃b sinh(2ξB ), (b) (a) b0 = p2P − a0 ξP , (a) b0 = (a0 − a0 )ξA + b0 + f̃b cosh(2ξA ), (b) (a) (b) (a) (b) f̃b a1 − sinh(2ξB ) )= ( (b) ( ), cosh(2(ξB − ξA )) cosh(2ξB ) b1 (a) (b) kb a1 a1 sinh(2ξA ) ( (a) )= {( (b) )} , ) + f̃b ( − cosh(2ξA ) ka b1 b1 θP = f̃b B.4 ⎞ kb ⎛ cosh(2(ξB − ξP )) − cosh(2(ξA − ξP )) . ka ⎝ cosh(2(ξB − ξA )) ⎠ Cover present: three sub-domains Case p2 (ξS , θ) = p2S + θS cos(2θ) ∶ ⎧ ⎪ ⎛ ⎪ (a) a0 = ⎨p2S − p2P + f̃b cosh(2ξB ) − cosh(2ξA ) − 2(ξB − ξA ) sinh(2ξA ) + 2(ξS − ξB )× ⎪ ⎝ ⎪ ⎩ ⎫ ⎫ ⎪ ⎪ ⎪ ⎧ kb ⎞⎪ ka ka ⎪ ⎪ ( sinh(2ξB ) − sinh(2ξA )) ⎬/⎨ξA − ξP + (ξB − ξA ) + (ξS − ξB )⎬, ⎪ ⎪ ⎪ ks ⎠⎪ ⎪ kb ks ⎪ ⎭ ⎭ ⎩ k a (a) (b) a0 = a + 2f̃b sinh(2ξA ), kb 0 a0 = (s) kb ka (a) a0 − 2f̃b ( sinh(2ξB ) − sinh(2ξA )), ks ks (a) b0 = p2P − a0 ξP , (a) b0 = b0 + ξA (a0 − a0 ) + f̃b cosh(2ξA ), (b) (a) (a) (b) b0 = p2S − a0 ξS , (s) (s) − cosh(2ξA )/ka a ), ( 1(s) ) = kb f̃b tanh(ξB − ξA ) ( cosh(2ξB )/ks a1 (a) sinh(2ξA )/ka b ), ( 1(s) ) = kb f̃b tanh(ξB − ξA ) ( − sinh(2ξB )/ks b1 (a) a1 − cosh(2ξA ) − sinh(2ξA ) ( (b) ) = f̃b {tanh(ξB − ξA ) ( )+( )} , sinh(2ξA ) cosh(2ξA ) b1 (b) k /k sinh(2(ξA − ξP )) θ ( P ) = f̃b tanh(ξB − ξA ) ( b a ). θS kb /ks sinh(2(ξS − ξB )) The case u ⋅ n̂ = 0 has no solution of the form (12b). 17 (15a) B.5 Tangential flow components The boundary conditions at ξP and ξB were allowed a tangential component to obtain the analytical solution of the form (12b). One must show that the tangential component is sufficiently small throughout the domain. Moreover, by (14a) for two sub-domains and pressure boundary condition a degree of freedom remains in the choice of θP or θB . Therefore it is desirable to minimise the resulting azimuthal velocity v or a related 1 quantity. Both these ends are attained by integrating either v 2 or the dynamical pressure ρv 2 . Note that 2 k 1 ∂p the integral thereof vanishes by v = − µ h ∂θ ξB 2π 2 ∫ ∫ vh dθdξ = 0. ξP Thus start with 0 ⎛ ∂p ⎞ 1 Iv = 2 ∫ ∫ k 2 dθdξ. µ ⎝ ∂θ ⎠ 2 ξB 2π ξP (16a) 0 The linear dependence in (14a) between θP and θB can be written for convenience as θB = aθ − bθ θP with aθ and bθ constants. This allows to write the pressure as p2 = A0 (ξ) + cos(2θ)(A1 (ξ)θP + B1 (ξ)). The functions A0 (ξ), A1 (ξ) and B1 (ξ) are continuous, but not differentiable at ξA , just like p. With this Iv can be shown to equal Iv = ξB π/2 2 2 ∫ k (A1 θP + B1 ) ∫ cos(2θ)× µ2 ξP (16b) 0 ⎧ ⎫ ⎪ ⎞ ⎛ ⎞⎪ ⎪ ⎪ ⎛ ⎨ ln A0 + A1 θP + B1 − 2(A1 θP + B1 ) sin2 θ − ln A0 − A1 θP − B1 + 2(A1 θP + B1 ) sin2 θ ⎬dθdξ. ⎪ ⎝ ⎠ ⎝ ⎠⎪ ⎪ ⎪ ⎩ ⎭ With the identity [9] √ π/2 1− 1+a 2 √ ln (1 + a sin x) cos(2x)dx = π ∫ 1 + 1+a 0 defining a1 = 2(A1 θP + B1 )/(A0 + A1 θP + B1 ) and a2 = 2(A1 θP + B1 )/(A0 − A1 θP − B1 ) leads to ⎛ 4π Iv = 2 ∫ k 2 A0 − µ ⎝ ξB √ ξP Differentiating with respect to θP , 2⎞ A20 − (A1 θP + B1 ) dξ. ⎠ ξB A1 (A1 θP + B1 )dξ dIv 4π = 2 k2 √ = 0. dθP µ ∫ 2 2 A0 − (A1 θP + B1 ) (16c) (17a) ξP If an extremum exists, it must be a minimum since ξB d2 Iv 4π = 2 2 dθP µ ∫ k 2 (A0 A1 ) dξ 2 2 3/2 2 − (A θ + B ) ) (A 1 P 1 0 ξP 18 > 0. (17b) RR R RR RR RR RR RR A1 θP + B1 RRRR RR A1 RR RR B1 RR RRR ≪ 1 and With the typical values in appendix A the ratios RRRR RRRR, RRRR RRRR are small, thus supposing RRRR A0 RRR RRR RRR A0 RRR RRR A0 RRR expanding in Taylor series gives ξB θmin ≈ − ∫ ξP ξB 2 2 A1 2 A1 B1 k dξ / dξ, k ∫ A20 A20 ξP (17c) which computation shows to be very small. The exact numerically obtained minimum is very close to this approximate value and indeed ∣θP ∣ ≪ 1. I/Imin 1.018 1 0 5e-08 1e-07 1.5e-07 θP Figure 12: Minimisation of square azimuthal velocity (Iv , grey) and dynamical pressure (Iq , black) to obtain an optimal value of θP for two laminae and pressure boundary condition. Respective approximate optima, (17c) and (19c), are marked in diamonds. Lateral extent: ℓX = rX . All other parameters listed in appendix A. Performing a similar analysis for the dynamical pressure is somewhat more complicated, as it is impossible to simplify the integral 1 1 2 ⎛ ∂p ⎞ 2 Iq = 2 dθdξ = 2 ∫ ∫ k p ∫ ∫ k (A1 θP + B1 ) cos(2θ) 2µ RT µ RT ⎝ ∂θ ⎠ 2 ξB 2π ξP 0 ξB 2π ξP 0 √ A0 + cos(2θ)(A1 θP + B1 ) dθdξ (18) that at first glance appears reducible to an elliptic integral, but for the compound (A1 θP + B1 ) changing sign. Nonetheless obtain ξB π ⎛ ⎞ dIq k 2 cos(2θ)A1 1 = 2 2A0 + 3 cos(2θ)(A1 θP + B1 ) dθdξ = 0. √ dθP µ RT ∫ ∫ ⎝ ⎠ A0 + cos(2θ)(A1 θP + B1 ) ξ P 0 19 (19a) If an extremum exists, it is a minimum, as ξB π d2 I q 1 = 2 dθP 2µ2 RT ∫ ∫ ξP 0 k 2 cos2 (2θ)A21 ⎛ ⎞ 4A0 + 3 cos(2θ)(A1 θP + B1 ) dθdξ > 0. ⎠ 3/2 ⎝ (A0 + cos(2θ)(A1 θP + B1 )) (19b) The inequality follows by the expression in the last set of parentheses equalling A0 + 3p2 and A0 > 0. RRR R A1 θP + B1 RRRR Expanding as before with RRRR RR ≪ 1 allows to complete the integration dθ, giving RRR RRR A0 R R ξB θmin ≈ − ∫ k ξP ξB A2 √ dξ / ∫ k 2 √ 1 dξ. A0 A0 2 A1 B1 (19c) ξP This approximate optimum is very close to the exact value obtained numerically, upholding the same conclusions as for the square tangential velocity. Figure 12 shows a typical variation Iv (θP ) and Iq (θP ) and the approximate minima based on the foregoing derivation. Both are small, close to each other and the respective exact minima. Moreover, the minimum is shallow, and for convenience in the numerical solution it is possible to adopt θP = 0, avoiding an iterative procedure. References [1] A.A. Al-Thani, R.P. Beaven, J.K. White. Modelling flow to leachate wells in landfills. Waste Management, 24(3): 217–276, 2004. [2] S.G. Arigala, T.T. Tsotsis, Y.C. Yortos, J.J. Kattapuram, I.A. Webster. Gas generation, transport and extraction in landfills. J. Environ. Eng. ASCE, 121(1): 33–44, 1995. [3] C. Berglund. Environmental hydraulics of sanitary landfills. J. Hydraulic Research, 36(1): 101–116, 1998. [4] P.C. Carman. Fluid flow through granular beds. Transactions of Institution of Chemical Engineers, 15: 150–167, 1937. [5] S-J. Feng, Q-T. Zheng, H-J. Xie. A gas flow model for layered landfills with vertical extraction wells. Waste Management, 66: 101–113, 2017. [6] S-J. Feng, Q-T. Zheng, H-J. Xie. A model for gas pressure in layered landfills with horizontal gas collection systems. Comput. Geotech., 68: 117–127, 2015. [7] W.B. Fulks, R.B. Guenther, E.L. Roetman. Equations of motion and continuity for fluid flow in a porous medium. Acta Mechanica, 12: 121–129, 1971. [8] D.K. Fytanidis, E.A. Voudrias. Numerical simulation of landfill aeration using computational fluid dynamics. Waste Management, 34: 804–816, 2014. [9] I.S. Gradshteyn, I.M. Ryzhik. Table of integrals, series and products. 7th edition, edited by Alan Jeffrey and Daniel Zwillinger, Academic Press, Elsevier; identities 4.399, 2007. [10] M. Hrad, M. Huber-Humer. Performance and completion assessment of an in-situ aerated municipal solid waste landfill – final scientific documentation of an Austrian case study. Waste Management, 63: 397–409, 2017. 20 [11] J.D. Hyman, S. Karra, N. Makedonska, C.W. Gable, S.L. Painter, H.S. Viswanathan. dfnWorks: a discrete fracture network framework for modeling subsurface flow and transport. Comp. Geosci., 84: 10–19, 2015. [12] P. Jain, T.G. Townsend, T.M. Tolaymat. Steady-state design of horizontal systems for liquids addition at bioreactor landfills. Waste Management, 30: 2560–2569, 2010. [13] L. Liu, J. Ma, Q. Xue, J. Shao, Y. Chen, G. Zeng. The in situ aeration in an old landfill in China: multi-wells optimization and application. Waste Management, 76: 641–620, 2018. [14] H-J. Kim, H. Yoshida, T. Matsuto, Y. Tojo, T. Matsuo. Air and landfill gas movement through passive gas vents installed in closed landfills. Waste Management, 30: 465–472, 2010. [15] K. Münich, J. Bauer, K. Fricke. New measurement techniques for the monitoring of landfills. Land Contamination and Reclamation, 13(1): 49–59, 2005. [16] Y. Nec, G. Huculak. Landfill gas flow: collection by horizontal wells. Preprint [17] M.N. Panda, L.W. Lake. Estimation of single-phase permeability from parameters of particle-size distribution. American Association of Petroleum Geologists Bulletin, 78(7): 1028–1039, 1994. [18] M.K. Singh, J.S. Sharma, I.R. Fleming. A design chart for estimation of horizontal displacement in municipal landfills. Waste Management, 29: 1577–1587, 2009. [19] R. Slimani, L. Oxarango, B. Sbartai, A-J. Tinet, F. Olivier, D. Dias. Leachate flow around a well in MSW landfill: analysis of field tests using Richards model. Waste Management, 63: 122–130, 2017. [20] S. Whitaker. Flow in porous media I: a theoretical derivation of Darcy’s law. Transport in Porous Media, 1: 3–25, 1986. [21] G.M. Williams, R.S. Ward, D.J. Noy. Dynamics of landfill gas migration in unconsolidated sands. Waste Manage. Res., 17(5): 327–342, 1999. [22] N.D. Woodman, T.C. Rees-White, A.M. Stringfellow, R.P. Beaven, A.P. Hudson. Investigating the effect of compression on solute transport through degrading municipal solid waste. Waste Management, 34: 2196–2208, 2014. [23] L. Yu, F. Batlle, J. Carrera, A. Lloret. Gas flow to a vertical gas extraction well in deformable MSW landfills. Journal of Hazardous Materials, 168: 1404–1416, 2009. [24] FlexPDE 7, PDE Solutions Inc. http://www.pdesolutions.com 21