Solution of weakly compressible isothermal flow in
landfill gas collection networks
Y Nec1 ‡ and G Huculak2

1

Thompson Rivers University, Kamloops, British Columbia, Canada

E-mail: cranberryana@gmail.com
2
GNH Consulting, Delta, British Columbia, Canada
E-mail: greg@gnhconsulting.ca
Abstract. Pipe networks collecting gas in sanitary landfills operate under the
regime of a weakly compressible isothermal flow of ideal gas. The effect of
compressibility has been traditionally neglected in this application in favour of
simplicity, thereby creating a conceptual incongruity between the flow equations
and thermodynamic equation of state. Here the flow is solved by generalisation
of the classic Darcy-Weisbach equation for an incompressible steady flow in a
pipe to an ordinary differential equation, permitting continuous variation of
density, viscosity and related fluid parameters, as well as head loss or gain due
to gravity, in isothermal flow. The differential equation is solved analytically in
the case of ideal gas for a single edge in the network. Thereafter the solution is
used in an algorithm developed to construct the flow equations automatically
for a network characterised by an incidence matrix, and determine pressure
distribution, flow rates and all associated parameters therein.

Keywords : flow rate, head loss, Darcy-Weisbach coefficient, pipe network,
weakly compressible flow

‡ corresponding author

2

Flow solution in landfill networks
1. Background
1.1. The landfill problem

Landfill collection networks manifest flow conditions unique both within the gas industry
and the field of hydraulics. Most industrial applications involving gas pipeline flow are
characterised by a fully compressible regime, whilst hydraulic networks usually conduct
incompressible fluid, whence the variation of temperature is of scarce interest. The flow
in a landfill network falls somewhere in between: the fluid is weakly compressible and the
flow cannot be solved without determination of the temperature variable. Moreover, this
particular type of network presents the additional complication of multiple source points
with streams of distinct chemical composition mixing at each node. Traditionally this
flow was approximated as incompressible for simplicity, albeit the mainstay complexity
did not stem from the choice of a head loss model. The current contribution sets forth
the theory that is conceptually consistent from the vantage point of fluid mechanics
and thermodynamics, and accounts for properties peculiar to this application, such as
differences in elevation between the network nodes and early onset of turbulence.
1.2. Classic Darcy-Weisbach equation
Historically Darcy-Weisbach equation was developed for steady incompressible flow upon
the underpinning of dimensional analysis combined with empirical measurements to
estimate the head loss |∆p| in a horizontal pipe of diameter d and length ℓ:
ℓ
1
|∆p| = ρu2 f ,
2
d

(1)

wherein ρ denotes the fluid density, u is flow velocity, and f represents a dimensionless
friction coefficient as given by the Colebrook equation (Colebrook, 1939; White, 1999)
for a fully turbulent regime


2.51
ε
1
√ = −2 log10
+ √
(2)
3.7 Re f
f
and involving the relative surface roughness ε and Reynolds number Re, both based
on diameter d. Some remarks regarding (2) are in order. The salient moment of this
phenomenological approach is a realistic estimate of |∆p| owing to the inclusion of
interaction with a non-smooth surface through f , whilst Navier-Stokes flow equations
entail solely fluid properties. That f is a function of Re and ε can be obtained by
dimensional analysis, cf. Benedict (1980). Full turbulence in the current problem of
landfill networks is obtained somewhat earlier than the commonly concomitant with (2)
threshold of Re > 3000 (Moody, 1944), the reason being twofold: one, Moody’s diagram
is for flow in a single pipe with no geometric disturbances, whilst in the network the
pipes’ junctures engender turbulence; and two, the fluid comes into the network from
vertical wells, wherein flow rate measuring devices are installed, e.g. orifice plates,
rendering the flow turbulent even at very low Re numbers. Characteristic Re numbers

3

Flow solution in landfill networks

√
in this application are 103 − 105 . Equation (2) is implicit in f and must be solved
numerically or equivalently through Lambert W function (Corless et al., 1996).
Forthright application to flows with inherent compressibility is impossible, however.
The current contribution is inspired by the unique flow regime occurring in gas collection
networks in landfills. Simplistically depicted, the landfill comprises an underground
organic mass, wherein anaerobic biodegradation takes place; a covering substratum with
an imbedded system of collection pipes; and an extraction facility. The biodegradation
of waste is exothermal, and the amount of gas generated is rarely uniform within the
landfill. The gas is conducted through a series of wells into the system of pipes, the
flow induced by suction created at the extraction facility. The landfill is not absolutely
impermeable to ambient air, but when the landfill pressure is properly monitored, there
is no exchange between the landfill and atmosphere. Two circumstances decide the
thermal regime: one, the pipe network is situated within the substratum sufficiently
far both from the landfill cavity, where the chemical reactions occur, and from the
surface; and two, the pipes are not insulated, allowing dissipation of friction generated
heat. A thermal equilibrium thus ensues within the substratum, midmost layers
thereof maintaining a constant temperature, virtually unaffected by seasonal changes.
Obviously, where the pipes are above ground, isothermal conditions are unlikely to
occur. Henceforth the flow within the collection network is taken isothermal. The
thermodynamic properties of the fluid render the assumption of ideal gas tenable,
forthwith imparting flow compressibility, as is seen from the state equation
p = ρRT,

