FLASH

FLASH (Flow-Log Analysis of Single Holes) is a spreadsheet-based graphical user interface (GUI) for a multi-layer Thiem solution for steady-state flow to a borehole. The software allows for analysis of data from boreholes in which step increases in flow (from discrete fractures or zones) or linear increases in flow (over multiple aquifer layers) are observed. Data are entered in the worksheet 'FIELD DATA.' Interpretations of ambient and stressed flow are enetered along with physical parameters and water-level measurements on the 'INPUTS' sheet. Visual Basic for Applications (VBA) routines provide an interface to the Excel Solver, which can be used for model calibration.

This spreadsheet uses a multi-layer Thiem analytical model to generate a borehole flow log representing the flow that would be measured in a borehole with a specified number of aquifers for a given set of values for the water levels in each of the aquifers, and the difference between the open-hole water level under static and steady pumping (or injection) conditions. The program is used in a forward modeling mode (that is, trial and error) until a satisfactory match between measured and computed flow is achieved for a gvien set of input parameters. In the GUI the input parameters, which generally includes transmissivity and head, can be changed; the simulated flow profiles are updated automatically. The GUI also supports use of the Excel Solver for automated calibration.

The equations we are fitting to find the flows are (Day-Lewis et al., 2011):

\begin{align} Q^a_i = \frac{2\pi T^{factor}_iT^{total}(h^a_w-h^0_i)}{ln(r_o/r_w)} \end{align}

\begin{align} Q^s_i = \frac{2\pi T^{factor}_iT^{total}(h^s_w-h^0_i)}{ln(r_o/r_w)} \end{align}

REQUIRED INPUT

Be careful to use English parameters -- feet, inches, gallons per minute
Elevation of measuring point - elevation used for calculation of water levels based on drawdown and depths entered below, in feet;
Number of flow Zones - number of transmissive zones (layers or fractures) up to a maximum of 10;
Well diameter - diameter of well in the interval over which drawdown occurs, in inches;
Drawdown - the distance between the ambient and the quasi-steady state water level measured under stressed (pumping or injection conditions), in feet; Depth to ambient water level - depth from measurement point to water level in well in absence of stress, in feet;
Depth at bottom of casing - depth at which the well is open to the formation, in feet. The hydraulic solution will be plotted from this depth down to the depth at the bottom of the well.
Depth at bottom of well - depth at the bottom of the well, in feet. The hydraulic solution will be plotted to this depth.
Radius of influence - the radial distance from the well at which all transsmive zones remains at equilibrium condition, in feet;
Total Transmissivity - the transmissivity based on an open-hole pump test, used for a starting value, in feet squared per day;

Import Require Modules



In [1]:

    
# Import the necessary Python modules
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
from math import *
from scipy.optimize import minimize

Initialize Parameters



In [2]:

    
# Initialize paramters
elev = 555.53           # Elevation of the measuring point (ft)
nzones = 3              # Number for flow zones
Dw = 6                  # Well diameter (inches)
drawdown = 3.5          # Drawdown (ft)
Hwa = 4.62              # Depth to ambient water level (ft)
bot_casing = 10         # Depth to bottom of casing (ft)
bot_well = 140          # Depth at bottom of the well (ft)
Ro = 95                 # Radius of influene (ft)
Ttot = 26               # Total Transmissivity (square ft/day)

FRACTURE INFORMATION



In [3]:

    
# The number of fractures can be changes by editing the lists of properties. 
# All lists must be the same length.
depth = np.array([35.0, 45.0, 65.0])                # Fracture depth (ft)
Qamb = np.array([0.0, 0.02, 0.02])                  # Ambient flow (GPM)
Qstress = np.array([0.5, 0.5, 0.23])                # Stressed flow (GPM)
ff_h = np.array([547.41, 550.91, 551.24])           # Farfield head (ft)

FIELD DATA

Flow normalization:

Assumes inflow to borehole is proportional to the total discharge or drawdown
If maximum (reliable flow measurement) undermeasures, the normalization can correct the raw data to reflect the total proportional flow.
If the pump cuts out or if the pumping rate cannot be held constant, normalize the flow at any location and divide by the total amount of flow at that time.

\begin{align} CF = Q_{raw}/Q_{raw-max} \end{align}\begin{align} Q_N = (CF)(Q_{raw}) \end{align}

Quasi-Steady State: Change in borehole storage converted to flow rate, which is less than the lower limit of the flowmeter.

6-inch borehole: $\Delta h$ = 0.01 ft/min corresponds to 0.015 gpm for HPFM (or 0.1 gpm for EMFM)
If $Q_{ss}$ is not reached, the change in storage can be removed from the stressed flow rate. To do this, use a transducer to measure water levels during the flowmeter measurements.

\begin{align} Q_{ss} = \Delta Q_{bs} = \pi r^2 \Delta h \end{align}

Remember to account for borehole storage and normalize the flow rate to the known or measured flow.



