Using the borehole permeameter to estimate saturated hydraulic conductivity for glacially over-consolidated soils

The borehole permeameter (BP) method was developed in the 1950s by the United States Bureau of Reclamation to estimate saturated soil hydraulic conductivity (KS) in shallow boreholes completed above the water table. The approach has been improved over the years, and now accounts for flow due to pressure, gravity and soil capillarity. However, the BP method is calibrated only for normally consolidated soils and ponding depth (H) versus borehole radius (r) ratios (H/r) ≤ 22. The primary objective of this study was to recalibrate the BP method for use in glacially over-consolidated soils with H/r ranging from 0.05 to 200. Recalibration consisted of using numerically simulated steady BP flow for five representative glacially over-consolidated soils to update the BP shape function fitting parameters (Z1, Z2, Z3) for nine specified KS values and 15 test pit and borehole configurations. Four sets of fitting parameters were determined, which apply for H/r ≤ 20, H/r ≥ 20, soil with <12% silt content, and soil with >12% silt content. Relative to specified KS, the updated shape function parameters yielded BP estimates of KS with a maximum error of 13% and an average error of 3%, whereas the original shape function parameters (developed for normally consolidated soils and H/r ≤ 22) produced a maximum KS error of 94% and an average error of 23%. The numerical simulations were also used to develop criteria for estimating time required to achieve steady BP flow, and for correcting BP estimates of KS where steady flow was not achieved.


Introduction
Structures for infiltrating stormwater runoff are now common in areas with declining groundwater resources or excessive surface erosion. Infiltration of stormwater is also a key aspect of low impact development (LID) and green stormwater infrastructure (GSI), which are now desired and/or required by many jurisdictions throughout the world. Structures for infiltrating stormwater include ponds or basins, gravel-filled trenches, bioretention swales, drywells or "soakaways", subsurface "leach fields", and infiltration wells or "underground injection control" (UIC) wells. In most cases, these facilities are installed above the water table, and thereby infiltrate into unsaturated soil, i.e. the soil "vadose" zone.
Historically, infiltration facilities were installed in highly permeable soils and sized according to past experience in similar soils (the empirical "trial-and-error" approach). Due to occasional failures, however, many jurisdictions now require site investigations and formal sizing calculations to obtain more reliable estimates of design infiltration rate and capacity.
A variety of approaches have been used to estimate the design infiltration rate and capacity for infiltration facilities including: "guesstimates" based on soil texture, empirical equations based on grain-size distributions, small-scale single-ring and double-ring infiltrometer tests, and larger-scale infiltration tests such as the pilot infiltration test (i.e. "PIT" test, WSDOE 2014). Some ring infiltrometer and most traditional infiltration tests (e.g. pits, drywells, etc.) assume only vertical water flow due to gravity, and do not account for flow due to hydrostatic pressure, horizontal flow, or flow due to the capillary suction (capillarity) of the unsaturated soil surrounding the test facility. Since pressure flow, horizontal flow and capillarity are often a large percentage of total flow, some traditional infiltration tests do not accurately represent the soil's actual infiltration rate and capacity.
A more accurate and reliable approach is to determine infiltration rate and capacity using measurements of soil saturated hydraulic conductivity (K s ) obtained from test methods that formally account for both flow directions (vertical, horizontal), and all three components of soil-water flow (pressure, gravity, capillarity). The constant head well or borehole permeameter (BP) is one such test which is simple to use and well-suited for site investigations (Reynolds et al. 1983Philip 1985;Stephens et al. 1987;Elrick et al. 1989;Reynolds 2008). However, current applications of the BP method focus on ponded head (H) versus borehole radius (r) ratios (H/r) ≤ 22; and it uses a quasi-empirical shape function, C(H/r), which has been calibrated only for normally consolidated soils (Zhang et al. 1998;Reynolds 2008).
In parts of the world that were covered by glaciers during the last ice age, stormwater infiltration is practiced in very dense glacial soils that were over-consolidated by hundreds of meters of glacial ice. These soils have smaller porosity and K S values than normally consolidated soils with similar grainsize distributions. In addition, stormwater infiltration into drilled wells can have H/r ratios as high as 200. BP tests were simulated numerically to determine if the original C(H/r) functions developed by Zhang et al. (1998) were sufficiently accurate for a wide range of H/r ratios (0.05-200) in glacially over-consolidated soils. It was found that the original C(H/r) functions produced a maximum K S error of 94% and an average error of 23% and that K S was over-estimated in 77% of the 150 test simulations. Hence, the C(H/r) functions presented by Zhang et al. (1998) can cause substantial error when applied to large H/r ratios in glacially over-consolidated soils.
The primary objective of this study was consequently to develop revised C(H/r) shape function fitting parameters that would improve the accuracy of the BP method for a broad range of test configurations (i.e. H/r ratios between 0.05 and 200) in glacially over-consolidated soils. Secondary objectives were to use numerical simulations of BP flow for a range of soil types and H/r ratios to: (1) demonstrate the relative importance of the three BP flow components (pressure, gravity, capillarity), (2) estimate the time required to approximate steady-state BP flow, and (3) determine the sensitivity of the BP analysis to background soil matric suction, the soil unsaturated water content function, and the soil unsaturated hydraulic conductivity function.