(3)

p denoting pressure, T – temperature and R – the gas constant, whereby pressure p and
density ρ must perforce vary proportionately since T and R are fixed.
Estimation of flow rates in the collection network is a basic step in landfill design,
whence a relation akin to (1) is of the essence. Surprisingly no adequate model is in use.
The sundry branches of the gas industry have developed empiric formulae galore, each
suitable to a particular flow regime and swiftly losing accuracy when applied elsewhere
(Menon, 2005). Some justification of the physical kernel of the aforesaid formulae can
be found in Coelho and Pinho (2007) and references therein. In addition to custom fit
exponents two features of import prevail: compressibility factor z < 1 and efficiency η,
or variants thereof. Introduction of the compressibility factor is equivalent to adoption
of a more general equation of state
p = zρRT,

z ⩽ 1,

(4)

however here z < 1 will be incongruent with the regime of weak compressibility. The
multiplicative efficiency factor η is a crude attempt to account for natural phenomena
pertaining to head loss that elude properly grounded, yet pragmatic modelling. By
contrast, Darcy friction factor f manifests noteworthy generality through its continuous
dependence on physical parameters of both fluid and conduit, namely ε and Re.

4

Flow solution in landfill networks

Changes in fluid density in the landfill collection network rarely are such as to render
the accuracy of (1) with a fixed value of ρ appalling. Characteristic Mach numbers range
below 0.1, therefore the flow is only weakly compressible. Nonetheless, application of
(1) is conceptually incorrect and entails significant error when pipe lengths are large.
In that light the current contribution generalises (1) for continuously varying density
ρ, yet preserves the valuable complexion of Darcy friction factor f . Affording accurate
estimates of head loss in landfill pipework is essential for determination of gas quantities
extracted, as will become apparent infra, and subsequent effectuation of sustainable
practices such as diversion of the gas for beneficial uses.
2. Continuous Darcy-Weisbach equation
Hereinafter integral mass and momentum equations are used to derive an ordinary
differential equation to supplant (1). The layout of landfill pipes follows local landscape,
therefore inclusion of gravity in the balance of momentum is imperative, as under low
flow conditions head loss or gain due to sloping terrain becomes non-negligible. Consider
a control volume of a straight pipe of length ℓ and radius r, inclined at angle α as in
figure 1. The flow is steady. Mass conservation equation reads (Batchelor, 1990)
I
ρu · ds = 0,
(5)
∂V

wherein ∂V is the surface of the control volume V , ρ is fluid density, u denotes velocity
vector and s – area vector with the normal directed outwards. Integration yields
ρin uin = ρout uout

(6)

with subscripts ( · )in and ( · )out referring to the respective cross-sections, and u = |u|
being the velocity magnitude. For consistence hereunder velocity and other quantities
inherently varying within the cross-section tacitly refer to equivalent uniform values.
Momentum conservation equation reads (Batchelor, 1990)
Z
Z
Z
u (ρu · ds) +
p ds =
ρ fbody dV + fsurf ,
(7)
∂V

∂V

V

wherein p is fluid pressure, fbody is the gravity vector, and the only surface force fsurf is
due to the shear stress τ within the boundary layer, via dimensional analysis given by a
friction coefficient cf and integration of dynamic pressure 21 ρu2 along the cylinder walls:
fsurf = πr

Z ℓ

cf ρu2 dx.

(8)

0

Applying (7) in the axial direction yields
g sin α
−ρin u2in + ρout u2out − pin + pout =
2
πr

Z

1
ρ dV −
r
V

Z ℓ
0

cf ρu2 dx.

(9)

5

Flow solution in landfill networks

in

α
α
out

Figure 1: Control volume schematic. Angle α is defined positive as shown, conforming to head gain due to gravity.
def

For incompressible flow, to wit density constant throughout ρin = ρout = ρ and implying
def
velocity similarly constant uin = uout = u, equation (9) becomes
!
Z
u2 ℓ
pout − pin = gℓ sin α −
cf dx ρ,
(10a)
r 0
recovering a gravity inclusive variant of (1)
pout − pin =

1 f
g sin α − u2
2 d

!

ρℓ

(10b)

upon setting d = 2r and
1
ℓ

Zℓ

1
cf dx = f.
4

(10c)

0

In this sense the classic Darcy-Weisbach friction coefficient f is a multiple of the average
friction coefficient, or inversely, thenceforth the friction coefficient is regarded as an
effective average value that up to a multiplicative constant locally equals the DarcyWeisbach coefficient f .
To continue, suppose a generic equation of state, incorporating both (3) and (4)
as particular cases (even if z < 1, the thermodynamic properties of the fluid are not
expected to vary insofar as to render z a function of x),
ρ = F(p, T ),

