Python Program for Analyzing Leaky Aquifer Pumping Test Data

Introduction to Python Leaky Aquifer Test Analysis Web Page

The purpose of this web page is to describe the application of Python code for pumping test analysis to actual data.  I developed the Python code by translating parts of Fortran code developed by Moench (1985) and adding (1) an input module and (2) output graphics.  The Python program consists of the two modules in Attachment 1 and Attachment 2.  The input module is PYLAQ_INP.py and the main module is PYLAQ.py.  Both modules must be altered when applied to different pumping tests.   Appropriate input values must be entered in PYLAQ_INP.py and the graph title must be entered in PYLAQ.py.  The program provides an efficient way to modify aquifer and aquitard parameters and instantaneously compare the resulting calculated time-drawdown curve with the field data curve.  This web page is technical and assumes knowledge of water well pumping test analysis and Python programming.  The actual data used to illustrate the application of PYLAQ is described by Green and others (1991). 

Description of the Leaky Aquifer Tested

The aquifer that was tested is composed of sand lenses interbedded with silt and clay.  Its thickness is 320 feet.  The aquifer is overlain by 320 feet of clay, siltstone, and very fine to fine grained sandstone, which is in turn overlain by a water table aquifer.  This superjacent 320 foot thick unit is treated as a leaky aquitard in this pumping test analysis.  The aquifer is underlain by 450 feet of sandy clay and siltstone which is in turn underlain by alternating layers of bentonitic clay and siltstone.  This subjacent 450 foot thick unit is treated as a leaky aquitard.

Pumped Well Description

The pumped well was a 20-inch diameter borehole.  The screened interval extended through the entire thickness of the aquifer, and consisted of alternating sections of 8-inch diameter casing and 8-inch diameter wire wrapped screen with 0.025-inch openings.   The total length of screen was 180 feet. The annulus between the alternating screen and casing and the wall of the hole was packed with gravel.

Observation Well Description

One observation well (E-2) was used.  It was cased with 2-inch PVC pipe.  The same intervals were screened as in the pumped well.  The distance from the pumped well to the observation well was 100 feet.    

Description of the Pumping Test and the Original Analysis

The original pumping test is described by Green and others (1991).  The well was pumped at a constant rate of 405 gpm for 98 hours with less than 3 percent discharge variation.  The original analysis used a type curve method developed by Papadopulos and Cooper and presented in Reed (1980). The analysis was applied to time-drawdown data from the the monitoring well.  This method of analysis assumes a fully penetrating well of finite diameter in a non-leaky aquifer.   It improves on the type curve method based on the Theis equation by considering the effect of depletion of water stored in the well resulting from drawdown during the test.

The result was as follows:

Transmissivity: 2,244 gpd/ft

Storativity: 3E-4 (dimensionless)

PYLAQ Constant Rate Pumping Test Analysis

PYLAQ computes curves for constant rate leaky aquifer tests  for three cases:

1. Aquitards above and below the aquifer have constant head boundaries on the top and bottom, respectively.

2. Aquitards above and below the aquifer have no-flow boundaries on the top and bottom, respectively.

3. The aquitard above has a constant head boundary and the aquitard below has a no-flow boundary.

One can also specify whether leakage is from (1) both aquitards, (2) only the superjacent aquitard or (3) only the subjacent aquitard . The pumped well may be line-source or finite-diameter.

PYLAQ calculates curves based on input that includes the following parameters:

• Thickness of the aquifer and aquitards,

• Specific storage of the aquifer and aquitards,

• Horizontal hydraulic conductivity of the aquifer,

• Vertical hydraulic conductivity of the aquitards,

• Distance to  the observation well,

• Effective radius of the pumped well,

• Radius of the pumped well casing where the water level is declining,

• Pumped well “skin” constant related to well loss.

Figure 1 shows the match between the calculated and observed curves.  The calculated curve that matches the observed curve is found by successive trials of values entered into PYLAQ_INP.py.

Figure 1.  Type curve match for observation well E-2.  Green is observed data.  Red is the calculated curve.

The values that produced the calculated curve are given in Attachment 2.  These values include the following:

Aquifer hydraulic conductivity: 0.9375 (ft/day)

Aquifer specific storage:  0.375e-7 (1/ft)

Aquitard vertical hydraulic conductivity:  0.000003 (ft/day)