In [4]:

    
# Edit lists to add your field data

# Ambient Flow

# Depth (ft)
Za = [21, 33.5, 37.0, 39.1, 50.0, 63.0, 69.5, 89.1, 100.0, 109.0]

# Discharge, Q, (GPM)
Qa = [0.0, 0.0, 0.019, 0.017, 0.02, 0.019, 0, 0, 0, 0]

# Pumped depth (ft)
Zs = [21.3, 35.5, 37, 50, 64, 66, 82, 100, 125]

# Pumped flow
Qs = [0.33, 0.50, 0.49, 0.25, 0.23, 0, 0, 0, 0]

# Check that the parameter lists are the same length (easy mistake to make)
assert len(Za) == len(Qa), "Ambient flow vectors not the same length"
assert len(Zs) == len(Qs), "Pumping flow vectors not the same length"

Forward Model Code



In [5]:

    
# Forward Model
import numpy as np
def thiem(Tfac, Ttot, del_h, Ro, Dw):
    '''This function implements equation 2 and 3 
    based on the Thiem Equation'''
    # Convert all quantities to mks
    del_h = del_h * 0.3048 # Convert ft to m
    Ro = float(Ro)*0.3048
    Rw = float(Dw)/2 # Convert well diameter to radius
    Rw = float(Rw)*(0.3048/12) # Convert inches to m
    Ttot = Ttot*0.3048*0.3048/(3600*24) # Convert to m2/sec
    Q = 2.0*np.pi*Tfac*Ttot*del_h/log(float(Ro)/Rw)
    Q = Q * 1000.0 # Convert cubic meters/sec to liters/sec
    Q = Q * 15.8503230745 # Convert L/s to gpm
    return Q



In [6]:

    
# Code to calculate the cumulative flow
def calculate_flow(Tfac,Ttot,del_h,Ro,Dw,ambient=True):
    # Loop over the number of fractures
    Q_sim = np.zeros(len(Tfac))
    Qfrac = 0.0
    for i in range(len(Tfac)):
        if ambient:
            Qfrac += thiem(Tfac[i],Ttot,del_h[i],Ro,Dw)
        else:
            Qfrac += thiem(Tfac[i],Ttot,del_h[i]+drawdown,Ro,Dw)
        Q_sim[i] = Qfrac
    Q_sim = np.flipud(Q_sim) #Reverse to go from top of well down
    return Q_sim

Model Parameters



In [7]:

    
Tfac = [0.42, 0.54, 0.04]           # Tfactor (square ft/day)
del_h = [0.33, 0.0, -3.5]           # Difference between far-field head and ambient water levels (ft)

To run the forward model mandually, adjust the parameters below and re-execute the parameter cell and the plotting cell below it.

Simulated Flows for Manually Selected Parameters



In [8]:

    
# Calculated Flows for given model parameters
Qamb_sim = calculate_flow(Tfac,Ttot,del_h,Ro,Dw,ambient=True)
print('Ambient Flow:', Qamb_sim)
# Pumped flow
Qstress_sim = calculate_flow(Tfac,Ttot,del_h,Ro,Dw,ambient=False)
print('Pumped Flow:', Qstress_sim)









    



('Ambient Flow:', array([-0.00020001,  0.01980098,  0.01980098]))
('Pumped Flow:', array([ 0.49982481,  0.49982481,  0.22981141]))



In [9]:

    
# Plot the Ambient Flow
# Interpolate flow values
z = arange(0,bot_well,.1)
qa_calc = zeros(len(z))
qs_calc = zeros(len(z))
d = [0, bot_well]
d[1:1] = depth
for i in range(len(depth)):
    qa_calc[np.logical_and(z >= d[i], z < d[i+1])] = Qamb[i]
    qs_calc[np.logical_and(z >= d[i], z < d[i+1])] = Qstress[i]

Automatic Model Calibration

Additional Parameters for Inversion



In [10]:

    
# Minimum allowable transmissivity
Tfac_min = 1.0e-9
# Maximum allowable absolute error in far-field head
del_h_max = 5.0

# Additional parameters that need to be passed but not altered
alpha = 1e-4  # Regularization parameter
beta = 1.0    # Can increase weight to ambient measurements (beta >0)

Initialize parameters values to starting guess



In [11]:

    
# Create array of initial guesses at the fitting parameter values
n = len(Tfac)
init_guess = np.zeros(2*n)

# Starting guess at model paramters based on manual calculation 
init_guess[0:n] = Tfac
init_guess[n:] = del_h

Implement parameter bounds and contraints



In [12]:

    
# Contraints and bounds
# Sum of Tfac = 1, or the difference with 1 is zero
# Define a function generator to create the constraint function because
# the number of fractures is not known in advance
def make_contraint(n):
    def f(x):
        return np.sum(x[0:n]) - 1
    return f
    
# Use the function generator to create the contraint function        
total_one = make_contraint(n)
cons = ({'type': 'eq', 'fun': total_one})