(11)

for isothermal flow giving together with (6)
ρin
F(pin , T )
uout
=
=
.
uin
ρout
F(pout , T )

(12)

6

Flow solution in landfill networks
Further combining with (6) and (9) gives
F(pin , T )
pout − pin
1
1−
−
=
2
F(pout , T )
ρin uin
ρin u2in

1
r

Z ℓ

g sin α
cf ρu dx −
πr2
0
2

Z

ρ dV
V

!

.

(13)

To facilitate a continuous variation of ρ along the control volume, map the finite length
ℓ to infinitesimally small δℓ, whereby ρin 7−→ ρ, uin 7−→ u, pin 7−→ p, pout 7−→ p + δp,
ℓ 7−→ δℓ in (10c) and
F(pout , T ) ∼ F(pin , T ) +

dF
dp

in



δp + O δp2 .

(14)

Then (13) becomes
1 dF
1
δp
− 2
F dp ρu

!

=

!
f
g sin α
δℓ
−
4r
u2

(15)

upon completion of integrals and omission of higher order terms. At the limit the
following ordinary differential equation ensues
!,
!
f
1 dF
dp
g sin α
1
δp
=
=
−
−
.
(16)
lim
δp−→0 δℓ
dℓ
4r
u2
F dp ρu2
δℓ−→0

Note that (16) was derived through integral mass and momentum conservation with a
generic equation of state and without making any assumptions regarding the variation of
the friction factor f . A less generic version can be found in classic literature, derived for
an ideal gas, horizontal flow (α = 0) and f gratuitously assumed constant throughout
(Benedict, 1983; Ward-Smith, 1980). The friction factor f depends on the Reynolds
number, thus introducing a dependence on the fluid viscosity, which is known to vary
with temperature. Hence implicitly (16) involves all properties sufficient to determine
the thermodynamic state of a substance and will require a specific equation of state for
solution.
The function F is readily identified as F = p/(zRT ), and (16) simplifies to
!,
!
dp
f
1
g sin α
1
=
−
−
.
(17)
dℓ
4r
u2
p ρu2
If the equivalence (10c) is premised, i.e. the constraint of Darcy-Weisbach relation
conforming solely to incompressible flow is relaxed ad hoc, equation (17) can be obtained
from the differential form of the mass and momentum conservation equations:
d(ρu) = 0
u du +

f
dp
− g sin α dℓ + u2 dℓ = 0.
ρ
4r

(18a)
(18b)

7

Flow solution in landfill networks

Before integration of (17) is attempted, it is necessary to determine how the friction
coefficient f varies. By (2) f is a function of ε, fixed for a given pipe, and Re, with the
expression of mass flow rate
ṁ = πr2 ρu,
(19)
to be written as
Re =

4ṁ
.
πµd

(20)

The mass rate ṁ is constant by (6). The viscosity µ depends on temperature.§ Since
the flow is isothermal, µ is constant, whereby so are Re and f .
Recollecting (3) and using (19), for α 6= 0 equation (17) is forthwith cast into the
more elegant form
g sin α p(p2 − p2eq )
dp
,
(21a)
=
dℓ
RT
p2 − p2∞
wherein

p2eq =

ṁRT
2πr2

!2

f
,
gr sin α

√
ṁ RT
p∞ =
.
πr2

(21b)

For α > 0, i.e. downward slope, the pressure peq entails a perfect balance of head gain
due to gravity and loss due to friction. For α < 0, i.e. upward slope, (21) possesses
2
no fixed
√ points as peq < 0. The singular point p∞ corresponds to the subsonic velocity
u = RT k and has been obtained heretofore for horizontal geometry (Benedict, 1983).
In the current system this singular point is never approached, the flow being weakly
compressible with Mach number below 0.1 and 1.2 < γ < 1.25.
2.1. Solutions
Both forms (17) and (21) are separable and easily integrated. A simple familiar case is
that of ṁ = 0, whereby the integration is immediate, recovering the result for hydrostatic
pressure change
gℓ sin α
pout
.
(22)
=
ln
pin
RT
For α = 0, i.e. horizontal flow, integration of (17) over a pipe of length ℓ gives

4r pout
2π 2 r5  2
2
ln
−
p
−
p
in = ℓ,
f
pin
f ṁ2 RT out

(23a)

 3/2
µ
T
To + s
§ Common models are Sutherland’s formula
=
(Crane, 1982) or exponential law
µo
To
T +s
 n
µ
T
=
(White, 1999) with µo , To , s, n experimentally determined constants.
µo
To
√
k For an adiabatic flow of ideal gas Mach number of unity is reached at u = γRT with γ > 1 being
the ratio of heat capacities cp /cv .

Flow solution in landfill networks

8