Aquitard specific storage:  0.0001 (1/ft) 

and the known values given above.

These values for aquifer hydraulic condictivity and specifec storage are the same as those obtained by Green.   The very small aquitard hydraulic conductivity is close to the zero value used by Green.  PYLAQ will not accept zero for aquitard hycraulic conductivity.

Green and others (1991) also provided drawdown data for the pumping well, but they did not analyze it.  Anaysis of this pumping well data using PYLAQ yielded Figure 2.  This result was obtained by increasing the aquifer hydraulic conductivity to 1.01 ft/day, which is close to the value obtained for the observation well 100 feet distant (0.9375 ft/day).  No other aquifer or aquitard parameters were changed.

Figure 2.  Type curve match for the pumping well.  Green is observed data.  Red is the calculated curve.

Attachment 1:  PYLAQ Python Program for Leaky Aquifer Pumping Test Analysis

'''

PYLAQ.py

Author Darrel Dunn

Generates type curves for leaky aquifers.

Well may be line source or large diameter with well bore storage and skin.

Based on DP_LAQ Fortran Program developed by A. F. Moench.  U.S.G.S.

Adjustment to Python zero-based indexing is by not using the zero index.

Input is from PYLAQ_INP.py

This version uses day and foot units.

Any use of this Python code or derivatives is the responsibility of the user.

Under no circumsances shall the author be liable for damages that result from

such use.  This code has not been tested extensively and may contain errors.

'''

import numpy as np

import math

import matplotlib.pyplot as plt

from PYLAQ_INP import  *

import matplotlib

print(matplotlib.__version__)

#Function for stehfest coeffiecients:

