Convective Acceleration in Porous Media

Convective acceleration occurs in porous media flows due to the spatial variations of the nonuniform flow channel geometry of natural pores. This article demonstrates that the influence of convective acceleration in a nonuniform a pore channel is analogous to that of a constricting pipe channel. Their fluid mechanical behaviour can be comparable, provided that their geometrical characteristics are described precisely in the same manner, and from the same point of reference with regards to the fluid velocity in the flow channels. The analogy of the dissipation mechanisms in nonlinear porous media flow to the "minor loss" approach in fluid mechanics of pipes is therefore appropriate. Conventional nonuniform pipe channel geometries obtain dissipation coefficients within the range 0 < CKL < 0.2. These pipe geometries are relevant reference points for natural porous media, and it is thus expected that most natural pore geometries will obtain values within this range. This assumption holds true for the nine different 3D porous media samples presented here. However, the results show that the rate of change in the pore geometry, and consequently the magnitude of induced convective acceleration, depends on: the area ratio a of the pore channel, the angle of approach θ and the rounding of the pore channel geometry. The rounding of the pore channel reduces the dissipation coefficient, as the rate of change becomes smoother along the channel length. The results also indicate that the pore tortuosity increase the magnitude of nonlinear dissipation, in good agreement with pipe flow behaviour. This knowledge can help improve our interpretation of experimental data and enhance the predictability of porous media equations that incorporate the appropriate dissipation coefficients CKL as a variable. The analogy of porous media flow to the "minor loss" approach in fluid mechanics of pipes is appropriate, and the angle of approach θ and the area ratio a of the pore channel govern the magnitude of induced convective acceleration in porous media The rounding of the pore channel geometry reduces the magnitude of induced convective acceleration The tortuosity of a pore influences the dissipation coefficient CKL and increase the magnitude of induced convective acceleration The analogy of porous media flow to the "minor loss" approach in fluid mechanics of pipes is appropriate, and the angle of approach θ and the area ratio a of the pore channel govern the magnitude of induced convective acceleration in porous media The rounding of the pore channel geometry reduces the magnitude of induced convective acceleration The tortuosity of a pore influences the dissipation coefficient CKL and increase the magnitude of induced convective acceleration