a transcendental equation that cannot be solved as is for an explicit expression akin to
(10). For α 6= 0 form (21) is preferable and upon integration yields ¶
!
p2out − p2eq
p2∞ pout 1
gℓ sin α
p2∞
ln
=
+
,
(23b)
1 − 2 ln 2
2
2
peq
pin
2
peq
pin − peq
RT
similarly insoluble for pout in closed form without further comparison of the various
quantities’ magnitude. Equations (23) pose no difficulty when solved numerically by
methods commonly used for non-linear equations.
To exemplify the differences effected by the continuous variation of density,
equations (23) and (10) were solved for two situations:

p
out
< 1 for horizontal flow
† downstream prediction: pin given, compute pout
p
 p in

in
† upstream prediction: pout given, compute pin
> 1 for horizontal flow .
pout
Albeit the former is the more commonly encountered case, the upstream prediction is
the inherently applicable case for the landfill networks, since the pressure maintained
at the extraction facility dictates the conditions up the stream. The magnitude of the
difference in head loss as predicted by the continuous model (23) versus the classic (10)
varies with the flow parameters. For a typical example see figure 2. The conceptual
disparity between (23) and (10) is summarily discernible: by (10) the density is constant
and pressure decreases linearly down the stream, whereas (23) imbeds a continuous
change of density through (3) and entails a non-linear diminution of pressure. When the
density is deemed constant and taken as the entry value, (10) overestimates the pressure
downstream (figure 2, left panel). Similarly, with an exit value (10) will overestimate the
pressure up the stream (figure 2, right panel). As expected, for small ℓ density changes
scarcely evolve, rendering the foregoing difference negligible. Over long distances the
changes due to compressibility are more pronounced. The continuous model (23) should
therefore be considered the correct one in the case of weakly compressible flow, where
usage of (10) is conceptually erroneous, albeit as a discrete approximation giving a close
estimate for sufficiently small ℓ.
3. Application to landfill gas collection network
The continuous Darcy-Weisbach equation (23) was used to solve the flow in a landfill
collection network. The foremost property that distinguishes this network from a
common hydraulic network is that each node serves as a source point of fluid with
distinct thermodynamic properties. The gas in each well flowing from the landfill cavity
upwards has its own composition that depends on the decay processes in that vicinity.
The stream of each well mixes with the gas incoming into the node by pipes collecting
gas from other parts of the network.
¶ The equilibrium solution pin = pout = peq for α > 0 does not require integration.

9

−6

137

p − patm [kPa]

p − patm [kPa]

Flow solution in landfill networks

89

−30

42

−54

0

ℓ [m]

5000

−6

0

ℓ [m]

Figure 2: Comparison of head loss prediction with correct inclusion of compressibility effect (equation (23), thick
black curve) versus classic Darcy-Weisbach formula neglecting density changes (equation (10), thin grey line). Left
panel: downstream prediction (pin given, compute pout ) results in overestimation of outlet pressure by (10) based
on fixed entry density. Right panel: upstream prediction (pout given, compute pin ) results in overestimation of inlet
pressure by (10) based on fixed exit density. Parameters used: relative roughness ε = 10−5 (absolute roughness of
−6 m for plastic and drawn tubing, consult Crane (1982, A-23), and diameter as below), T = 13o C, molar fractions
1.5×10
o
n
o
n
xCH4 , xO2 , xCO2 , xN2

=

0.5, 0.01, 0.45, 0.04 , α = arctan(0.05), ṁ = 0.2kg/s, d = 0.15m, µ = 1.4 × 10−5 Pa · s

(corresponding Re = 120, 000), reference pressure p = 6kPa below atmospheric (pin for left panel and pout for right
panel). All numeric values given up to two significant digits.

A generic solution to such a system is conceptually impossible without a given
network configuration. Hereunder the problem is posed for a general network with the
aid of fundamentals of graph theory, and a solution algorithm is developed. Thereafter
a solution is presented for a specific example of a network of sufficient difficulty
to be representative with the purpose to establish that implementation of (23) was
computationally feasible on par with the classic model (10).
The generic geometry of the system is as follows. The flow network is a graph with
n nodes corresponding to extraction wells and e edges representing collection pipes. An
example of typical complexity is given in figure 3. One of the nodes is designated as the
flare, where sub-atmospheric pressure is imposed, inducing flow from the wells that act as
sources for the network. Two parameters control the collection: suction strength at the
flare as dictated by considerations of efficiency, and a maximal permissible pressure (subatmospheric) throughout the system as set by environmental and safety regulations+
(Conestoga-Rovers, 2010). The well flow rates are deemed known input.∗
To determine the flow rates in the collection pipes one is required to solve a system
+

Excessive pressure bespeaks insufficient gas extraction and is bound to beget leaks through the landfill
cavity boundary into the atmosphere.
∗ During landfill design phase well flow is estimated based on the amount and nature of refuse to
be buried at that locality. When fully operational, each well’s flow rate is measured. Certain field
conditions necessitate additional care to achieve adequate accuracy (Nec and Huculak, 2015).

5000

10

Flow solution in landfill networks
15
14

14
21
23

22

23

17
16

22

24

21
27

13
25

18
16

18

19

17

13

20
19

20

26
15