def stehfest_coefficients(N):

  """

  Calculates the Stehfest coefficients for a given number of terms N.

  Args:

    N: An even integer representing the number of terms.

  Returns:

    A NumPy array containing the Stehfest coefficients (v_j) for j = 1 to N.

  """

  if N % 2 != 0:

    raise ValueError("N must be an even integer.")

  v = np.zeros(N + 1)

  for j in range(1, N + 1):

    sum_val = 0

    for k in range(int((j + 1) / 2), min(j, N // 2) + 1):

      numerator = k**(N // 2) * math.factorial(2 * k)

      denominator = math.factorial(N // 2 - k) * math.factorial(k) *\

       math.factorial(k - 1) * math.factorial(j - k) * math.factorial(2 * k - j)

      sum_val += numerator / denominator

    v[j] = (-1)**(j + N // 2) * sum_val

  return v[1:] # Return coefficients from v[1] to v[N]

def kzeone(x, y, re0, im0, re1, im1):

    """

    Translates the Fortran subroutine KZEONE to a Python function.

    Args:

        x (float): Double precision real part of the argument.

        y (float): Double precision imaginary part of the argument.

        re0 (float): Output for the real part of K0.

        im0 (float): Output for the imaginary part of K0.

        re1 (float): Output for the real part of K1.

        im1 (float): Output for the imaginary part of K1.

    Returns:

        tuple: A tuple containing (re0, im0, re1, im1).

    """

    tsq = [0.0, 3.19303633920635e-1, 1.29075862295915e0,

           2.95837445869665e0, 5.40903159724444e0, 8.80407957805676e0,

           1.34685357432515e1, 2.02499163658709e1]

    exsq = [0.5641003087264e0, 0.4120286874989e0, 0.1584889157959e0,

            0.3078003387255e-1, 0.2778068842913e-2, 0.1000044412325e-3,

            0.1059115547711e-5, 0.1522475804254e-8]

    r2 = x * x + y * y

    if x > 0.0 or r2 != 0.0:

        if r2 >= 1.96e2:

            rterm = 1.0

            iterm = 0.0

            re0_val = 1.0

            im0_val = 0.0

            re1_val = 1.0

            im1_val = 0.0

            p1 = 8.0 * r2

            p2 = math.sqrt(r2)

            l = 3.91 + 81.2 / p2

            r1 = 1.0

            r2_val = 1.0

            m = -8

            k = 3

            for n in range(1, int(l) + 1):

                m += 8

                k -= m

                r1 *= (k - 4)

                r2_val *= k

                t1 = n * p1

                t2 = rterm

                rterm = (t2 * x + iterm * y) / t1

                iterm = (-t2 * y + iterm * x) / t1

                re0_val += r1 * rterm

                im0_val += r1 * iterm

                re1_val += r2_val * rterm

                im1_val += r2_val * iterm

            t1 = math.sqrt(p2 + x)

            t2 = -y / t1

            p1 = 0.886226925452758 / p2

            rterm = p1 * math.cos(y)

            iterm = -p1 * math.sin(y)

            r1 = re0_val * rterm - im0_val * iterm

            r2_val = re0_val * iterm + im0_val * rterm

            re0 = t1 * r1 - t2 * r2_val

            im0 = t1 * r2_val + t2 * r1

            r1 = re1_val * rterm - im1_val * iterm

            r2_val = re1_val * iterm + im1_val * rterm

            re1 = t1 * r1 - t2 * r2_val

            im1 = t1 * r2_val + t2 * r1

            return re0, im0, re1, im1

        elif r2 >= 18.49:

            x2 = 2.0 * x

            y2 = 2.0 * y

            r1 = y2 * y2

            p1 = math.sqrt(x2 * x2 + r1)

            p2 = math.sqrt(p1 + x2)

            t1 = exsq[1] / (2.0 * p1)

            re0_val = t1 * p2

            im0_val = t1 / p2

            re1_val = 0.0

            im1_val = 0.0

            for n in range(2, 8):

                t2 = x2 + tsq[n]

                p1 = math.sqrt(t2 * t2 + r1)

                p2 = math.sqrt(p1 + t2)

                t1 = exsq[n] / p1

                re0_val += t1 * p2

                im0_val += t1 / p2

                t1 = exsq[n] * tsq[n]

                re1_val += t1 * p2

                im1_val += t1 / p2

            t2 = -y2 * im0_val

            re1 = re1_val / r2

            r2_val = y2 * im1_val / r2

            rterm = 1.41421356237309 * math.cos(y)

            iterm = -1.41421356237309 * math.sin(y)

            im0 = re0_val * iterm + t2 * rterm

            re0 = re0_val * rterm - t2 * iterm

            t1 = re1 * rterm - r2_val * iterm

            t2 = re1 * iterm + r2_val * rterm

            re1 = t1 * x + t2 * y

            im1 = -t1 * y + t2 * x

            return re0, im0, re1, im1

        else:

            x2 = x / 2.0

            y2 = y / 2.0

            p1 = x2 * x2

            p2 = y2 * y2

            t1 = -(math.log(p1 + p2) / 2.0 + 0.5772156649015329)

            t2 = -math.atan2(y, x)

            x2_val = p1 - p2

            y2_val = x * y2

            rterm = 1.0

            iterm = 0.0

            re0_val = t1

            im0_val = t2

            t1_val = t1 + 0.5

            re1_val = t1_val

            im1_val = t2

            p2_val = math.sqrt(r2)

            l = 2.106 * p2_val + 4.4

            if p2_val < 0.8:

                l = 2.129 * p2_val + 4.0

            for n in range(1, int(l) + 1):

                p1_val = float(n)

                p2_val = n * n

                r1 = rterm

                rterm = (r1 * x2_val - iterm * y2_val) / p2_val

                iterm = (r1 * y2_val + iterm * x2_val) / p2_val

                t1_val += 0.5 / p1_val

                re0_val += t1_val * rterm - t2 * iterm

                im0_val += t1_val * iterm + t2 * rterm

                p1_val += 1.0

                t1_val += 0.5 / p1_val

                re1_val += (t1_val * rterm - t2 * iterm) / p1_val

                im1_val += (t1_val * iterm + t2 * rterm) / p1_val

            r1 = x / r2 - 0.5 * (x * re1_val - y * im1_val)

            r2_val = -y / r2 - 0.5 * (x * im1_val + y * re1_val)

            p1_val = math.exp(x)

            re0 = p1_val * re0_val

            im0 = p1_val * im0_val

            re1 = p1_val * r1

            im1 = p1_val * r2_val

            return re0, im0, re1, im1

    else:

        print("ARGUMENT OF THE BESSEL FUNCTIONS IS ZERO OR LIES IN LEFT HALF\

               COMPLEX PLANE")

        return float('nan'), float('nan'), float('nan'), float('nan')

def beskz(z):

    """

    Args:

        z (complex): A complex number.

    Returns:

        tuple: A tuple containing two complex numbers (k0, k1).

    """

    x = z.real

    y = z.imag

    re0, im0, re1, im1 = kzeone(x, y, 0.0, 0.0, 0.0, 0.0)

    k0 = complex(math.exp(-x) * re0, math.exp(-x) * im0)

    k1 = complex(math.exp(-x) * re1, math.exp(-x) * im1)

    return k0, k1

###################################

#Main execution block begins here.

sigp = (Ssover/Ss) * (thicksup/thick) # =(Ss'/Ss)*(b'/b) (1)

print('sigp =',sigp)

gammp = rw * math.sqrt((Ksup/Kaq) * (1/(thicksup*thick))) # (1)

# = rw * sqrt((K'/K)*(1/bb'))

print('gammp =', gammp)

sigpp = (Ssunder/Ss) * (thicksub/thick) # (Ss''/Ss)*(b''/b)  (1)

print('sigpp = ',sigpp)

gammpp = rw * math.sqrt((Ksub/Kaq) * (1/(thicksub*thick))) # (1)

# =rw * sqrt((K''/K)*(1/bb''))

print('gammpp =', gammpp)

rd = r/rw #(1)

print('rD =', rd)

wd = 3.1416*rc**2/(6.2832*rw**2*Ss*thick) # Well bore storage parameter (1)

print('WD =', wd)

v = stehfest_coefficients(n) #Gets v from stehfest_coefficients function. (1)

xln2 = math.log(2.0) # (1)

td = []#To initiate td.

hdplot = [] #To inititate hdplot.

adjtd = []

adjhdplot = []

print("  TD            HD")

x1p = gammp * math.sqrt(sigp) #(1)

x1pp = gammpp * math.sqrt(sigpp) #(1)

argp = math.sqrt(sigp) / gammp #(1)

argpp = math.sqrt(sigpp) / gammpp #(1)

nts = nox * nlc + 1 # Number of dimensionless times for type curve calcs (1)

ded = 1.0 / nox #(1)

deldt = 10.0**ded # /dimensionless time interval for type curve calcs. (1)

time = deldt * rd * rd * tdlast # Initiats time. (1)      

for j in range(1, nts+1): #Number of hd values in a column

    time /= deldt #(1)

    xt = 0.0 #(1)

    xp = 0.0 #(1)

    xpp = 0.0 #(1)

    xppp = 0.0 #(1)

    for i in range(n): 

        pp = xln2 * (i+1) / time #Changed from i to i+1 so treating 

#                                 like Stefest coeficient 1-based. {1)

        eta = math.sqrt(pp) # (1)

        agp = eta * argp # (1)

        agpp = eta * argpp # (1)

        ap = agp * zdp # (1)

        app = agpp * zdpp # (1)

        if ap > 100.0:

                    ap = 100.0

        if app > 100.0:

                    app = 100.0

        if agp > 100.0:

                    agp = 100.0

        if agpp > 100.0:

                    agpp = 100.0

        qp = 0.0

        qpp = 0.0

        if lcase == 1:

                    qp = eta * x1p / math.tanh(agp)# (1)

                    qpp = eta * x1pp / math.tanh(agpp)# (1)

        elif lcase == 2:

                    qp = eta * x1p * math.tanh(agp)# (1)

                    qpp = eta * x1pp * math.tanh(agpp)# (1)

        elif lcase == 3:

                    qp = eta * x1p / math.tanh(agp)# (1)

                    qpp = eta * x1pp / math.tanh(agpp)# (1)

        if iq == 0:

                    qp = 0.0

                    qpp = 0.0

        if iq == 2:

                    qpp = 0.0

        if iq == 3:

                    qp = 0.0

        x0 = eta * eta + qp + qpp # (1)

        x0 = rd * math.sqrt(x0)# (1)

        x0t = rd * eta# (1)

        if ld == 2:

                    z = complex(x0, 0.0) #creates a complex number 

                    re0, re1 = beskz(z) # Multi-variable assignment. (1,1)

                    pdl = re0.real / pp #(1)

        else:

                    xw = x0 / rd # (1)

                    z = complex(xw, 0.0)

                    re0_beskzw, re1_beskzw = beskz(z)

                    xaa = re0_beskzw.real

                    xbb = xw * re1_beskzw.real

                    xcc = wskin * xbb # (1)

                    z = complex(x0, 0.0)

                    re0_beskz0, re1_beskz0 = beskz(z)

                    xdd = re0_beskz0.real

                    xnum = xdd # (1)

                    if rd == 1.0:

                        xnum = xaa + xcc

                    xden = wd * (eta) * (eta) * (xaa + xcc) + xbb #(1)

                    pdl = xnum / (pp * xden) #(1)

        z = complex(x0t, 0.0)

        re0_t, re1_t = beskz(z)

        pdlt = re0_t.real / pp

        pdlp = 0.0

        pdlpp = 0.0

        if iq == 0:

                    pdlp = 0.0

                    pdlpp = 0.0

        if iq == 2:

                    pdlpp = 0.0

        if iq == 3:

                    pdlp = 0.0

        xt += v[i] * pdlt

        xp += v[i] * pdl  

        xpp += v[i] * pdlp 

        xppp += v[i] * pdlpp 

    hd = 2.0 * xp * xln2 / time #(1)

    hdp = 2.0 * xpp * xln2 / time #(1) 

    hdpp = 2.0 * xppp * xln2 / time #(1)

    td.append(time / (rd * rd)) #(1)

    hdplot.append(hd) #(1) 

    print(f"{td[-1]:13.6e} {hdplot[-1]:13.6e}")

adjtd = [item * (Ss/Kaq) for item in td] #(day/(ft*ft))

adjhdplot = [item * (Qwell/(4.*3.1416*thick*Kaq)) for item in hdplot] # (ft)

print()

print('Adjusted td   Adjusted hd')

print('\n'.join(

    (f"{adjtd:13.6e} {adjhdplot:13.6e}"

     for adjtd, adjhdplot in zip(adjtd, adjhdplot))

))

adjt = [item / (1440. * r *r) for item in obst] # (days/(feet*feet))

print()

print('Adjusted t (days/(ft*ft))    Observed drawdown (ft)')

print('\n'.join(f"{adjt:13.6e}{obsdd:13.6e}" for adjt,obsdd in zip(adjt,obsdd)))

#Begin plotting.

plt.figure(figsize=(8.0,6.0)) #Adjust figure size (decimal inches) as needed.

plt.xscale('log')

plt.yscale('log')

plt.minorticks_on()

plt.tick_params(which='minor',length=10,width=2,colors="k",direction='inout')

#Get the current axes object

ax = plt.gca() #Sets ax as name of the axis object of the current axes.

#Set the linewidth for each of the spines (top, bottom, left, right)

#Adjust the linewidth to desired thickness

ax.spines['top'].set_linewidth(4.0)

ax.spines['bottom'].set_linewidth(4.0)

ax.spines['left'].set_linewidth(4.0)

ax.spines['right'].set_linewidth(4.0)

#Set the color for each of the spines.

ax.spines['top'].set_color('k')    # Color for the top spine

ax.spines['bottom'].set_color('k') # Color for the bottom spine

ax.spines['left'].set_color('k')   # Color for the left spine

ax.spines['right'].set_color('k')  # Color for the right spine

#Modify this block to get the labels and title wanted.

ax.set_xlabel('TIME DIVIDED BY RADIUS SQUARED (DAYS/(FT*FT)', fontsize = 16)

ax.set_ylabel('DRAWDOWN (FT)', fontsize = 16)

ax.set_title('LEAKY AQUIFER TEST (OBSERVATION WELL)')

ax.grid()

#Optional: Tick marks thickness.

#plt.tick_params(axis='x', which='both', width=1.5, length=7) 

#plt.tick_params(axis='y', which='both', width=1.5, length=7)

plt.plot(adjtd,adjhdplot,linewidth=8, color ='red') #adjtd is on x-axis.

plt.plot(adjt,obsdd, linewidth=4, color = 'green') # Plots on same grapth.

plt.savefig('adjhd_vs_adjtd.png')

print('Close the show window to complete execution')

plt.show()

print('END of PYLAQ')

Attachment 2:  Input Module for PYLAQ

'''

PYLAQ_INP.py

The program provides input for PYLAQ.py.

It assumes units of feet and days.  Other units would requre a different treatment of units in PYLAQ.py.

''

import numpy as np

n = 8 #Number of terms in the Stehfest algarithm (8 to 12, 8 is usual).

ld = 1 #1 means large diameter well.  2 means line source well.

lcase = 1 

#    1 means upper boundary of the superjacent aquitard and the lower

#    boundary of the subjacent aqitard are constant head boundaries.    

#    2 means the two boundaries indicated above are no-flow boundaries.

#    3 means the upper boundary of the superjacent aquitard is a constant head

#    boundary and the lower boundary of the subjacent aquitard is no-flow.

iq = 1

#    0 means no leakage.  

#    1 means leakage from both aquiards.

#    2 means leakage for overling aquitard only.0

#    3 means leakage from underlying aquitard only.

Ssover = 0.0001 # Specific storage of superjacent aquitard

Ss = 9.375e-7 #specific storage of aquifer.

thicksup =320. # Thickness of superjacent aquitard

thick = 320. # Thickness of aquifer. 

rw = 0.5 # Effective radius of the well.

Ksup =0.000003 # Hydraulic conductivity of superjacent aquitard.

Kaq = 0.9375 # Hydraulic conductivity of aquifer

zdp = 0.0 # Vertical distance above the base of the superjacent aquitard divided

#    by thickness of the superjacent aquitard. (Used when you want to calculate

#    head in the aquitard.)

Ssunder = 0.0001 # Specific storage of subjacent aquitard.

Ksub = 0.000003

thicksub = 450. #Thickness of subjacent aquitard.

zdpp = 0.0     

r= 100. # Distance to opservation point from center of pumped well.

#         Equal to rw for pumped well.

rc = 0.5 # Casing radius.

wskin = 1.0 # Constant of proportionality that related the flux of water through

#             the well skin to head difference across the skin.

tdlast = 1.0e3  # Largest value of dimensionless time that is desired.

nlc = 3 #Total number of logarithmic cyceles on the time scale. 

#        Stehfest algorithm may become unstable is this number is too large.

nox = 3 # Number of equally spaced times per logarithmic cycle.

#xmult = 3.162 #Multiplier for rw. #I do not see this in the program.

Qwell = 77962.5 # Water well pumping rate in field pumping test.

#with open('Pumping_well.txt', 'r') as f:

#    Drawdown = f.read()

#print('Drawdown', Drawdown)

#The following lines loads *.txt and converts columns to individual arrays.

# The text should be space separated columns of time (minutes) and drawdown (ft).

#The following line loads *.txt and converts columns to individual arrays.

column1, column2 = np.loadtxt('Obs_Well_E2.txt', unpack=True) 

obst = column1.tolist() #Converts the numpy array to a list.

obsdd = column2.tolist()

print('TIME, minutes')

print(obst)

print('Drawdown, feet')

print(obsdd)

print('End PYLAQ_INP')

Attachment 3: Drawdown Text File for PYLAQ_INP

2     0.10

3     0.62

4     1.48

5     2.57

6     3.84

7     5.12

8     6.54

9     7.94

10    9.38

12   11.94

15   15.68

20   20.87

30   28.31

40   33.53

50   38.29

60   42.03

70   45.38

80   48.08

90   50.38

100  52.59

120  56.81

150  61.42

200  66.95

300  74.58

400  79.86

500  83.50

600  87.11

700  89.96

800  91.67

900  93.13

1000 94.43

1200 99.12

1500 103.65

2000 108.90

3000 117.16

4000 123.40

5000 127.35

5500 128.35

Bibliography for PYLAQ Leaky Aquifer Pumping Test 

Green, Earl A., Mark T. Anderson, and Darrel D. Sipe (1991): Aquifer Tests and Water-Quality Analyses of the Arikaree Formation Near Pine Ridge, South Dakota; U.S. Geological Survey Water-Resources Investigations Report 91-4005. 🔗

Moench, Allen F. (1985): Transient Flow to a Large-Diameter Well in an Aquifer With Storative Semiconfining Layers; Water Resources Research, Vol. 21, No. 8.

Prickett, T. A. (1965): Type-Curve Solution to Aquifer Tests under Water-Table Conditions; Ground Water, Volume 3, Number 3.

Reed, J. E. (1980): Type Curves for Selected Problems of Flow to Wells in Confined Aquifers; Techniques of Water-Resources Investigations of the United States Geological Survey, Chapter B3. 🔗

Posted 12/3/25