Borehole permeameter equation
Considerable research has been conducted regarding analytical methods for estimating K S from borehole infiltration tests in the unsaturated zone (see e.g. Reynolds 2010Reynolds , 2011Reynolds , 2013Reynolds , 2015. These methods generally assume a flat-bottom cylindrical test facility (e.g. borehole or pit excavation), isotropic and homogeneous soil, and no water table effects. Work in the 1950's by the United States Bureau of Reclamation provided the Glover analysis (Zangar 1953): where: which for H/r > 10 reduces to: In Eqs.
(1)-(3), K S is soil field-saturated hydraulic conductivity (m/day), Q is the steady-state flow rate out of the test facility (m 3 /day), H is the steady head (ponded water depth) in the test facility (m), r is the radius of the test facility (m), and C g is the Glover C(H/r) shape function (dimensionless). These equations consider only flow caused by the hydrostatic pressure of the ponded water, and do not account for gravitational flow or flow caused by the capillary suction (capillarity) of the surrounding unsaturated soil. In addition, the Glover shape function, C g (Eqs. 2 and 3), considers only test configuration (H/r ratio) and does not account for variations in soil texture or structure. Using numerical simulations of steady BP flow, Reynolds (2013) found that Glover estimates of K S (i.e. Eqs. 1-3) could be in error by as much as 25-1,800% when soil capillarity was strong to very strong, and H/r ratios were moderate to small (i.e. H/r ≤ 10). Stephens and Neuman (1982) recognized the importance of capillarity, and developed quasi-empirical regression relationships that accounted for both pressure flow and capillarity flow: where: C u is the Stephens and Neuman (1982) C(H/r) shape function (dimensionless), and α 1 is a soil capillarity factor (m −1 ). The impact of soil capillarity on steady BP flow rate was found (via numerical simulations) to be a function of soil texture and H/r ratio. Stephens et al. (1987) provided improved regression solutions that used more convenient capillarity parameters, and also extended the analysis to a broader range of soil types. ,  and Philip (1985) developed approximate analytical BP equations that formally account for pressure, gravity and capillarity flow. The Reynolds analysis, which has been tested extensively over the years, has the form: α* is the soil sorptive number (m −1 ), C is the BP shape function (dimensionless), and Z 1 , Z 2 , and Z 3 are the shape function fitting parameters (dimensionless). Equation (6) assumes that H is less than the uncased or screened portion of the test facility, while other constant head and falling head analyses assume that H is greater than the uncased or screened portion of the test facility (see e.g. Reynolds 2010Reynolds , 2011. The three terms in the denominator of Eq. (6) account, respectively, for flow through the wall and base of the test facility due to the hydrostatic pressure of the ponded water, gravity flow through the base of the test facility, and capillarity flow through the wall and base of the test facility due to the surrounding unsaturated porous material. Flow due to hydrostatic pressure accounts for most of the flow out of the test facility when H >> r, while gravity flow and capillarity flow often dominate when H < r (Reynolds 2008;Elrick and Reynolds 1992). The relative importance of the three flow components is discussed in section 'Relative importance of pressure flow, gravity flow and capillarity flow'.
As described in  and Elrick et al. (1989), soil capillarity can be parameterized using the "alpha" relationship, α (m −1 ), which is defined by: where the matric flux potential, ∅ m (m 2 /day), is given by: K(ψ) is soil hydraulic conductivity as a function of soil matric suction ψ (m), and ψ i (m) is the antecedent matric suction in the background soil at the time of the test. (Note that matric suction is defined here as the negative of soil pore water pressure head, so ψ is positive when the soil is unsaturated). In effect, ∅ m equals the area under the K(ψ) curve between the matric suction in the background soil (ψ = ψ i ) and the matric suction at the leading edge of the saturated bulb surrounding the BP injection zone (ψ = 0; . This means that ∅ m is a maximum when the background soil is dry (ψ i large), small when the background soil is wet (ψ i close to zero), and zero when the background soil is saturated (ψ i = 0).
When the unsaturated soil is at field capacity or drier, K(ψ i ) << K S ; and hence, α in Eq. (8) can often be simplified to: where α* (m −1 ) is known as the "sorptive number". Field studies have shown that α* is relatively constant for a broad range of porous materials and can therefore be estimated using a lookup table based on soil texture and structure (Table 1; Elrick et al. 1989;Reynolds 2008Reynolds , 2013. "Structured soil" (Table 1) refers to soil with cracks and/or biopores (e.g., root holes) that can increase bulk soil α* and K S . Generally speaking, target soils for stormwater infiltration are well below the plant root zone, and therefore have few cracks and biopores. Although not explicitly stated by Elrick et al. (1989) and Reynolds (2008Reynolds ( , 2013, the α* values and shape function parameters (Z 1 , Z 2 , Z 3 ) in Table 1 apply for near-surface, normally consolidated soils and H/r ratios ≤ 22.

Borehole permeameter calibration
The calibration procedure involved calculating K S via the BP equation (Eq. 6) using steady flow rate values (Q) generated by numerical simulations of BP flow for 150 test scenarios. The test scenarios included all combinations of five "representative" glacially over-consolidated soils, two plausible K S values for each soil type, and 15 BP test configurations where H/r ratio varied from 0.05 to 200. Calibration was conducted using the Solver optimization algorithm in Excel, which Table 1 Sorptive number (α*) and shape function (C) parameters (Z 1 , Z 2 , Z 3 ) for a range of normally consolidated soils where ratio of steady ponding depth (H) to borehole radius (r) ≤ 20 (Adapted from Elrick et al. 1989;Zhang et al. 1998) changes user-selected variables until a specified objective is minimized, maximized or becomes equal to a specified value. In this study, the C(H/r) shape function fitting parameters (Z 1 , Z 2 , and Z 3 in Eq. 7) were varied by Solver (using the generalized nonlinear reduced gradient method) until the maximum error between the BP-calculated K S values and the specified K S values was minimized for the 150 test scenarios. Borehole permeameter flow was estimated using SEEP/W, a finite element numerical model that can simulate multidimensional and axisymmetric flow in saturated and unsaturated porous media (GEOSLOPE International Ltd., Calgary, Alberta, Canada). Unsaturated flow simulation requires specifying the unsaturated volumetric soil-water content function θ(ψ) (soil-water content as a function of soil matric suction) and the unsaturated hydraulic conductivity function K(ψ) (hydraulic conductivity as a function of soil matric suction). These functions are described in section 'Soil-water content function θ(ψ) and unsaturated hydraulic conductivity function K(ψ)'.

Model domains and test configurations
The SEEP/W numerical flow domains for the three test configurations are shown in Fig. 1 and summarized in Table 2. The simulations assumed axisymmetric flow, with no-flow boundaries along the top and outside radius of the flow domain, and a unit hydraulic head gradient boundary at the bottom of the flow domain. The test facilities were cylindrical test pits or boreholes with constant hydraulic head boundaries along the base and submerged portion of the test facility wall. The facility radius (r), constant ponded head (H), and H/r ratios for each of the 15 test scenarios are provided in Table 3. The simulations used graded meshes of rectangular and triangular elements. For pit simulations, element size increased in steps from 0.05 m × 0.05 m along the pit wall to 0.1 m × 0.1 m at distance; and for borehole simulations, element size increased in steps from 0.025 m × 0.025 m along the borehole wall to 0.1 m × 0.1 m at distance (Fig. 1). Test simulations showed that larger flow domains and finer element sizes had minimal impact on Q and K S (data not shown).
The numerical simulations calculated water flow rate or discharge, Q (m 3 /day), versus time, t (days), out of the test facilities over a 24-h period. As transient flow was simulated, true steady flow rate (constant Q) was usually closely approached but not truly achieved. True steady flow requires steady-state simulations that are often impractical because they require very large flow domains to avoid external boundary effects, as well as very large run times to achieve convergence. As a result, longer time transient flow simulations generally result in slightly lower Q values, which in turn result in slightly lower K S estimates by the BP approach. The Appendix illustrates the time required to obtain Q values that approximate steady state for the different soil types and test configurations. The "fixed head" boundary condition applies to the base and submerged portion of the test facility wall, and it refers to specified hydraulic head that is constant in space and time

Representative soil types
Five "representative" glacially over-consolidated soils (based on data from the Puget Sound region of Washington State) were defined for this analysis, including four derived from glacial advance outwash materials (silty Qva, fine Qva, finemedium Qva, and fine-coarse Qva) and one glacial till (Qvt). Glacially over-consolidated soils are common throughout the world, and they are often used for stormwater infiltration. Figure 2 shows grain size distributions and Table 4 summarizes key parameters for the representative soils.
The four Qva soils represent materials that were deposited by streams in front of an advancing glacier, and then overrun and consolidated by the glacier. The silt content (particles passing through a 0.075-mm sieve) ranges from 3 to 17%. Using the Unified Soil Classification System (USCS, per ASTM D 2487), the four soils are classified as SP (poorly graded sand with less than 5% silt), SW (well-graded sand with less than 5% silt), SP-SM (sand with 5-12% silt), or SM (sand with greater than 12% silt). The silty Qva, fine Qva and fine-medium Qva soils are relatively poorly graded, while the fine-coarse Qva soil is well graded ( Fig. 2; Table 4). The Qvt represents glacial till soils which typically have minimal sorting and little or no layering. Qvt is well graded (Fig. 2) with 20% silt and a USCS classification of SM.
Because of glacial over-consolidation, Qva soils range from medium dense to very dense, with standard penetration test (SPT, ASTM D1586 / D1586M-18) blow counts typically ranging from 20 to 80 blows per 30 cm. Qvt is usually very dense with SPT blow counts greater than 100 blows per 30 cm. The SPT is a standard geotechnical method used to document relative soil density during drilling and soil sampling. The SPT is conducted by driving a 51-mm-diameter sampler using a 63.5-kg slide-hammer dropped from a height of 76 cm.
Two "typical" K S values were assigned to each soil type (Table 4), based on the aforementioned soil characteristics and hundreds of K S field measurements in the Puget Sound area. The remaining parameters in Table 4 (porosity, liquid limit, residual soil-water content, background soil matric suction, α′, and n) are inputs for either the modified Kovacs θ(ψ) function (Aubertin et al. 2003) or the van Genuchten (1980) θ(ψ) function (section 'Soil-water content function θ(ψ) and   unsaturated hydraulic conductivity function K(ψ)'). Most of the θ(ψ) input parameters are difficult to measure in situ and not readily available in the literature for glacially overconsolidated soils; hence, "representative" values are given in Table 4. Assumed porosities θ S of the five soils were 17% for Qvt, 25% for silty Qva, and 30% for the fine, fine-medium and fine-coarse Qva soils. These porosities are less than typically measured for normally consolidated soils and are intended to reflect the effects of extreme glacial compaction. The Atterberg liquid limit (gravimetric water content at which soil behavior transitions from plastic to liquid, ASTM method D 4318) is controlled primarily by the soil clay/silt ratio on a weight percent basis. Advance outwash and glacial till soils usually contain very little clay, hence an Atterberg liquid limit of zero was assigned to all five soils. Residual soil-water content θ r (α′) and n values were estimated during development of the θ(ψ) functions, as described in the following section. The background soil matric suction ψ i was estimated from the θ(ψ) functions using an assumed background soil-water content of 10% (dry soil).
Soil-water content function θ(ψ) and unsaturated hydraulic conductivity function K(ψ) The unsaturated volumetric soil-water content function θ(ψ) was initially defined using the modified Kovacs model (Aubertin et al. 2003). This θ(ψ) model uses soil porosity, grain diameters representing 10% passing and 60% passing on the grain size distribution curve, and the Atterberg liquid limit. Table 4 provides the soil parameters used in the modified Kovacs model and Fig. 3 shows soil-water content versus matric suction (negative pore pressure) for each of the five soil types using the modified Kovacs model. During sensitivity analyses it was observed that the numerical model results were unrealistic for certain scenarios and it was determined that de-coupling of the soil-water content function θ(ψ) (based on the Modified Kovac's model) and the unsaturated hydraulic conductivity function K(ψ) (based on van Genuchten) was the likely cause. Therefore, van Genuchten θ(ψ) curves, which are coupled to the van Genuchten K(ψ) curves, were fit to the modified Kovacs θ(ψ) curves as shown on Fig. 3. The van Genuchten θ(ψ) equation is: where θ S is porosity, θ r is residual soil-water content, α′ (m −1 ), n (−) and m (−) are empirical fitting parameters, and the Mualem (1976) pore continuity model is assumed, i.e. m = 1 -(1/n). Fitting of the van Genuchten θ(ψ) curves to the Modified Kovac's curves assumed that θ S was fixed while θ r , α′, and n were allowed to vary. The fitted parameters are provided in Table 4 and the results of the curve fitting are shown in Fig. 3. Note that the van Genuchten θ(ψ) curves tend to systematically overestimate the modified Kovacs curves in the dry soil end (ψ > 10 m, Fig. 3); however, this is well beyond the maximum ψ (ψ i ) used for the simulations (Tables 4 and 8). The SEEPW simulations were performed using the van Genuchten θ(ψ) curves.
The unsaturated hydraulic conductivity function K(ψ) was defined using the van Genuchten (1980) equation: where K S (m day −1 ) is saturated hydraulic conductivity, and α′ (m −1 ), n (−) and m (−) are the same fitting parameters used in Table 4 Characteristics of representative glacially over-consolidated soil types and baseline SEEP/W parameters used in the soil-water content and hydraulic conductivity models. Qva is advance outwash, Qvt is glacial till, D 60 and D 10 are grain diameters corresponding, respectively, to 60% passing and 10% passing on the grain size distribution curve, and USCS is Unified Soil Classification System. Background soil matric suction is based on an assumed background volumetric soil-water content of 10%. The van Genuchten fitting parameters apply for Eqs. (5) and (6) permeameter Eq. (11). The SEEPW model develops the van Genuchten K(ψ) curves based on K S and the van Genuchten θ(ψ) curves. As shown in Fig. 4, K(ψ) for fine-coarse Qva begins to decrease at a much lower matric suction than for the finergrained soils.

Sorptive number (α*) calculations
The sorptive number (α*) can be estimated for the five soil types using Eq. (10) and the K(ψ) functions generated by SEEP/W (Table 4; Fig. 4). The matric flux potential (∅ m ) was calculated by numerically integrating under the K(ψ) curves, and α* was calculated for background soilwater content (θ b ) ranging from dry soil to virtual saturation. Numerical integration used the trapezoidal rule and the matric suction/water content intervals generated by SEEP/W. Plots of α* versus background soil-water content, θ b (Fig. 5), reveals that α* is relatively constant in dry and moist soil, but increases dramatically as the background soil approaches saturation. The sudden and rapid increase occurs because the highly non-linear ∅ m relationship suddenly decreases towards zero at near-saturation, signifying that capillarity is negligible in all near-saturated soils (and zero in all saturated soils). The soil α* values for dry/moist soil increased from 1.17 m −1 in Qvt, to 1.33 m −1 in silty Qva, to 2.5 m −1 in fine Qva, to 3.9 m −1 in fine-medium Qva, and to 25 m −1 in finecoarse Qva (Table 5), which reflects the fact that soil capillarity is often substantial in silty soils (e.g. Qvt, silty Qva), but can decrease significantly with decreasing silt content (e.g. fine-coarse Qva). Fig. 4 Example unsaturated hydraulic conductivity curves K(ψ) used for SEEP/W simulation of BP flow in four advance outwash soils (Qva) and one glacial till soil (Qvt) ( Table 4) Legend Fig. 3 Volumetric soil-water content curves used for SEEP/W simulation of BP flow in four advance outwash soils (Qva) and one glacial till soil (Qvt) ( Table 4). This figure illustrates the match between the modified Kovacs model (solid blue line) and the van Genuchten model (dashed orange line)

Sensitivity analyses
Calibration of the shape function C(H/r) curves was conducted using deterministic values of the soil parameters (θ S , θ r , α′, and n) used to develop the θ(ψ) and K(ψ) curves. Sensitivity of the calibrated C(H/r) fitting parameters to variations in the θ(ψ) and K(ψ) curves was tested by modifying the underlying parameters used to develop the θ(ψ) and K(ψ) curves (θ S , θ r , α′, and n). Numerical sensitivity runs were performed for a single test configuration: the shallow borehole with H = 2 m and r = 0.25 m, and the following changes to the θ(ψ) and K(ψ) parameters: porosity θ S was increased by 5%; θ r and α′ were increased by 50%, and n was decreased by 0.5. Calibration of the shape function C(H/r) curves assumes a constant value for α* based on a background soil-water content θ b of 10%. As shown in Fig. 5, α* was constant for all soil types at this soil-water content and did not change for dryer soil conditions. However, as discussed in section 'Sorptive number (α*) calculations', capillarity flow does decrease as soil-water content approaches full saturation and BP flow is therefore expected to decrease as well. Sensitivity of the numerical flow results to θ b is evaluated for a single test configuration (the shallow borehole with H = 2 m and r = 0.25 m) using wetter θ b values, as discussed in section 'Numerical model sensitivity analysis'.
Only one parameter was changed for each sensitivity run and the percent difference between baseline Q and revised Q was calculated using:

Numerical simulations
SEEP/W simulations of BP flow were conducted for five soil types, two K S values for each soil type, and 15 test configurations, for a total of 150 simulations (Tables 3 and 4). The simulations were run for 24 h, except for the fine-medium Qva and fine-coarse Qva soils, which were terminated after 6 h because the wetted zone started to impinge on flow domain boundaries. As discussed in the Appendix, 6 h was still sufficient to achieve approximate steady-state flow in these coarse textured soils. Zero matric suction and water content contours are provided in Fig. 6 for each soil type and one test configuration (H = 2 m, r = 0.25 m, H/r = 8) after 6 h of flow. As shown in the figure, zero matric suction (dashed blue contour line) extends deeper below the borehole as K S and α* increase. Borehole flux reached the bottom flow domain boundary for finemedium Qva and fine-coarse Qva, which are the soil types with the largest K S and α* values. Because unit hydraulic gradient was specified on this boundary, presence of borehole Legend Fig. 5 Calculated values of α* as a function of background soil-water content (θ b ) for the four advance outwash soils (Qva) and one glacial till soil (Qvt) ( Table 4). Note that α* is relatively constant until the soil approaches full saturation for all the soils except fine-coarse Qva, which begins to increase at a background soil-water content of 15% Table 5 Sorptive Number (α*) and C(H/r) shape function (Eq. 7) parameters (Z 1 , Z 2 , Z 3 ) for five representative glacially over-consolidated soils, including four types of advance outwash (Qva) and one glacial till (Qvt). Different shape function parameters are developed for test configurations where ponded head (H) to radius (r) ratio was H/r ≤ 20 or H/r ≥ 20, and for soils with > 12% silt (USCS soil type SM) or < 12% silt (USCS soil types SP, SW, or SP-SM) Soil type α* (m −1 ) Low ponded head (H/r ≤ 20) High ponded head (H/r > 20) flux does not affect simulated flow appreciably as long as the boundary is not contacted by the zero matric suction contour.
The Appendix provides an analysis of the approximate simulation time required to approximate steady-state BP flow for the different soil types and test configurations. "Approximate steady state" was defined as the time required to achieve a BP flow rate that was within 5% of the BP flow rate after 24 h. As discussed in the Appendix, approximate steady-state BP flow was achieved within 0.5−4 h for the fine-medium Qva and finecoarse Qva; within 3-6.5 h for fine Qva; and within 11-17 h for the fine-grained soils (Qvt and silty Qva). This information could be used to estimate suitable/likely test durations for particular BP test configurations and soil types.

Calibrations
As discussed in section 'Borehole permeameter calibration', a spreadsheet was set up to estimate the K S specified in the SEEP/W simulations using the BP equation (Eq. 6) and recalibrated C(H/r) shape function fitting parameters (Z 1 , Z 2 , Z 3 in Eq. 7). The calibration process was designed to minimize the maximum individual error in the BP estimate of K S across the 150 test scenarios. The maximum individual K S error was minimized (instead of minimizing the average K S error) to ensure that the BP analysis always met or exceeded a known degree of accuracy.
Initial recalibration tests indicated that a single set of C(H/ r) shape function fitting parameters would not provide sufficiently accurate K S determinations across the full range of soil types and test configurations. The next phase of calibration therefore extended the approach of Reynolds et al. (1983, and assumed that the shape function parameters (Z 1 , Z 2 , Z 3 ) depended on both soil type (i.e. α* value) and H/r ratio. This resulted in four sets of C(H/r) shape function parameters that applied separately for fine-grained soil (>12% silt), coarse-grained soil (<12% silt), small H/r ratio (≤ 20), and large H/r ratio (≥20). These parameters provided BP estimates of K S with a maximum error of 13% and an average error of 3%, relative to the K S specified in the numerical simulations. Maximum error for each soil type and H/r range are summarized in Table 6.
Plots of the four recalibrated C(H/r) shape functions are given in Fig. 7, and the corresponding Z 1 , Z 2 and Z 3 fitting parameters are given in Table 5. Figure 7a illustrates relatively minor differences in the shape functions when H/r ≤ 20, while Fig. 7b illustrates substantial differences when H/r > 20. Figure 7 also shows the Zhang et al. (1998) C(H/r) shape function for α* = 4 m −1 , which is seen to overestimate the shape functions for over-consolidated soils when H/r < 30 and underestimate the shape functions for over-consolidated soils when H/r > 30.

Relative importance of pressure flow, gravity flow and capillarity flow
The BP flow equation can be re-arranged to: where Q T is total flow from the test facility (borehole or pit), Q P refers to the first term on the left and represents threedimensional (3D) pressure flow through the sides and base of the test facility due to the hydrostatic pressure of the ponded water; Q G refers to the second term on the left and represents vertical gravity flow through the test facility base; and the third Elevation (m)

Distance (m)
Dashed line indicates the zero-matric suction contour (ψ = 0). Qvt -Glacial till Qva -Advance Outwash Legend Fig. 6 Zero matric suction and water content contours after 6 h of SEEP/W simulated flow from the BP into four advance outwash soils (Qva) and one glacial till soil (Qvt) ( Table 4). Borehole configuration was H = 2 m and r = 0.25 m, H/r = 8 term on the left (Q C ) represents 3D capillarity flow through the facility walls and base due to the capillary suction of the background unsaturated soil. The relative importance of Q P , Q G and Q C varies with soil type and H/r ratio and is illustrated in Fig. 8 using the five representative soils (Table 4) and recalibrated shape function fitting parameters (Table 5).
For all soils, pressure flow (Q P ) increases as H/r ratio increases, accounting for less than 30% of total flow (Q T ) when H/r is less than 0.1, but at least 90% of Q T when H/r exceeds 10 (Fig. 8). For H/r less than 1.0, capillarity flow (Q C ) is relatively important for fine-grained soils (Qvt and silty Qva) and less important for coarse-grained soils (fine-coarse Qva). When H/r is less than 0.1, gravity flow (Q G ) dominates in coarse-grained soils such as fine-coarse Qva. The maximum error in K S estimates (Table 6) suggests that the BP method tends to be more accurate for larger H/r ratios (i.e. ≥ 20) and soils with lower silt contents (i.e. α* > 1.33 m −1 ).

Numerical model sensitivity analysis
The five "representative" soil types were defined using parameters with significant variability and uncertainty. This section examines how modifying some parameters impacts the numerically simulated Q values. The sensitivity analyses were conducted for all five soil types and a single test configuration (H = 2 m, r = 0.25 m, H/r = 8). Baseline and revised Q values were obtained using the baseline and revised parameters provided in Tables 7 and 8. As evident in Eqs. (6) and (14), K S is linearly related to Q so variation in Q produces the same variation in K S .
As discussed in section 'Soil-water content function θ(ψ) and unsaturated hydraulic conductivity function K(ψ)', the van Genuchten (1980) equations were used to define the K(ψ) and θ(ψ) functions. As discussed in section 'Sensitivity analyses', sensitivity to changes in these functions was tested by changing the van Genuchten (1980) parameters (θ S , θ r , α′, and n) as shown on Table 7. Comparison of the revised Q values with the baseline Q values indicate that the results were virtually unchanged for all the soil types except Qvt, which changed by ±4% or less (Table 7). SEEP/W simulation of BP flow requires an activation pressure, which is the specified background matric suction ψ i , and corresponding background soil-water content θ b of the soil surrounding the test facility. As shown on Table 8, sensitivity to θ b was determined by specifying a lower ψ i (and thereby wetter θ b ). Increasing background soil-water content resulted in less than 2% reduction in steady Q for the sandy soils (fine Qva, fine-medium Qva, fine-coarse Qva), an 11% reduction for silty Qva, and a 16% reduction for Qvt (Table 8). The greater impact for silty soils occurred because: (1) the capillarity of the Legend Fig. 7 Recalibrated shape functions (C) for soils with <12% silt (blue lines), soils with >12% silt (red lines), H/r < 20 (solid lines), and H/r > 20 (dashed lines). Panel a gives shape functions for H/r between 0 and 20; panel b gives shape functions for H/r between 0 and 200. Zhang et al. (1998) shape function for α* = 4.1 m −1 is provided for comparison. H is borehole ponding depth, r is test facility (pit or borehole) radius  Table 7 Comparison of SEEP/W-simulated steady flow rates Q for different values of parameters used by the van Genuchten (1980) θ(ψ) and K(ψ) models, including porosity θ S , residual soil-water content θ r , and the fitting parameters α′. and n. Baseline and revised steady flow rates are shown for four types of advance outwash soil (Qva) and one glacial till soil (Qvt background soil (as represented by ∅ m , Eq. 9) decreases as ψ i decreases and θ b increases; and (2) sensitivity to capillarity increases with increasing silt content (Fig. 8). In theory, this sensitivity could be mitigated by measuring θ b in the field and recalculating α*. Field measurements of θ b are rare, however.

Uncertainty in the BP analysis
The recalibrated BP analysis provides (as discussed in section 'Calibrations') K S estimates with maximum error of 13% and average error of 3% for a selected range of representative soils (Table 4) and H/r ratios (Table 3). Figure 9 shows the 95% confidence limits for K S for each shape function. As shown on the figure, the 95% confidence bands are less than 6.5% for all H/r ratios and all soil types except for soils with silt >12% and H/r ratios <1. For these silty soils and low H/r ratios, the 95% confidence bands are as high as 12%. Although these errors might be reduced by defining more than four shape functions (Table 5), this is likely not warranted as actual field soils have unknown and uncontrollable degrees of heterogeneity and anisotropy. The sorptive number, α* (m −1 ), represents the capillarity of the soil, and is determined for each soil using Eq. (10).  However, the parameters required to solve Eq. (10) will generally not be available when testing actual field soils, as standard geotechnical exploration methods usually provide only grainsize distribution, soil-water content, and a qualitative measure of soil density. Therefore, α* must be estimated using limited soil texture and structure information such as provided in Table 1 for normally consolidated soils, and in Table 5 for glacially over-consolidated soils. The relative proportion of capillarity flow usually increases with increasing fines (e.g. increasing silt content), and this is represented in the BP analysis (Eq. 14) by smaller α* values (i.e. capillarity flow increases as α* decreases). Hence, uncertainty in α* has proportionately greater impact on K S accuracy in fine soils than in coarse soils. This can always be mitigated, however, by increasing the H/r ratio, which decreases the proportion of capillarity flow relative to pressure flow (Fig.  8, see also Reynolds 2008).

Limitations of the BP method
As discussed in Reynolds (2008), Archer et al. (2014) and others, the BP approach does not account for the following factors: & Entrapped or encapsulated air. Rapidly infiltrating water from the BP test facility can entrap/encapsulate air in soil pores, which may decrease flow and reduce the effective K S value. This air often dissolves gradually, resulting in an increase in effective K S over time. & Siltation and/or drill-induced smearing and compaction along the test facility wall and base. Infiltration surfaces that are substantially silted, smeared or compacted have lower K S values than the background soil. & Proximity of a water table. A regional or perched water table that intersects or occurs just below a test facility can reduce hydraulic gradients, producing flow and K S estimates that are artificially low. & Heterogeneity and/or anisotropy. In heterogeneous soils, the BP method provides a bulk average K S of the soil volume wetted by the test facility (Fig. 6). In materials with vertical-horizontal anisotropy, the BP tends to yield a K S value that falls between the vertical K S and the horizontal K S (e.g. .
As one or more of these factors are likely present in all borehole infiltration tests, K S estimates from field testing will inevitably reflect the aggregate effects of manmade and natural porous medium conditions within the wetted soil surrounding the pit/borehole injection zone. One might argue, however, that the same porous medium conditions will also exist in production-scale stormwater infiltration facilities; and hence K S estimates perturbed by the preceding factors are appropriate, and perhaps even preferred, for feasibility assessments and facility design.

Applicability of glacially over-consolidated BP fitting parameters
The C(H/r) shape function parameters in Table 5 were developed using soils in the Puget Sound region (as described in section 'Representative soil types'), but are likely also valid in other parts of the world for glacially over-consolidated soils with similar grain size distributions. Glacially overconsolidated soils are also found across most of Canada, southern Alaska, the mid-western and northeast portions of the United States, Scandinavia, the northern portions of the British Isles, the northern portions of eastern Europe, portions of Russia, and within the world's major mountain ranges (Ehlers et al. 2011). The shape function parameters in Table 5 require further validation, however, before use in normally consolidated soils, or in soils with different structure or grainsize distributions. The shape function parameters developed by Zhang et al. (1998) (Table 1) are advisable for structured and normally consolidated soils when H/r ≤ 22.

Conclusions
Numerical simulations of BP tests were used to develop four sets of recalibrated C(H/r) shape function fitting parameters (Z 1 , Z 2 , Z 3 , Eq. 7; Table 5) for use in glacially overconsolidated advance outwash and till materials. The parameters were developed for H/r ratios between 0.05 and 200 and apply to both excavated pits and boreholes completed above the water table. The parameters provided BP estimates of soil K S with a maximum error of 13% and an average error of 3%, which is more than accurate enough for feasibility assessment and design of stormwater infiltration facilities. The original C(H/r) fitting parameters which were developed for normally consolidated soils and H/r ≤ 20 caused unacceptably large K S errors up to 94%. The numerical simulations were also used to develop criteria for estimating time required to achieve steady BP flow, and for correcting K S estimates when steady BP flow was not achieved (Appendix).
It was concluded that the BP method using recalibrated C(H/r) shape functions is suitable for estimating K S in glaciated soils within the Puget Sound region of western Washington State (United States), and in other parts of the world with similar soils. The calibration approach developed here could be used to develop BP shape functions for other soil types and specialized borehole configurations. early stages of this work from Salini Sasidharan, Agricultural Research Service (USDA), and Daniel B. Stephens, Daniel B. Stephens & Associates Inc., are greatly appreciated. The authors also thank associate editor Andrew Frampton, reviewer Robson A. Armindo, and an anonymous reviewer for their detailed and constructive comments on the initial manuscript.
Funding information This work was supported in part by the Science and Technology Branch of Agriculture and Agri-Food Canada (Study Number 2380).
Appendix: approximation of steady-state borehole flow rate The borehole permeameter equation (Eq. 6) assumes steady-state flow rate (Q) within an infinite flow domain. This is difficult to simulate numerically, as it requires very long run times to achieve numerical convergence of the steady-state flow equation, and very large numbers of nodes to place the radial and bottom flow domain boundaries at "numerical infinity". It is possible, however, to conduct transient flow simulations in much smaller flow domains (e.g. Fig. 1) where the radial and bottom boundaries are just far enough way to allow near-steady BP flow before boundary effects occur (e.g. Fig. 6). Therefore, Q vs. t was simulated for test pit, shallow borehole and deep borehole configurations, and Q after 't' hours (Q t ) was compared with Q after 24 h (Q 24 ). The simulations (detailed in the following) indicated that flow was effectively steady after 24 h in most scenarios, and Q t /Q 24 ≤ 1.05 was considered indicative of effective steady flow after 't' h of infiltration. Determining if and where Q 6 /Q 24 ≤ 1.05 occurs was of particular interest because most field infiltration tests are terminated after about 6 h (due to cost).
As in actual field tests, simulated BP flow rate (Q) decreased with time to become effectively constant (Fig. 10). Only two scenarios (deep borehole with H = 20 m, finemedium and fine-coarse Qva) exhibited boundary interference effects, as evidenced by an abrupt change in slope of Q vs. t at about 12 h (Fig. 10b). Generally speaking, the Q 6 /Q 24 ratio decreased and the percentage of silt decreased from Qvt, to silty Qva, to fine Qva, to fine-medium Qva, to fine-coarse Qva ( Fig. 11; Table 9). Q 6 /Q 24 ratios ≤ 1.05 occurred for H < 10 m in fine Qva, and for all test configurations in fine-medium Qva and fine-coarse Qva, but not for the other scenarios (Table 9). Hence, 6 h was not enough time to achieve steady BP flow rate for H ≥ 10 m in fine Qva, and for all configurations in silty Qva and Qvt (Table 9). This information may be useful for correcting K S estimates when BP test duration was not long enough to approximate steady-state flow.
Time to achieve near-steady BP flow rate (i.e. Q t /Q 24 ≤ 1.05) varied substantially among test configurations and soil types, requiring 0.5-4 h in fine-medium Qva and fine-coarse Qva, 3-6.5 h in fine Qva, and 11-17 h in silty Qva and Qvt (Table 10). Planning BP test durations must consequently account for both soil type and test facility configuration, and might be estimated using Table 10. and H/r ratio for four advance outwash soils (Qva) and one glacial till soil (Qvt) ( Table 4). H is steady ponding depth, r is pit/ borehole radius. The lines are linear regressions for each soil type. Note that the Q 6 /Q 24 ratio changes with test type, H/r ratio, and soil type  Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.