9

12

28
9

8
10

8

11

12

10

11

7
7
6
1

6

1

5
2

5

2

4
3

4
3

Figure 3: Example graph: vertices correspond to wells (numbered in diamonds), edges correspond to pipes (numbered in
circles). Flare is located at node #1, where sub-atmospheric pressure is imposed. The general collection progress is from
the remote wells toward the flare.

of non-linear equations. A network flow is governed by two principles: conservation of
mass at the nodes and continuity of pressure over any closed path. A graph of n nodes
and e edges contains e − n + 1 fundamental cycles (Godsil and Royle, 2001), i.e. number
of independent closed paths, conforming to the number of pressure continuity equations
required. Perforce the number of independent mass conservation equations will be n− 1.
Hence define Jiin and Jiout respectively as the sets of incoming and outgoing pipes at the
ith node. With ṁj denoting the flow rate through the j th pipe, mass conservation
equations read
X
X
ṁj + ṁ(i) =
ṁj , i = 1, . . . , n − 1,
(24a)
j∈Jiin

j∈Jiout

wherein ṁ(i) is the mass flow rate input by the ith well. Construct the network’s incidence
matrix as follows (Godsil and Royle, 2001). The matrix is of size n × e. Each edge is
assigned an arbitrary flow direction. In each column, corresponding to an edge, entry
ij is set to equal (−1) if the flow in pipe j proceeds from node i, unity if the flow enters
node i, and zero otherwise. Each row is then used to determine the set Jiin ∪ Jiout and

Flow solution in landfill networks

11

write (24a) in matrix form.
Define Jkc as the set of pipes forming the k th fundamental cycle. With ∆pj being the
pressure drop over the j th pipe (sign determined by flow direction), pressure continuity
equations read
X
∆pj = 0, k = 1, . . . , e − n + 1.
(24b)
j∈Jkc

The fundamental cycles are obtained by finding the null space of the incidence matrix
(Godsil and Royle, 2001). The dimension of that space gives the number of required
pressure continuity equations. The vectors spanning the space allow to write (24b) in
matrix form. Together (24a) and (24b) give a system of e equations in e unknown mass

j=e
flow rates ṁj j=1 , with (24a) being linear, whereas (24b) are non-linear by (23). It
is possible to construct the incidence matrix conforming to an acyclic sub-graph of the
network, sized n × (n − 1), such that the inverse of a square block (n − 1) × (n − 1)
thereof, computed only once, enables the determination of n − 1 out of e unknowns via
(24a). Thence only the non-linear system (24b) in e − n + 1 unknowns requires iterative
solution.
Solution of (24) is somewhat convoluted for the particular application of a landfill
network, because the gas differs in composition between the different wells. Simple
causality implies that properties of the mixture can be computed only in the downstream
direction, whereas head loss must be computed starting with the prescribed flare pressure
and proceeding upstream. Note that, identically to (10), both (23a) and (23b) accord the
determination of pout (pin ) (downstream prediction) and pin (pout ) (upstream prediction)
with equal ease. This attribute proves of import in the solution of (24b), since in a
network of that ilk well flow rates might fluctuate significantly over prolonged periods
of time as a cause of natural changes in chemical reactions within the landfill, bringing
about flow reversal in particular parts of the network, an occurrence all the more likely
when the terrain manifests elevation changes. In light of the above, when equation (23)
supplants (10) in the computation of pressure, the impact on the overall computational
complexity of (24) is peripheral.
3.1. Mixing
A further peculiarity of the landfill network is that each node is a conflux of streams
of different chemical compositions. With the depth of description accorded by the
integral mass and momentum conservation alone it is consistent to model the mixing
as instantaneous. The prevalent components are CH4 , O2 , CO2 , N2 . The molar
fractions
of the mixture
o input by each well are measured and recorded as a quadruple
n
xCH4 , xO2 , xCO2 , xN2 . By definition
xCH4 + xO2 + xCO2 + xN2 = 1.

(25)

12

Flow solution in landfill networks

Upon completion of mixing at the ith node the molar fraction xY of component Y of
the mixture entering the set of outgoing pipes Jiout is given by
!,
!
X
X
(i)
xY =
ṁj xY j + ṁ(i) xY
ṁj + ṁ(i) ,
(26)
j∈Jiin

j∈Jiin

(i)

wherein xY j and xY are the fractions in j th pipe and ith well respectively. The
verification that the newly obtained fractions sum to unity is immediate. Naturally
the computation of gas composition throughout the network is causal and must proceed
from the remotest well towards the flare. Given that the solution to the non-linear
system (24) is perforce iterative, the mixing is ineluctably performed at every iteration.
The gas constant R is then computed for each pipe
,
X
xY MY , j = 1, . . . , e
(27)
Rj = Ro
Y