Introduction
Inertial forces and convective acceleration have a significant influence on the dissipation of mechanical energy in fluid flow regimes characterized by relatively high flow velocities. These forces induce nonlinear pressure responses that are not accounted for by the linear Darcy-equation (Darcy 1856). The second-order Forchheimer-equation is needed to explain these effects on the pressure gradient (Forchheimer 1901). Such flow conditions in porous media are relevant in many practical engineering problems. Particularly for the determination of the performance and the capacity of production wells and injection wells, i.e. in the petroleum industry (Huang and Ayoub 2008;Barree and Conway 2004), in groundwater wells utilized for domestic and irrigational purposes (Houben et al. 2018;Wen et al. 2011;van Lopik et al. 2019), and for groundwater heat pump applications (Gjengedal et al. 2021).
The inertial resistance factor β of the Forchheimer equation is a geometrical property that describes how the pore channel shape influence the nonlinear flow characteristics. It is challenging to determine this factor without performing field and laboratory tests. In-situ or experimental tests are therefore often routinely performed to determine well performances in the petroleum industry (Huang and Ayoub 2008). However, in many civil engineering projects there is limited time and funding available for detailed in-situ testing and laboratory tests, and many engineering projects therefore depend on empirical equations for prediction purposes. Most empirical equations used in soil engineering rely on soil characterization rather than on pore space characterization (van Lopik et al. 2019(van Lopik et al. , 2017Chapuis 2012). The Forchheimer equation is also frequently utilized in numerical modelling and computational fluid mechanics (Guo et al. 2019). Continuous use of these equations argues that further study on this topic is of scientific interest and is useful for many industries.
Convective acceleration occurs due to spatial variations of the flow channel geometry. This paper explores how geometrical properties of a pore channel induced convective acceleration and how this affects the inertial losses. Recent laboratory studies on 3D fabricated particle beds with predefined pore geometries enables more detailed and controlled analysis of the geometrical properties that are important for porous media flows (Huang, et al. 2013;Gjengedal et al. 2020). Gjengedal et al. (2020) has hypothesized that porous media flows exhibit similar behavioural traits to that of pipe flows in situations where the pipe channels are nonuniform, the so called minor losses approach in fluid mechanics. A method of comparison with the pipe flow equations is here presented, and the similarities between the porous media flows and pipe flows are further assessed.

Theory
The general Forchheimer equation corresponds to Eq. (1) for horizontally directed, one-phase, incompressible fluid flow in pore channels. Here the pressure dissipation ΔP (Pa) across the superficial length of the porous medium L (m) is proportional to a first-order viscous term and a second-order inertial term. The viscous resistance is proportional to the superficial velocity of the flow u s (m/s), the viscosity of the fluid μ (Pa·s) and is inversely proportional to the permeability of the porous medium k (m 2 ). The inertial resistance is proportional the inertial resistance factor β (m −1 ), the density ρ (kg/m 3 ) of the fluid, and the superficial velocity squared. The superficial fluid velocity is here defined as the bulk fluid flow rate Q (m 3 /s) divided by the bulk cross-sectional area A (m 2 ) of the porous media oriented perpendicular to the superficial flow direction.
The permeability k and the inertial resistance factor β represent intrinsic geometrical properties that describe the size and shape of the porous media pore channels (Ruth and Ma 1992;Yazdchi and Luding 2012;Whitaker 1996;Gjengedal et al. 2020). These geometrical properties are typically represented by a diameter of the solids that surround the pores (Chukwudozie and Tyagi 2013;van Lopik et al. 2019). However, the permeability and the inertial resistance factor should be identified with respect to the shape and size of the pores within the pore matrix. A new semi-analytical approach is presented in detail by Gjengedal et al. (2020) and provides Eq. (2). This equation describes the fluid flow through a single pore in relation to the average flow velocity in the pore body region. The permeability and the inertial resistance factor are here represented by Eqs. (3) and (4), respectively.
The geometry of the pore is described by the porosity n (m 3 /m 3 ), the specific surface of the pore S (m 2 /m 3 ), the tortuosity factor τ (-), and the pore body shape factor k 0 (-). The dissipation coefficient C d is not yet clearly defined in the porous media literature. It is suggested that the C d depends on the constriction ratio and the streamlining shape of the pore channel, like the "minor loss" coefficient in the literature of fluid mechanics of pipes, Fig. 1 (Gjengedal et al. 2020).
The concept of minor losses in pipes is an experimentally developed approach for determining the hydraulic loss contribution of various pipe components that have an obstructive effect on the flow, such as pipe bends, diffusers and flow controlling valves. These head loss components are termed "minor" because, in typical piping systems with long sections of pipes, they occur relatively infrequent in the pipeline and these losses are consequently minor compared to the head losses induced by the long straight pipe section itself (consequently the "major" losses). However, in piping systems with numerous turns, bends and valves over short distances, the minor losses can be far greater than the major losses. The corresponding equation for pipes is that of Eq. (5) (Çengel and Cimbala 2010). (1) (2) The total losses of hydraulic head h L (m) through a pipe segment of length L (m), and diameter D (m) is a sum of two terms that are proportional to the average velocity squared V s (m/s) and inversely proportional to the acceleration due to gravity g (m/s 2 ). The liner term, represented by the Darcy-Weisbach friction factor component f, corresponds to the viscous friction in straight and uniform pipe channels. The second-order term, represented by the loss coefficient K L , express the minor losses of each individual obstructive component in the piping layout. The sum of these coefficients account for the total inertial losses that are induced by convective acceleration through the entire pipeline. An obstruction in the flow path of the fluid might e.g. be the contraction and expansion of the pipe channel geometry seen in Fig. 1.
Pipe manufacturers readily provide experimental data for the K L of various pipe components. The K L of a single pipe component is calculated from experimental data via Eq. (6) (Crane Company, 2009).
The loss coefficient is treated as a product of an inertial resistance factor f T times the ratio of the length of the pipe segment L and the diameter D of the chosen representative pipe size within this segment. The smaller pipe size is typically selected as the reference size and the average flow velocity V s in Eq. (5) is based on the flow velocity in the smaller pipe segment. The letter "T" infers that the inertial resistance factor f T represents turbulent flow conditions. Including Eq. (6) into Eq. (5) produces Eq. (7) with the addition of the Darcy-Weisbach friction factor for circular pipes (f = 64/Re). The hydraulic head is now represented in terms of the corresponding pressure differential ∆P (Pa), instead of heads.
The Reynolds number (Re) for pipe flow in Eq. (8) relate to the pipe diameter D (m), which is the scientific convention in pipe flow engineering problems. However, if a pipe has a nonuniform and non-circular geometry the characteristic shape of the pipe is not related to the diameter of the pipe alone. Therefore, to evaluate nonuniform pipe channels, the fundamental characteristic shape definition must be applied instead (Schiller 1923), corresponding to Eq. (10).
This necessitates a correction to the conventional loss coefficient K L of Eq. (6). It is evident from Eq. (10) that to utilize the pipe diameter D as the characteristic length unit must produce a loss coefficient four times greater than would otherwise occur should the m pipe unit be applied instead. This numerical difference is not obvious because it is imbedded in the experimentally derived K L value. The corrected pipe flow equation therefore becomes Eq. (11), where the characteristic length unit m pipe describes the shape of the pipe channel along the pipe length L. Equation (2) for porous media and Eq. (11) for pipes display apparent similarities now that they are expressed in terms of a single pore channel and pipe channel. However, a direct comparison of Eqs. (2) and (11) necessitates that these two equations are organized in the same manner, where the fluid flow is evaluated from the same point of reference within the flow field in the pore channel and pipe channel.
Equation (2) must also be expressed by the characteristic length unit m in the most fundamental form through Eq. (12), and the average flow velocity must be expressing by Eq. (13). The letter "l" specifies that the average velocity V l (m/s) is the velocity of the flow with respect to the "larger" channel opening within the pore, corresponding to the pore body region of the pore. This provides Eq. (14), which is more fitting for the problem at hand.
(12) m = porosity of the pore unit volume spesific surface of the pore unit volume = n S .
The resemblance of the second-order terms of Eqs. (11) and (14) is now apparent. The comparison is particularly fitting in the special case with the tortuosity of the pore equal to unity (τ = 1.0). However, a direct comparison requires that a correction must be performed, to either one of Eq. (11) or Eq. (14), since the reference velocity for both equations is not yet the same. The reference velocity (V s ) in piping system is based on the velocity of the smaller pipe segment, while the reference velocity (V l ) in the pore channel is based on the fluid velocity in the pore body region of the pore. Gjengedal et al. (2020) therefore suggested that a correction should be performed on Eq. (14) for the inertial term according to Eq. (15), where the cross-sectional area ratio a (-) is applied in compliance with the conventions in fluid mechanics of pipes (Crane Company 2009;Idelchik 1994). The comparable porous media equation is therefore Eq. (16).
The dissipation coefficient C KL of porous media is deemed similar to the inertial resistance factor f T for pipes and should therefore display similar behavioural traits with respect to the constriction ratio and the streamlining shape of the flow channel. Both factors should also display similar traits with respect to the flow characteristics that dominate in laminar and turbulent flow conditions. It is known that the minor loss coefficients of pipes display an increased numerical value in laminar f L flow conditions versus turbulent f T flow conditions, f L > f T (Idelchik 1994;Çengel and Cimbala 2010). The flow characteristics in Forchheimer flow is laminar, and it is therefore the laminar inertial resistance factor (f L ) that should be similar to the dissipation coefficient C KL .
The suggested analogy for convective acceleration in porous media flows and laminar flow in nonuniform pipes can be verified through evaluating the legitimacy of Eq. (17) with published experimental data.

Method
A comparison is performed on known C KL and f L coefficients derived from experimental data found in the literature. Two independent experimental works on different 3D fabricated porous media provide the C KL data (Huang et al. 2013;Gjengedal et al. 2020). These porous media were constructed by stacking spherical and octahedra particles in the predefined packing configurations shown in Fig. 2. These packing arrangements produce a pore in between the solids that correspond to the shape of the pores shown in Fig. 3. In the experiments of Huang et al. (2013) the spherical balls in the Samples T1-T4 were fabricated by gluing acrylic ball together in cubic containers. These balls were smooth, and the corresponding pore geometries in their experiments have a similar shape as shown in Fig. 3.
Unlike the smooth acryl balls T1-T4, the Sample A-D particles of Gjengedal et al. (2020) are very rough due to the 3D printing technology employed in the fabrication. This caused the actual shape of the pores to be slightly different than the ideal and smooth configurations seen in Fig. 3. The rough samples have lower porosities, larger degree of contraction, and rougher surfaces than shown in the ideal sketches in Fig. 2. The estimated geometrical properties and the measured geometrical properties are reproduced in Table 1 for comparison, along with the calculated dissipation coefficient.
The spherical particle configuration has the most severe constriction ratio of all the packing configurations. It is assumed that the smooth acrylic balls T1-T4 have the same area ratio as the ideal, a = 0.215, while the A3 sample has a much more severe constriction ratio of a = 0.128, due to the rough particle surface. Still, the A3 pore channel has a similar curved pore shape to that of T1-T4. Curved and rounded channels tend to reduce the "minor losses" observed in pipe flow experiments (Idelchik 1994;Crane Company 2009;Çengel and Cimbala 2010). This geometrical trait is reflected in the dissipation coefficient, which is lowest for the spherical samples A3 C KL ≈ 0.03 and T1-T4 C KL ≈ 0.033-0,048. Note that the T1-T4 and A3 samples have similar numerical C KL values, regardless of the particle roughness.

Table 1
The measured geometrical properties of the 3D porous media A3, B3, B6, C1, and D3 of Gjengedal et al. (2020)   It is conventional in pipe flow experiments it to assign an angle of approach θ that describes the rate with which the contracting geometry changes towards the narrowest section of the flow channel (Crane Company 2009;Gibson 1910;Gibson 1911Gibson -1913. For these 3D geometries the narrowing does not occur evenly along the pore channels, and there will not be a single angle of approach that uniquely describes the rate of contraction. However, an apparent angle of approach θ* can be assigned based on the uniform angle that would result from an inscribed circle within the pore throat cross section and the pore body cross section. This angle allows for an intuitive comparison between the 3D sample types. The apparent angle of approach for A3 is approximately θ* ≈ 82° and for the T1-T4 samples the θ* ≈ 75°.
The octahedron Samples B, C and D are also fabricated with a rough surface and the actual shape of the pores are slightly different than the ideal and smooth configurations in Fig. 3. The ideal configuration of the B and C samples should produce the same area ratio of a = 0.500 regardless of the particle size. However, the C1 sample is relatively rougher than the B and D samples, which produce a more severe constriction ratio a = 0.290 than the B3 and B6 samples a ≈ 0.42. Consequently, the apparent angle of approach is more abrupt for the C1 θ* ≈ 71° than for B3 and B6 θ* ≈ 61°. A slightly higher dissipation coefficient of C1 C KL = 0.11 than the B3 and B6 samples C KL = 0.09 also arise.
The D3 sample has the least constricting area ratio a = 0.640 and the least abrupt apparent angle of approach θ* ≈ 38°. Yet, the channel pathway of sample D3 is offset by an interfixed octahedron particle and the pores have a tortuosity of τ = 1.2247. An alteration of the predominant flow direction tends to increase the "minor losses" observed in pipe flow experiments (Idelchik 1994;Crane Company 2009;Çengel and Cimbala 2010). Gjengedal et al. (2020) suggest that this trait is reflected in the D3 sample having the largest C KL = 0.13.
For comparison with pipe flow equations it is necessary to evaluate a pipe geometry that has similar geometric qualities to these pores. This is a pipe with a constricting segment along the flow channel and must consist of a contracting part and an expanding part, like the symmetrical half-cone pipe shown in Fig. 1A. Crane Company (2009) and Rennels and Hudson (2012) have summarized experimental data for the loss coefficients of such halfconed pipe segments and provide semi-empirical equations. The area ratio a and the angle of approach θ determine the shape and length of the pipe channel. The minor loss coefficients of pipes can be determined via Eq. (18-22). The minor losses depend on the flow direction through the half-cone, where Eq. (18-20) apply for the expanding cone segment (Rennels and Hudson 2012;Crane Company 2009).
Note that the equations are rendered in terms of the standardized area ratio a, while most pipe literature present these equations in terms of pipe diameter ratio, corresponding to the root of the standardized area ratio. Equation (21-22) apply for the contracting cone segment (Rennels and Hudson 2012).
where The minor loss coefficients of Eqs. (18-22) are typically determined in fully developed turbulent flow conditions, where the inertial friction factor (f T ) is independent of the Reynolds number and where the viscous contribution to the dissipation is negligible. However, in laminar flow conditions the minor loss coefficients show an increased numerical value (Idelchik 1994). A correction of the minor loss coefficient is typically applied by introducing the coefficient α L,T (-), which is close to unity for fully turbulent flow (α T ≈ 1.05) but is two for fully developed laminar flows (α L ≈ 2.0) (Çengel and Cimbala 2010). This phenomenon also occurs in porous media when the onset of turbulence is observed in experimental data. The data typically display a reduction in β of the Forchheimer equation as the flow velocity increases beyond the Forchheimer flow regime (Barree and Conway 2004;Fand et al. 1987;Skjetne and Auriault 1999).
Equations (18)(19)(20) are plotted in Fig. 4 and Eq. (21) in Fig. 5 and show the trends of the minor loss coefficients K L with respect to the contraction ratio and the angle of approach. The contraction losses and the expansion losses increase when the θ increase and when the a becomes progressively more constricting. The expansion losses are typically larger than the contraction losses for the same geometrical shape but are of similar magnitude for a 0,5 > 0.60. This demonstrate that the direction of flow through a constricting pipe segment has an influence on the induced convective acceleration.
Equations (18-21) are used to calculate the corresponding dissipation coefficient C KL and the inertial resistance factor f L for these pipe geometries. The geometrical properties of a half-coned pipe segment: the pipe channel length L pipe , the pipe channel volume n pipe and the pipe channel internal surface area S pipe can be expressed by the diameter size of either one of the pipe entrances. The larger pipe size is here selected as the reference size (A large in Fig. 1) and calculated via Eq. (23, 24, and 25). The numerical value of the pipe size D large does not affect the outcome of the calculations, because the ratio of m/L in Eq. (17) render the outcome independent of the pipe size. The only geometrical variables that affect the outcome of Eq. (17) is the area ratio a and the angle of approach θ, which demonstrate that the inertial friction factors f L and f T are a function of the shape of the flow channel, but not the size of the flow channel.
If the contracting segment and the expanding segment of the pipe are not symmetrical the length units, L and m, must be calculated separately for both the contraction and the expansion segment of the channel. However, for a symmetrical half-cone pipe the angle of approach θ is equal in both segments. The comparable pipe flow expression then becomes the sum of the contraction loss and the expansion loss of the pipe, shown in Eq. (26).
This equation is utilized to determine the pipe dissipation coefficients (C KL-pipe ) for comparison with the porous media data.

Results
The C KL-exp due to expansion of a pipe is presented in Fig. 6 for the range of area ratios 0.01 < a < 0.98 and the range of angle of approaches 0° < θ < 170°. The figure demonstrates that a large angle of approach induces larger dissipation coefficients than smaller angles, in good agreement with the general behaviour of the minor loss coefficient (K L-exp ). However, unlike the K L-exp data that plot almost equally for all angles above 45°, the C KL-exp provide individual curves for each individual angle of approach. The C KL-exp does not depend on the area ratio a in the same manner to that of K L-exp in Fig. 4. The C KL-exp display a peak value within the 0.2 < a < 0.3 range (slightly dependent on the angle of approach). When the area ratio approaches zero the dissipation coefficients trends towards a unique value depending on the angle of approach. On the other hand, the C KL-exp drops sharply towards zero when the area ratio increases and approach unity. In the special case with a = 1.0 the C KL-exp reaches a critical point and becomes zero. This demonstrate that the C KL-exp will always affect the fluid flow through the channel except when a = 1.0 and there is no change in the flow channel geometry.
The C KL-con due to contraction of a pipe is presented in Fig. 7 for the range 0.01 < a < 0.98 and the range of angle of approaches 0° < θ < 170°. The figure demonstrates that larger angles of approach induce larger dissipation coefficients than smaller angles, which is in good agreement with the minor loss coefficient K L-con . However, opposite to the K L-con data that have largest values at small area ratios, the C KL-con display a reversed trend where the peak values occur for relatively small constricting area ratios (0.6 < a < 0.7).
The dissipation coefficients for the symmetrical half-cone pipe geometry are obtained by combining the data of Figs. 6 and 7. This is presented in Fig. 8, in addition to the dissipation coefficients of the T1, T2, T3, T4, A3, B3, B6, C1 and D3 samples. The combined effect of the C KL-exp data and the C KL-con data, produce a "bell-shaped" curve.  demonstrates that a contracting and expanding pipe segment can obtain dissipation coefficients as high as C KL ≈ 0.2. The porous media dissipation coefficients plot within the same region as those of the calculated pipe data.
However, the porous media data plot differently than what the apparent angle of approach (θ*) would suggest. The C KL of T1, T2, T3, T4 with θ* ≈ 75°and A3 sample with θ* ≈ 82° is here shown to plot equal to a comparable pipe segment with θ° = 30-40°, arguably much lower than what would be expected. The opposite is observed for the octahedron samples C1, B3, and B6 whom all show higher C KL values than for comparable angle of approaches in pipes θ° = 85-95°. The differences are particularly profound for the D3 sample with θ* ≈ 38°, where the C KL is shown to plot equal to a pipe segment with θ° = 125°. The larger C KL value of C1 than for B3 and B6 is consistent with C1 (θ* ≈ 71°) having a more abrupt constriction ratio than the other two (θ* ≈ 61°). The same phenomenon occurs for the pipe data where larger θ° consistently produce larger C KL values for the same area ratio.
The data in Figs. 6, 7 and 8 are calculated with the smaller channel velocity V s as the reference velocity. This point of reference visualizes how the angle of approach affect the C KL , more so than the area ratio. The data of Fig. 8 demonstrate this particularly well for very small angles of approach, e.g. θ° = 10°, which shows that the dissipation coefficient will be very small, e.g. C KL ≤ 0.001, regardless of the area ratio a.
The larger channel velocity V l provides another point of reference that, in contrast to the V s , visualize more clearly how the area ratio affect the magnitude of dissipation. The corresponding dissipation coefficient C d calculated from Eqs. (17) and (26) for the symmetrically half-cone pipe segments is presented in Fig. 9 along with the corresponding porous media data.
It is evident that the C d values are much larger than the corresponding C KL values and Fig. 9 demonstrates more clearly how the constricting area ratio a influences the flow characteristics. For instance, even for very small angles of approach, e.g. 10°, Fig. 9 shows that the C d will have a significant influence on the convective acceleration term if the area ratio is small enough, e.g. when a < 0.1. This is not as clear in the C KL data in Fig. 8 where the corresponding C KL values are very small. This effect is reflected in the porous media data, where the A3 sample now has the largest C d and the D3 sample has the lowest C d of all the samples, in good agreement with the data response of the experimental flow tests of Gjengedal et al. (2020) and Huang et al. (2013).

Discussion
The presented work compares pipe flow equations and porous media flow equations. It is demonstrated that the influence of convective acceleration in pipe channels and pore channels are comparable if their geometrical aspects are described precisely in the same manner, and from the same point of reference with regards to the fluid velocity in the channels. A contracting and expanding pipe segments obtain dissipation coefficients that plot within the same region as dissipation coefficients obtained from similar contracting and expanding porous media geometries. This is logical because both approaches deal with the same geometrical aspects, namely the expansion and contraction of a flow channel.
The comparison of the pipe flow data and the porous media data reveals some important geometrical characteristics and their influence on convective acceleration. From Fig. 8 it is shown that the dissipation coefficient C KL is only zero in the special case of a uniform flow channel geometry, when a = 1.0. This infers that the convective acceleration term can be ignored when the pore or pipe geometry is equal to a straight and uniform flow channel but must be accounted for in all other geometrical cases. For all other geometrical cases it can be expected that many relevant pore channel geometries in natural porous media will 1 3 obtain dissipation coefficients lower than C KL < 0.2, as is observed for these pipe geometries in Fig. 8.
The rate of change in the geometry along the channel length governs the magnitude of convective acceleration that develops in the flow. The rate of change is here expressed as a function of the angle of approach θ and the area ratio a of the flow channel. The pipe data and the porous media data in Figs. 8 and 9 show dissipation coefficients C KL and C d that have similar numerical values and trends with respect to the angle of approach θ and the area ratio a of the pipe flow channels. An abrupt constriction in pipes, with a large angle of approach θ induce more convective acceleration than does a gradual constriction where the angle of approach θ is small. The same trend occurs in the porous media data and can be observed by comparing the C KL data of C1 θ* ≈ 71° versus B3 and B6 θ* ≈ 61°. This is reasonable because the rate of change in the geometry will then be greater over shorter distances.
The T1-T4 and A3 samples achieve lower C KL values than the apparent angle of approach θ* would suggest when compared to the pipe C KL data. The opposite is seen in the C1, B3 and B6 samples with relatively higher C KL value than the pipe data. This indicate that the rounding of the channel geometry play an important role in reducing the numerical values of the C KL coefficient. These observations are in good agreement with the observed trends in other pipe flow experiments (Idelchik 1994;Crane Company 2009;Çengel and Cimbala 2010;Rennels and Hudson 2012). This implies that an abrupt and sharp change in the flow cross-sectional area induce more convective acceleration, than does a rounded and gradual change in the flow channel. This is a logical response because the rate of change in the geometry is greater over shorter distances for abrupt and sharp corners.
From Fig. 8 it can be expected that the most relevant straight, nonuniform flow channel geometries in natural porous media will obtain dissipation coefficients within the range 0 < C KL < 0.2. However, the influence of pore tortuosity is not directly accounted for in this study because the presented pipe data apply to straight pipes only. It is recognized from pipe flow experiments that an alteration of the predominant flow direction tends to increase the "minor losses" in pipe bends, even without contraction or expansion of the pipe channel (Crane Company 2009;Çengel and Cimbala 2010;Idelchik 1994). By comparing the D3 data with the pipe data in Fig. 8 it is evident that the C KL value of D3 is higher than the C KL of straight pipes with a similar angle of approach. The D3 angle θ* ≈ 38° is here shown to plot equal to a pipe segment with θ = 125°. This suggests that the induced convective acceleration in D3 originate from other sources than the contraction and expansion of the flow channel alone, and one such source might be the influence of the tortuosity. Further studies should evaluate this phenomenon by including a comparison with pipe bend data so that the influence of the tortuosity on the C KL can be accounted for.
An interesting characteristic can be observed in Eq. 26 with regards to the role of symmetry. If the pore channel is symmetrical the sum of the dissipation coefficients from the expanding and contracting segments will be independent of the direction of flow. However, if the pore channel is asymmetrical on either side of the obstruction, the sum of the dissipation coefficients will depend on the direction of flow because the sum of the two terms might be directionally dependant. This is shown in Eq. 26 where, for an asymmetrical pore channel, the length units L and m must be calculated separately for both the contraction and the expansion segment of the channel, which in turn will have different values of K L-exp and K L-con depending on the direction of flow. The resulting over-all effect on the dissipation coefficient C KL and C d will thus be larger in one flow direction than in the other flow direction through the pore.
Naturally occurring pores are rarely symmetrical, e.g. since the pore geometry in soils is determined by the packing of irregular particle grains. The role of symmetry therefore has practical relevance in many engineering problems that deal with non-Darcy flow conductions, particularly with regards to well performance testing of pumping wells and injection wells. One relevant example is the Aquifer Thermal Energy Storage (ATES) system where a groundwater well can be used for both groundwater production and groundwater injection (Abuasbeh and Acuna 2018). If an ATES well is operated under non-Darcy flow conditions, it might have a different hydraulic capacity depending on the direction of flow, e.g. into or out from the well, and it is thus not necessarily reasonable to assign the same well capacity in pumping mode as in injection mode. Similarly, conducting a pumping test in a well is not necessarily the appropriate testing method for evaluating its injection capacity. This might play a role e.g. when old oil wells, originally used for oil production, will be used as CO 2 injection wells in the future.
It is a scientific challenge to determine the shape of the pore channel in naturally occurring porous media. This requires extensive laboratory work, especially to evaluate the angle of approach θ* and the area ratio a of the flow channel, which arguably makes it challenging to utilize the C KL coefficient for prediction purposes or for comparison with other experimental data where these geometrical traits are unidentified. For practical applications it is probably more convenient to evaluate the combined influence of θ and a from the C d data of Fig. 9 rather than the C KL of Fig. 8. Still, the presented results might serve as useful reference points for numerical modelling, especially when more experimental data becomes available for other pore geometries in the future. Consulting Eq. 2 to evaluate the geometry of the pore space in e.g. numerical models might therefore be useful from an interpretational point of view.

Conclusion
This work demonstrates that the influence of convective acceleration in nonuniform pore channels is analogues to convective acceleration in nonuniform pipe channels. The two phenomena are comparable if their geometrical characteristics are described precisely in the same manner, and from the same point of reference with regards to the fluid velocity in the channels. Studies on pipe geometries are consequently relevant for comparison with porous media studies. Data interpretation can then be improved as more experimental data and knowledge regarding the influence of the channel geometry become available from future studies and this will help our data interpretation and improve the predictability of porous media equations.
The presented results suggest that most nonuniform straight pipe flow channel geometries obtain dissipation coefficients within 0 < C KL < 0.2. These geometries are relevant for natural porous media and it is expected that most pore geometries will obtain values within this range. However, the dissipation coefficient is a function of the rate of change of the geometry along the channel length. The results show that the rate of change depends on: • The angle of approach θ and the area ratio a of the pore channel. The C d data of Fig. 9 and the C KL of Fig. 8 demonstrate how these geometrical traits influence the magnitude of induced convective acceleration • The rounding of the pore channel geometry reduces the magnitude of induced convective acceleration • The tortuosity of a pore influences the C KL and increases the magnitude of induced convective acceleration Further work should investigate the influence of other pore geometries and, especially the influence of the pore tortuosity. A comparison with pipe bend data might reveal the influence of the pore tortuosity on the C KL and C d .

Appendix A
The experiments of Huang et al. (2013) were conducted on acrylic balls of uniform size and shape bonded together by chloroform to the produce a porous media of cubic packing of smooth spheres. Four different diameters d of acrylic balls were chosen for these tests; the T1 = 3.0 mm, T2 = 5.0 mm, T3 = 8.0 mm, and T4 = 10.0 mm. This data is here evaluated via the methods suggested by Gjengedal et al. (2020) to produce the dissipation coefficient C KL for each sample.
For the evaluation of their data, it would have been preferable if the porosity n, specific surface area S and area ratio a were measured directly, e.g. with a CT scanner and a 3D visualization software. However, in their paper the porosity n of the samples is not measured but is assumed to be equal to the ideal cubic sphere packing porosity, with n = 0.476, regardless of the sphere diameter. The specific surface area S of the pores is calculated indirectly via Eq. (27). This geometric relation is correct for the T1-T4 samples if the spheres are smooth and if the fabrication has ensured no alteration of the sphere shape and the shape of the corresponding pore geometry.
In ideal conditions these assumptions are reasonable, but it is evident from the pictures of the spheres (Fig. 2 of their paper) that the melting of the acrylic balls have caused a slight merging at the contact points of the spheres. These geometric discrepancies can lead to an erroneous correlation with the experimental results. However, without supplementary measurements of the pore structures it is currently not easy to evaluate the data further.
The data of Huang et al (2013) are presented in Table A1 and are used to calculate the corresponding dissipation coefficient C KL for samples T1-T4. The fluid properties of water are found via reversed calculation of the data found in Table 1 of their paper and reveal that  (28) and (29), where the linear coefficient α and the polynomial coefficients β are adjusted to best fit the experimental data. Visual comparison with other empirical equations in Fig. 9 of their paper show that the linear coefficient should be α < 144 and the coefficients β should be close to 0.6 or slightly higher.
To calculate the dissipation coefficient C KL one must first calculate the pore shape factor k 0 from Eq. (30). It is here assumed that the Stokes constant (3π) is a good approximation for the spatial strain ratio within the fluid flow field of the pores, as suggested by Gjengedal et al. (2020). The tortuosity is unity for the cubic packing of spheres and does not affect the calculations (τ = 1.0).
If the pores within samples T1-T4 are equal in shape to the ideal cubic packing of spheres, they should ideally produce equal α-values and β-values regardless of sphere size. In Table A1 it is shown that the estimated linear coefficients α agrees well with the plotted figures and graphs of Huang et al. (2013) where the α < 144. However, the largest sphere samples T3 and T4 attain smaller coefficient α than does the T3, and especially the T4 sample. The slight variation in the coefficients α of the four samples suggest that the shape of the pores in the four samples are slightly different and are altered compared to the ideal case.
A qualitative evaluation of these results suggest that this is due to the fabrication method. The bonding with chloroform melts the acrylic balls together, altering the ball's shape slightly at the contact points with the neighbouring spheres. This melting causes a slight, but observable distortion of the pore geometry compared to the ideal cubic packing of loose individual spheres. The size of the merged section between the acrylic balls will obviously depend on the relative amount of chloroform used in the fabrication. In practise, it is likely it was utilized relatively more amount of chloroform per individual sphere in the fabrication of the small spheres (T1 and T2) than for larger spheres (T3 and T4), at least if the same equipment was used for the fabrication of all four sample sizes. The size of the melted section between the acrylic balls should thus depend on the ball diameter, and the bonded section would be relatively larger for the small sphere sizes. In Fig. 2 of Huang et al. (2013) the melted section of acryl can be observed in several pictures of one sample. It is not specified which sphere size these pictures portray, but it is likely the T3 = 8.0 mm or the T4 = 10.0 mm sample size. If the alteration at the contact points is slightly different for the T1 and T2 it might explain the larger linear coefficient α. However, it is not easy to correct for such effects since both the porosity n and the specific surface area S will be reduced from such a merging, but with variable magnitude depending on the size of the melted section. The contraction ratio a will also be affected and increased compared to the ideal cubic case depending on the size of the melted section. However, since the pore shapes are not described in this manner in their paper it is not attempted to perform such corrections here.
The dissipation coefficient C KL is here estimated from Eq. (31) based on the assumption that the area ratio is equal to the ideal cubic ratio (a = 0.215). In Table A1 it is shown that the samples T2, T3, and T4 attain coefficient β that are close to 0.6 or slightly higher, which agrees well with the plotted figures and graphs of Huang et al. (2013). However, a slight increase in the coefficients β of T1 further suggest that the shape of the pores in this samples, more so than T2, is slightly altered compared to the T3 and T4 data.
The C KL of the T1-T4 samples attain numerical values C KL = 0.033-0.048, which is in relatively close agreement with the A3 sample with C KL = 0.03 of Gjengedal et al. (2020), but somewhat larger in magnitude, especially for T3 and T4.