# Create parameter bounds
# Tfac > Tfac_min
# ABS(del_h) < del_h_max
bndT = tuple((Tfac_min,1) for x in Tfac)
bndh = tuple((-del_h_max, del_h_max) for x in del_h)
bnds = bndT + bndh

Define the Objective Function

The objective function involves two terms: the mean-square error for the predicted flow measurements, and the mismatch between the borehole's water level and the far-field heads:

\begin{align} min(F) = \frac{1}{n}\left[\sum_{i=1}^{n} (Q_i^{s,sim} - Q_i^{s,int})^2 + (Q_i^{a,sim} - Q_i^{a,int})^2\right] + \alpha \sum_{i=1}^{n}(\Delta h_i )^2 \end{align}

subject to the constraints \begin{align} T_i^{factor} \geq T_{min}^{factor} \end{align}

\begin{align} abs(h_w^a - h_i^0) \leq \Delta h_{max} \end{align}

for all i.



In [13]:

    
def objective(x):
    '''This is the objective function for the minimization.  I have added a beta parameter to
    allow more or less weight to be placed on the ambient flow data relative to the pumped.'''
    # Unpack vector of parameters (x)
    Tfac = x[0:n]
    del_h = x[n:]    
    # Simulate ambient flow
    Qamb_sim = calculate_flow(Tfac,Ttot,del_h,Ro,Dw,ambient=True)
    # Simulate pumped flow
    Qstress_sim = calculate_flow(Tfac,Ttot,del_h,Ro,Dw,ambient=False)
    F = (np.sum((Qstress_sim - Qstress)**2) + beta*np.sum((Qamb_sim - Qamb)**2))/n
    F += alpha * np.sum(del_h**2)
    return F

Minimize the objective function



In [14]:

    
# Minimize the objective function
# Uses Sequential Quadratric Programs (Kraft, 1998)
res = minimize(objective, init_guess, bounds=bnds, constraints=cons, method='SLSQP')
# Calculate flow from fitted transmissivities and delta heads
Qamb_fit = calculate_flow(res.x[0:n],Ttot,res.x[n:],Ro,Dw,ambient=True)
Qstress_fit = calculate_flow(res.x[0:n],Ttot,res.x[n:],Ro,Dw,ambient=False)

RESULTS

Calculate flow values for a suite of points in the well



In [15]:

    
# Interpolate flow values
z = np.arange(0,bot_well,.1)
qa_calc = np.zeros(len(z))
qs_calc = np.zeros(len(z))
qafit_calc = np.zeros(len(z))
qsfit_calc = np.zeros(len(z))
d = [0, bot_well]
d[1:1] = depth
for i in range(len(depth)):
    qa_calc[np.logical_and(z >= d[i], z < d[i+1])] = Qamb_sim[i]
    qs_calc[np.logical_and(z >= d[i], z < d[i+1])] = Qstress_sim[i]
    qafit_calc[np.logical_and(z >= d[i], z < d[i+1])] = Qamb_fit[i]
    qsfit_calc[np.logical_and(z >= d[i], z < d[i+1])] = Qstress_fit[i]

Plot Data



In [16]:

    
# Plot the results
# Axis limits can be adjusted below with xlim and ylim
import matplotlib.pyplot as plt
fig = plt.figure()
fig.set_size_inches(8,10)
fig.subplots_adjust(wspace=0.5)
ax1 = fig.add_subplot(1,2,1)
ax1.set_xlim(-0.025,0.025)
ax1.set_ylim(bot_casing,bot_well)
#ax1.xaxis.tick_top()
plt.plot(Qa,Za,'ob')
plt.plot(qa_calc,z,'b--', linewidth=3)
plt.plot(qafit_calc,z,'b-')
plt.title('Ambient Flow')
plt.gca().invert_yaxis()
plt.xlabel('Upward Flow (gpm)')
plt.ylabel('Depth in Feet')
ax2 = fig.add_subplot(1,2,2)
plt.plot(Qs,Zs,'or')
plt.plot(qs_calc,z,'r--', linewidth=3)
plt.plot(qsfit_calc,z,'r-')
plt.title('Pumped Flow')
ax2.set_xlim(-0.5,1.0)
ax2.set_ylim(bot_casing,bot_well)
#ax2.xaxis.tick_top()
plt.gca().invert_yaxis()
plt.xlabel('Upward Flow (gpm)')
plt.ylabel('Depth in Feet')
plt.show()

References

Day-Lewis, F. D., Johnson, C. D., Paillet, F. L., & Halford, K. J. (2011). A Computer Program for Flow-Log Analysis of Single Holes (FLASH). Ground Water, 49, 926–931. doi:10.1111/j.1745-6584.2011.00798.x

Kraft D (1988) A software package for sequential quadratic programming. Tech. Rep. DFVLR-FB 88-28, DLR German Aerospace Center — Institute for Flight Mechanics, Koln, Germany.