with Ro being the universal gas constant, Y running over all components of the mixture
and MY denoting the corresponding molar weights.
3.2. Pressure distribution
The head loss across the network must be similarly computed at each iteration of the
numerical solution of (24). Starting at the flare, where a desired (design phase) or actual
(operational landfill) sub-atmospheric pressure is prescribed, the pressure is computed
sequentially upstream across the network with (23). Upon convergence the solution is
examined for compliance with the safety criteria. If the landfill is still under design, the
analysis determines where in the system the margin to maximal permissible pressure is
too small. Pipe geometry is consequently to be adjusted (larger diameter or additional
edges) until admissible margin is obtained. If the landfill is operational, when the
pressure limit is approached at any point in the network, the operator will adjust the
flow to the relevant wellheads, or if the problem appears to be of a more global nature,
modify the suction strength at the flare.
3.3. Example network
System (24) was solved with the underlying geometry of the network in figure 3.
The example features all properties pertinent to a landfill, rendering its generality
sufficient: (a) substantial number of nodes n = 23, edges e = 28 and fundamental
cycles e − n + 1 = 6; (b) exitence of isolated wells, namely Jiin = ∅ for i = 15, 17; (c)
a node designated as flare, J1out = ∅; (d) existence of crucial edges (j = 8, 15), whose
contingent disconnexion entails flow reversal in several parts of the network and possible
violation of the safety constraint delimiting the pressure to a maximal value throughout.
The following nominal real network geometry and landfill parameters were used:
 e
 e
80m ⩽ ℓj j=1 ⩽ 450m, 0.15m ⩽ dj j=1 ⩽ 0.2m, individual well volumetric flow

13

Flow solution in landfill networks
(i)

(i)

up to 0.04m3 /s, and gas composition with 0.45 ⩽ xCH4 ⩽ 0.55, 0 ⩽ xO2 ⩽ 0.02,
(i)
(i)
0.3 ⩽ xCO2 ⩽ 0.4, the fraction xN2 in accord with (25), i = 1, . . . , n. The suction at the
flare was set to 6kPa < pflare < 25kPa below atmospheric pressure, a feasible range for
the landfill in question. One should bear in mind that the solution network is a model
of the real system, thus the dimensions do not directly conform to the field installation.
For instance, the pipe lengths are effective values adjusted to account for losses due
to various valves, fittings, condensate drain devices etc. Moreover, each node might
represent a cluster of wells. Being an unalienable part of landfill flow computation,
such equivalences, extraneous to the specific head loss equation employed, are based on
experience of the design and construction engineers. Hence the graph layout ought to
be regarded as a framework, where numerous real pipework detail has been absorbed
into effective dimensions. For this reason all values are given in precision of no more
than two significant digits.
With the foregoing base parameters an increase in the following quantities was
tested: (a) pipe length; (b) elevation difference between nodes (equivalent to increase in
incline angle |α|); (c) well input flow rate. It was found that formula (10) consistently
led to overestimation of pressure at the nodes. For clarity of exposition the effect is
illustrated for the higher values of extraction vacuum at the flare.
Figure 4 shows the variation of pressure within the system relatively to the flare
along a path from the flare at node 1, corresponding to the lowest pressure, to well
15, the remotest point and highest pressure in the network. Discernible changes of
slope occur at wells with significant mass input. The full potential of extraction is
achieved when the maximal pressure is only slightly below atmospheric, i.e. the fluid
at the remote nodes recovers most of the vacuum imposed at the flare. As exemplified
below, this is where the use of (23) instead of (10) can make a difference, mitigating the
conservative prediction by the latter approach.
The solution of (24) with the head loss correctly given by (23) was compared
with the case, where (23) was replaced by the discrete approximation (10), the
implementation thereof based on the fluid density at the exit of the pipe and kept
fixed (for the sake of comparison only, in congruence with (10), albeit erroneously in
the current application). Figure 5 shows the resulting overestimation of pressure by the
traditional approach (10). Sufficiently far from the flare this effect is nearly uniform.
Observe the non-linear dependence between extracted mass and said overestimation:
increasing the well input by a factor of 1.5 resulted in a fourfold prediction error,
whilst a doubled input incurred a tenfold error. With the maximal permissible pressure
commonly set slightly below atmospheric, e.g. 1250Pa, overestimation of up to 625Pa
is significant. Whilst design guidelines usually suggest node vacuum of about 5kPa as a
reasonable safety margin, in reality the operators strive to maintain the limit pressure
as close to atmospheric as possible, forbearing to exert unnecessary suction liable to
draw unsustainable amounts of landfill gas, beget air infiltration and decimation of
the microorganisms responsible for biodegradation. The limit pressure is not set in
stone and depends on the operator’s experience with the fluctuations observed in each

14

Flow solution in landfill networks

pj − pflare [kPa]

25

20

15

10

5

0

0

3000
ℓ [m]

Figure 4: Pressure relatively to the flare along the path formed by pipes 1 to 14 (figure 3) for a total collection of 0.2 kg/s
(dotted), 0.3 kg/s (dashed)nand 0.4 kg/s (solid). Parameters
used: flare node
o n
o i = 1, pflare = 25kPa below atmospheric,
T = 13o C, molar fractions

xCH4 , xO2 , xCO2 , xN2

=

0.5, 0.01, 0.45, 0.04 .

locality. The extraction vacuum at the flare is usually significantly higher, although
that too might vary widely during the lifetime of a landfill, when for instance the cover
thickness or permeability is changed to accommodate addition of waste, landscaping
etc. The ultimate goal in controlling these parameters is to retain a productive colony
of anaerobic microorganisms.
From the vantage point of fluid mechanics the aforementioned findings are in
congruence with the essence of error introduced due to misapplication of (10), since
the computation is performed upstream from the flare (consult figure 2). Whilst the
result accords well with intuition, it falls far from triteness, as owing to the presence of
cycles in the network not every edge has the head change computed up the stream. In
fact the flow direction along each edge was determined automatically every iteration of
the solution. In a different application, where the reference pressure point is situated
upstream, or where several reference points exist, a result of a different nature might
ensue.
From the aspect of landfill design and operation the usage of (10) will result in more
conservative predictions than necessary. The main complexity of the numerical solution
of (24) and networks of this ilk stems from the non-linearity of pressure continuity
equations (24b) as well as the relative number of fundamental cycles and edges in the

15

Flow solution in landfill networks

∆pj [Pa]

625

500

375

250

125

0

0

3000
ℓ [m]

Figure 5: Variation between the classic and proposed methods: overestimation of pressure by the classic DarcyWeisbach model (equation (10), constant density) versus the continuous model (equation (23), correctly varying density)
along the path formed by pipes 1 to 14 (figure 3) for a total collection 0.2 kg/s (dotted), 0.3 kg/s (dashed) and
0.4
i = 1, pflare = 25kPa below atmospheric, T = 13o C, molar fractions
n kg/s (solid). Parameters
o
n used: flare node o
xCH4 , xO2 , xCO2 , xN2 = 0.5, 0.01, 0.45, 0.04 . For reference, maximal permissible pressure is commonly set slightly
below atmospheric, e.g. 1250Pa.

graph. In practice the replacement of (10) with the more complicated at first glance (23)
encumbers the solution but little (observe that a numerical solution of (2), for instance,
and (23) demands the selfsame computational tools), however brings forth the advantage
of a more accurate estimate. In particular, in order to direct the collected gas to
beneficial uses the landfill operator is under an obligation to supply designated amounts
of gas. As the maintained flow regime must comply with the maximal permissible
pressure limit, overestimated pressure values will result in superfluous restrictions on
the amount of gas extracted, undermining optimal collection and efficiency of the entire
facility.
Figure 6 depicts the variation of Mach number along the same path as in figures 4
and 5. In this example, as in the vast majority of cases, the flare node would give the
upper bound for the entire network. As expected, the Mach number is very low and
nearly uniform in the part of the network containing interconnected cycles, climbing to
the upper bound only in the acyclic part leading to the flare.

16

Flow solution in landfill networks

Mach number

0.08

0.06

0.04

0.02

0

0

3000
ℓ [m]

Figure 6: Mach number along the path formed by pipes 1 to 14 (figure 3) for a total collection of 0.2 kg/s (dotted),
o
0.3 kg/s (dashed)
used: flare node
n and 0.4 kg/s (solid). Parameters
o n
o i = 1, pflare = 25kPa below atmospheric, T = 13 C,

molar fractions

xCH4 , xO2 , xCO2 , xN2

=

0.5, 0.01, 0.45, 0.04 .

3.4. Kinked pipes
In landfill construction the occurrence of a pipe being kinked in the course of imbedment
within the substratum covering the refuse cavity is a contingency to contend with.
Mathematically a kinked pipe translates into elimination of the corresponding graph
edge. Depending on the location, this modification might precipitate momentous
changes for the entire network, possibly begetting flow reversal and flow rates manyfold
larger than initially expected. One must bear in mind that when heavy construction
machinery riding over the substratum damages a pipe laid within, the likelihood of
repairs being effectuated is minuscule, rendering the resultant arrangement a permanent
working regime. A gratuitously restrictive prediction of pressure distribution in the
reconfigured system due to misapplication of (10) instead of (23) will necessitate an
operator induced flow reduction at the wells in adherence to the limit of maximal
permissible pressure in the system. Even small changes in the total flow rate annually
amount to significant quantities of gas that might make a difference between the
feasibility of a profitable alternative fuel or methane content enhancement facility, and
the environmentally detrimental destruction by flaring with no useful work.

Flow solution in landfill networks

17

4. Discussion
A problem of flow in landfill gas collection network was considered. The fluid is
an ideal gas and the regime isothermal, thereby the flow is weakly compressible.
The traditional Darcy-Weisbach formula (1) developed for incompressible flows, whilst
providing a close estimate in a certain part of the parameter space, is conceptually
erroneous, the use thereof in conjunction with the state equation of ideal gas (3) setting
forth an inherent incongruence. On the other hand, related formulae in wide use in
the gas pipeline industry take into account compressibility effects insofar as to have
recourse to multifarious empiric loss coefficients. The suggested method generalises the
classic Darcy-Weisbach formula to include a continuous variation of density, yielding an
ordinary differential equation. The approach is based on integral mass and momentum
conservation equations, thence the result is obtained without undue premise on the
constancy of the Darcy friction coefficient.
In the case of ideal gas the differential equation is soluble in closed form comprising
non-linear algebraic equations (23), separate for the instances of horizontal and nonhorizontal flow. The equation involves the traditional Darcy friction factor. Its use was
justified by establishing an equivalence between a generic friction factor in the integral
momentum conservation equation and empirically sound as well as easily accessible in
practice Darcy friction factor. Any numerical algorithm suitable for solution of the
Colebrook equation (2) to obtain f is equally capable of solving (23).
The network was characterised by an incidence matrix, permitting to implement
the mass conservation equations at the nodes and pressure continuity equations in
matrix form. The feasibility of (23) supplanting (10) was confirmed by solution of a
flow in network geometry from a real landfill in the process of redesign and expansion.
It was seen that little additional computational complexity arose due to the change
of the head loss model, whereby it follows that the misapplication of the classic
Darcy-Weisbach formula (10) for isothermal, compressible flow cannot be justified by
the simplicity thereof. The solution stemming from the continuous approach herein
properly accounts for compressibility, whose effect is non-negligible over significant pipe
lengths. Furthermore, beyond the conceptually erroneous complexion of (10) in the
context of isothermal gas flow, twofold pragmatic implications ensue. First, from the
aspect of landfill operation the error in the estimation of pressure distribution in the
network entails consistently more conservative working regime with reduced flow rates.
Application of the conceptually correct (23) permits utilising the safety limits to their
full potential. Second, extraction of all available gas and accurate estimate of the mass
thereof allows for planning of environmentally more beneficial alternatives than flaring.
In summary, for a weakly compressible isothermal flow the Darcy-Weisbach formula
is inapplicable. With the correct, although more involved, non-linear equation, within
the context of pipe network analysis the additional computational complexity is minor.
The suggested approach enables more accurate engineering decision making as regards
efficient landfill operation and optimisation of gas collection.

18

Flow solution in landfill networks

pipe diameter
dj [m]

length
ℓj [m]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
16
16
17
18
19
20
21
22
23
24
25
26
27
28

200
200
300
300
175
150
80
350
120
330
200
180
450
300
220
100
200
150
150
150
210
160
200
185
185
80
185
180

0.15
0.15
0.15
0.2
0.2
0.2
0.2
0.2
0.2
0.2
0.2
0.2
0.2
0.2
0.2
0.2
0.15
0.15
0.15
0.15
0.15
0.15
0.15
0.15
0.15
0.15
0.15
0.15

node elevation
hi [m]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

10
10
10
10
10
10
10
10
15
15
15
12
15
8
15
12
12
15
15
15
15
15
15

flow rate
ṁ(i) [kg/s]
0
0.005
0.007
0.012
0.002
0.007
0.02
0.02
0.016
0.005
0.005
0.005
0.005
0.01
0.007
0.01
0.01
0.01
0.01
0.01
0.005
0.005
0.01

Table A1: Network details for the example in figure 3. Pipe data (left): internal pipe diameter dj and effective pipe
length ℓj , j = {1, . . . , 28}. Node data (right): elevation hi and input mass flow ṁ(i) , i = {1, . . . , 23}.

Appendix A. Example network details
Table A1 lists the pipe and node data used for the example network given in figure 3.

Flow solution in landfill networks

19

Acknowledgement
The authors are grateful to Richard Brewster (Thompson Rivers University) for
insightful discussions.
References
Batchelor G 1990 Introduction to fluid dynamics Cambridge.
Benedict R 1980 Fundamentals of pipe flow John Whiley & sons.
Benedict R 1983 Fundamentals of gas dynamics John Whiley & sons.
Coelho P and Pinho C 2007 J. Braz. Soc. Mech. Sci. Eng. 29, 262–273.
Colebrook C 1939 Journal of the Institution of Civil Engineers 11, 133–156.
Conestoga-Rovers A 2010 Landfill gas management facilities design guidelines Technical
report British Columbia Ministry of Environment.
Corless R, Gonnet G, Hare D, Jeffrey D and Knuth D 1996 Adv. Comp. Math. 5, 329–359.
Crane C 1982 Flow of fluids through valves, fittings and pipe. Technical paper No. 410M
300 Park Avenue, New York, NY 10022, appendix A-5.
Godsil C and Royle G 2001 Algebraic graph theory Springer, New York.
Menon E S 2005 Gas pipeline hydraulics CRC Press.
Moody L 1944 Transactions of ASME 66, 671–684.
Nec Y and Huculak G 2015 Can.J.Civ.Eng. 42, 563–569.
Ward-Smith A 1980 Internal fluid flow Clarendon Press.
White F 1999 Fluid Mechanics McGraw Hill 4th edition.