Kansas Geological Survey, Ground Water Series 3, originally published in 1980
Prev Page--Start
C PROGRAM: THEIS
C PURPOSE: CALCULATES THE DRAWDOWN FOR A WELL PUMPING AT
C CONSTANT DISCHARGE FROM AN INFINITE AQUIFER BY THE
C THEIS EQUATION
C
C ARGUMENTS (L IS ARBITRARY LENGTH, T IS ARBITRARY TIME )
C RS CONSTANT OBSERVATION SWITCH; YES OR NO
C TS CONSTANT OBSERVATION TIME SWITCH; YES OR NO
C US UNIT SWITCH, USE GAL, DAYS, FT; YES OR NO
C PRTSW PRINT SWITCH FOR U, W, AND SENSITIVITY; YES OR NO
C SC STORAGE COEFFICIENT FOR AQUIFER (UNITLESS)
C KB TRANSMISSIVITY OF AQUIFER (GAL/DAY/FT) OR (L**2/T)
C Q CONSTANT PUMPAGE OF WELL (GAL/DAY) OR (L**3/T)
C R OBSERVATION DISTANCE FROM WELL (FT) OR (L)
C T OBSERVATION TIME (DAYS) OR (T)
C S DRAWDOWN (FT) OR (L)
C DSDT SENSITIVITY W.R.T. TRANSMISSIVITY (FT**2*DAY/GAL) OR (T/L)
C DSDSC SENSITIVITY W.R.T. STORAGE (FT) OR (L)
C -------------------------------------------------------------
C
C READ IN ARGUMENTS
CHARACTER RS, TS, US, PRTSW
REAL KB
ICOUNT =0
PRINT 1
PRINT 15
15 FORMAT(' DO YOU WANT U, W, AND SENSITIVITY VALUES PRINTED ? ')
READ:PRTSW
1 FORMAT(1H0)
PRINT 14
14 FORMAT(' DO YOU WANT TO USE THE GAL, DAY, FT SYSTEM ? ')
READ:US
PRINT 2
2 FORMAT(' STORAGE COEFFICIENT (UNITLESS)')
READ:SC
PRINT 3
3 FORMAT(' TRANSMISSIVITY (GAL/DAY/FT) OR (L**2/T)')
READ:KB
PRINT 4
4 FORMAT(' CONSTANT PUMPAGE (GAL/DAY) OR (L**3/T)')
READ:Q
PRINT 10
10 FORMAT (' CONSTANT OBSERVATION POINT?
READ:RS
PRINT 11
11 FORMAT (' CONSTANT OBSERVATION TIME? READ:TS
80 IF (ICOUNT .GT. 0 AND. RS EQ. 3HYES) GO TO 70
PRINT 5
5 FORMAT(' OBSERVATION POINT (FT) OR (L)')
READ:R
70 IF (ICOUNT .GT. 0 AND. TS EQ. 3HYES) GO TO 95
PRINT 6
6 FORMAT(' OBSERVATION TIME (DAYS) OR (T)')
READ:T
C CALCULATE THEIS SOLUTION
C THEIS EQUATION FROM 'INTRO TO HYDROLOGY' BY VIESSMAN PAGE 257
C RATIONAL APPROXIMATION FROM 'HANDBOOK OF MATH FUNCTIONS'
C BY ABRAMOWITZ AND STEGUN PAGE 231
95 CONTINUE
U = R*R*SC/(4.0*KB*T)
IF (US .EQ. 3HYES ) U = 7.48*U
IF (V .LE. 0.0) GO TO 96
IF (U .GT. 1.0) GO TO 86
W= -ALOG(U)-.57721566+.99999193*U-.24991055*U*U+
& .05519968*U*U*U-.00976004*U**4+.00107857*U**5
GO TO 87
86 W (EXP(-U)/U)*(U*U+2.334733*U+.250621)/
& (U*U+3.330657*U + 1.681534)
87 S = (Q/(4.0*3.14159*KB))*W
DSDT = (Q/(4.0*3.14159*KB**2))*(-W + EXP(-U))
DSDSC=-(Q/(4.0*3.14159*KB*SC))*EXP(-U)
C OUTPUT RESULTS
90 PRINT 7, S
7 FORMAT(1H0, 'THE DRAWDOWN IS ', F20.8)
IF (PRTSW EQ. 3HYES ) PRINT 8,U
8 FORMAT (' U = ', E20.8)
IF (PRTSW EQ. 3HYES ) PRINT 9,W
9 FORMAT (' W = ', E20.8 )
IF (PRTSW EQ. 3HYES ) PRINT 12,DSDT
12 FORMAT ( ' SENSITIVITY WITH RESPECT TO TRANSMISSIVITY = ', E20.8 )
IF (PRTSW EQ. 3HYES ) PRINT 21,DSDSC
21 FORMAT (' SENSITIVITY WITH RESPECT TO STORAGE = ', E20.8 )
GO TO 97
96 PRINT 13
13 FORMAT (' DRAWDOWN UNDEFINED U TOO SMALL ')
97 PRINT 1
ICOUNT ICOUNT +1
IF (RS EQ. 3HYES AND. TS EQ. 3HYES ) GO TO 99
GO TO 80
99 STOP
END
SUBROUTINE THEIS(SC)KB,Q,R,T,S,DSDT,DSDSC,UNIT)
C PURPOSE CALCULATES THE DRAWDOWN FOR A WELL PUMPING AT
C CONSTANT DISCHARGE FROM AN INFINITE AQUIFER BY THE
C THEIS EQUATION
C
C ARGUMENTS (L IS ARBITRARY LENGTH, T IS ARBITRARY TIME)
C SC STORAGE COEFFICIENT FOR AQUIFER (UNITLESS)
C KB TRANSMISSIVITY OF AQUIFER (GAL/DAY/FT) OF (L**2/T)
C Q CONSTANT PUMPAGE OF WELL (GAL/DAY) OR (L**3/T)
C R OBSERVATION DISTANCE FROM WELL (FT) OR (L)
C T OBSERVATION TIME (DAYS) OR (T)
C S DRAWDOWN (FT) OR (L)
C DSDT SENSITIVITY W.R.T. TRANSMISSIVITY (FT**2*DAY/GAL) OR (T/L)
C DSDSC SENSITIVITY W.R.T. STORAGE (FT) OR (L)
C UNIT SYSTEM OF UNITS; YES FOR GAL, DAY, FT; OTHERWISE
C CONSISTENT UNITS ARE ASSUMED
C -------------------------------------------------------------
C
C CALCULATE THEIS SOLUTION
C THEIS EQUATION FROM 'INTRO TO HYDROLOGY' BY VIESSMAN PAGE 257
C RATIONAL APPROXIMATION FROM 'HANDBOOK OF MATH FUNCTIONS'
C BY ABRAMOWITZ AND STEGUN PAGE 231
REAL KB
CHARACTER UNIT
U = R*R*SC/(4.0*KB*T)
IF (UNIT EQ. 3HYES U = 7.48*U
IF (U .LE. 0.0) GO TO 96
IF (U .GT. 1.0) GO TO 86
W= -ALOG(U)-.57721566+.99999193*U-.24991055*U*U+
& .05519968*U*U*U-.00976004*U**4+.00107857*U**5
GO TO 87
86 W = (EXP(-U)/U)*(U*U+2.334733*U+.250621)/
& (U*U+3.330657*U + 1.681534)
87 S = (Q/(4.0*3.14159*KB))*W
DSDT = (Q/(4.0*3.14159*KB**2))*(-W + EXP(-U))
DSDSC=-(Q/(4.0*3.14159*KB*SC))*EXP(-U)
RETURN
96 PRINT 13
13 FORMAT DRAWDOWN UNDEFINED U TOO SMALL
99 STOP
END
C PROGRAM THEISFIT (TIME SHARING VERSION)
C PURPOSE: TO CALCULATE THE BEST FIT STORAGE AND TRANSMISSIVITY
C BY FITTING THE THEIS EQUATION TO EXPERIMENTAL PUMPTEST DATA IN
C A LEAST SQUARES SENSE.
C
C IMPORTANT VARIABLES: (L IS ARBITRARY LENGTH, T IS ARBITRARY TIME)
C UNIT SYSTEM OF UNITS; YES FOR GAL, DAY, FT; OTHERWISE
C CONSISTENT UNITS ARE ASSUMED.
C GUESS DO YOU WANT TO ENTER AN INITIAL GUESS FOR SC AND KB.?
C YES OR NO IS THE APPROPRIATE RESPONSE.
C SC STORAGE COEFFICIENT OF AQUIFER (UNITLESS)
C KB TRANSMISSIVITY OF AQUIFER (GAL/DAY/FT) OR (L"*2/T)
C Q CONSTANT PUMPAGE OF WELL (GAL/DAY) OR (L*113/T)
C R OBSERVATION DISTANCE FROM WELL (FT) OR (L)
C N NUMBER OF DRAWDOWN-TIME PAIRS TO BE READ
C S(I) EXPERIMENTAL DRAWDOWN AT TIME I (FT) OR (L)
C T(I) THE ITH TIME AT WHICH AN EXPERIMENTAL MEASUREMENT
C FOR THE DRAWDOWN IS MADE (DAYS) OR (T)
C SIGMA THE RMS ERROR IN DRAWDOWN AFTER THE BEST
C FIT HAS BEEN OBTAINED
C ---------------------------------------------------------------------------
DIMENSION S(100),T(100),SP(100),DSDT(100),DSDSC(100)
REAL KB
CHARACTER GUESS, UNIT
C READ IN THE INITIAL DATA
PRINT 801
801 FORMAT (' DO YOU WANT TO USE GAL, DAY, FT SYSTEM ? ')
READ:UNIT
PRINT 105
105 FORMAT (' GUESS FOR STORAGE AND TRANSMISSIVITY ')
READ: GUESS
IF (GUESS NE. 3HYES ) GO TO 114
PRINT 112
112 FORMAT (' ESTIMATE FOR STORAGE ? ')
READ: SC
PRINT 113
113 FORMAT(' ESTIMATE FOR TRANSMISSIVITY ? GAL/DAY/FT OR L**2/T ')
READ: KB
114 PRINT 115
115 FORMAT (' CONSTANT PUMPAGE RATE ? GAL/DAY OR L**3/T ')
READ: Q
PRINT 116
116 FORMAT (' OBSERVATION DISTANCE FROM PUMPING WELL ? FT OR L ')
READ: R
PRINT 117
117 FORMAT(' NUMBER OF DRAWDOWN-TIME PAIRS TO BE READ ? ')
READ: N
C ECHO PRINT THE INITIAL DATA
WRITE (6, 106 ) SC, KB, Q, R, N
106 FORMAT ( ' ECHO THE INITIAL DATA ' / 5H SC =, E15.8, 5H KB =, E15.8
& /4H Q =, E15.8, 4H R =, E15.8, 4H N =, I5 )
C TYPE IN THE DRAWDOWN-TIME PAIRS IN ORDER OF INCREASING TIME
PRINT 118
118 FORMAT (' TYPE IN DRAWDOWN-TIME PAIRS IN ORDER OF INCREASING TIME. ')
DO 120 1 = 1,N
120 READ: S(I),T(I)
C ECHO PRINT THE DRAWDOWN-TIME PAIRS
WRITE (6, 107) (S(I),T(I), I = 1, N)
107 FORMAT (' THE PUMP TEST DATA IN DRAWDOWN-TIME PAIRS '/(2E20.8))
101 FORMAT ( 8F10.0 )
100 FORMAT ( 4F10.0, I10 )
C CALCULATE THE INITIAL GUESS FOR KB AND SC IF NOT GIVEN
NUM = 4
IF (GUESS EQ. 6HYES GO TO 1
SUMLNT = 0.0
SMLNT2 = 0.0
SUMSE = 0.0
SSELNT = 0.0
DO 50 1 = 1,NUM
II = N +1-I
ALOGT = ALOG(T(II))
SUMLNT = SUMLNT + ALOGT
SMLNT2 = SMLNT2 + ALOGT**2
SUMSE = SUMSE + S(II)
50 SSELNT = SSELNT + S(II)*ALOGT
SLOPE = (NUM*IRSSELNT-SUMSE*SUMLNT)/(NUM*SMLNT2-SUMLNT**2)
CON =-(SLOPE*SUMLNT - SUMSE)/NUM
KB = Q/(SLOPE*4.0*3.14159)
ALNTO = -CON/SLOPE
SC = 4.0*KB/( R**C*2*EXP(.5772-ALNTO))
IF (UNIT EQ. 3HYES ) SC = SC/7.48
C PRINT THE CALCULATED GUESS FOR KB AND SC
WRITE (6,110) KB, SC
110 FORMAT (' THE CALCULATED GUESS FOR KB AND SC IS ', 2E16.8 )
1 ICOUNT = 0
5 TEMPKB = KB
TEMPSC = SC
C USE SENSITIVITY ANALYSIS AND LEAST SQUARES TO FIND A BETTER KB AND SC
DO 10 1 = 1, N
10 CALL THEIS(SC,KB,Q,R,T(I),SP(I),DSDT(I),DSDSC(I),UNIT)
SSUS = 0.0
SSUT = 0.0
SUTUS = 0.0
SUSDIF = 0.0
SUTDIF = 0.0
DO 15 1 = 1, N
SSUS = DSDSC(I)11*2+SSUS
SSUT = DSDT(I)**2+SSUT
SUTUS = DSDSC(I)*DSDT(I) + SUTUS
SUSDIF = DSDSC(I)*(S(I)-SP(I)) + SUSDIF
15 SUTDIF = DSDT(I)*(S(I)-SP(I)) + SUTDIF
DELTSC = (SSUT*SUSDIF-SUTUS*SUTDIF)/(SSUS*SSUT-SUTUS**2)
IF (DELTSC .LT. -.95*SC ) DELTSC = -.95*SC
IF (DELTSC .GT. 2.0*SC ) DELTSC = 2.0*SC
SC = SC + DELTSC
IF (SC .GT. 1.0 ) SC = 1.0
DELTKB = (SUTDIF-DELTSC*SUTUS)/SSIJT
IF (DELTKB .LT. -.95*KB ) DELTKB = -.95*KB
IF (DELTKB .GT. 2.0*KB ) DELTKB = 2.0*KB
KB = KB + DELTKB
WRITE (6,104) KB, SC
104 FORMAT( BEST FIT KB AND SC THIS ITERATION ARE ', 2E16.8)
ICOUNT ICOUNT + I
C IF IT DID NOT CONVERGE PRINT MESSAGE
IF ( ICOUNT .GT. 20) GO TO 40
C DO ANOTHER ITERATION IF IT HASNT CONVERGED AND ITERATION LIMIT NOT
C EXCEEDED.
IF (ABS((TEMPKB-KB)/KB) .GT. .001 .OR. ABS((TEMPSC-SC)/SC) .GT. .001)
& GO TO 5
C PRINT BEST FIT DRAWDOWN-TIME PAIRS IF CONVERGENCE OBTAINED
108 FORMAT (' THE BEST FIT DRAWDOWN-TIME PAIRS ')
WRITE (6,108)
DO 20 1 = 1, N
CALL THEIS(SC,KB,Q,R,T(I),SP(I),DSDT(I),DSDSC(I),UNIT)
WRITE (6,109) SP(I), T(I)
109 FORMAT(2E20.8)
20 CONTINUE
C CALCULATE THE RMS ERROR IN DRAWDOWN AFTER BEST FIT
SUM = 0.0
DO 30 1 = 1, N
30 SUM = SUM + (SP(I)-S(I))**2
SIGMA = SQRT(SUM/N)
WRITE(6,111) SIGMA
111 FORMAT(' THE RMS ERROR FOR DRAWDOWN IS ', E20.8)
STOP
40 WRITE (6,125)
125 FORMAT (' DID NOT CONVERGE IN 20 ITERATIONS)
STOP
END
SUBROUTINE THEIS(SC,KB,Q,R,T,S,DSDT,DSDSC$UNIT)
C PURPOSE CALCULATES THE DRAWDOWN FOR A WELL PUMPING AT
C CONSTANT DISCHARGE FROM AN INFINITE AQUIFER BY THE
C THEIS EQUATION
C
C ARGUMENTS (L IS ARBITRARY LENGTH, T IS ARBITRARY TIME)
C SC STORAGE COEFFICIENT FOR AQUIFER (UNITLESS)
C KB TRANSMISSIVITY OF AQUIFER (GAL/DAY/FT) OF (L**2/T)
C Q CONSTANT PUMPAGE OF WELL (GAL/DAY) OR (L**3/T)
C R OBSERVATION DISTANCE FROM WELL (FT) OR (L)
C T OBSERVATION TIME (DAYS) OR (T)
C S DRAWDOWN (FT) OR (L)
C DSDT SENSITIVITY W.R.T. TRANSMISSIVITY (FT**2*DAY/GAL) OR (T/L)
C DSDSC SENSITIVITY W.R.T. STORAGE (FT) OR (L)
C UNIT SYSTEM OF UNITS; YES FOR GAL, DAY, FT; OTHERWISE
C CONSISTENT UNITS ARE ASSUMED
C-------------------------------------------------------------
C
C CALCULATE THEIS SOLUTION
C THEIS EQUATION FROM 'INTRO TO HYDROLOGY' BY VIESSMAN PAGE 257
C RATIONAL APPROXIMATION FROM 'HANDBOOK OF MATH FUNCTIONS'
C BY ABRAMOWITZ AND STEGUN PAGE 231
REAL KB
CHARACTER UNIT
U = R*R*SC/(4.O*KB*T)
IF (UNIT EQ. 3HYES) U = 7.48*U
IF (U .LE. 0.0) GO TO 96
IF (U .GT. 1.0) GO TO 86
W= -ALOG(U)-.57721566+.99999193*U-.24991055*U*U+
& .05519968*U*U*U-.00976004*U**4+.00107857*U**5
GO TO 87
86 W = (EXP(-U)/U)*(U*U+2.334733*U+.250621)/
& (U*U+3.330657*U + 1.681534)
87 S = (Q/(4.0*3.14159*KB))*W
DSDT = (Q/(4.0*3.14159*KB**2))*(-W + EXP(-U))
DSDSC=-(Q/(4.0*3.14159"KB*SC))*EXP(-U)
RETURN
96 PRINT 13
13 FORMAT DRAWDOWN UNDEFINED U TOO SMALL
99 STOP
END
C PROGRAM THEISFIT (FORTRAN IV BATCH VERSION)
C PURPOSE: TO CALCULATE THE BEST FIT STORAGE AND TRANSMISSIVITY
C BY FITTING THE THEIS EQUATION TO EXPERIMENTAL PUMPTEST DATA
C IN A LEAST SQUARES SENSE.
C
C IMPORTANT VARIABLES: (L IS ARBITRARY LENGTH, T IS ARBITRARY TIME)
C UNIT SYSTEM OF UNITS; YES FOR GAL, DAY, FT; OTHERWISE
C CONSISTENT UNITS ARE ASSUMED.
C GUESS DO YOU WANT TO ENTER AN INITIAL GUESS FOR SC AND KB?
C YES OR NO IS THE APPROPRIATE RESPONSE.
C SC STORAGE COEFFICIENT OF AQUIFER (UNITLESS)
C KB TRANSMISSIVITY OF AQUIFER (GAL/DAY/FT) OR (L**2/T)
C Q CONSTANT PUMPAGE OF WELL (GAL/DAY) OR (L**3/T)
C R OBSERVATION DISTANCE FROM WELL (FT) OR (L)
C N NUMBER OF DRAWDOWN-TIME PAIRS TO BE READ
C S(I) EXPERIMENTAL DRAWDOWN AT TIME I (FT) OR (L)
C T(I) THE ITH TIME AT WHICH AN EXPERIMENTAL MEASUREMENT
C FOR THE DRAWDOWN IS MADE (DAYS) OR (T)
C SIGMA THE RMS ERROR IN DRAWDOWN AFTER THE BEST
C FIT HAS BEEN OBTAINED
C ---------------------------------------------------------------------------
DIMENSION S(100), T(100),SP(100),DSDT(100),DSDSC(100)
REAL KB
CHARACTER GUESS, UNIT
C READ IN INITIAL DATA
READ (5,105) GUESS, UNIT
105 FORMAT ( 2A6 )
READ (5,100) SC, KB, Q, R, N
C ECHO PRINT THE INITIAL DATA
WRITE (6, 106 ) SC, KB, Q, R, N
106 FORMAT ( ' ECHO THE INITIAL DATA ' / 5H SC =, E15.8, 5H KB =, E15.8
+ /4H Q =, E15.8, 4H R =, E15.8, 4H N =, I5 )
C READ IN DRAWDOWN-TIME PAIRS IN ORDER OF INCREASING TIME
READ (5, 101) (S(I),T(I), I = 1, N)
C ECHO PRINT THE DRAWDOWN-TIME PAIRS
WRITE (6, 107) (S(I),T(I), I = 1, N)
107 FORMAT THE PUMP TEST DATA IN DRAWDOWN-TIME PAIRS '/(2E20.8))
101 FORMAT 8F10.0 )
100 FORMAT 4F10.0, I10 )
C CALCULATE THE INITIAL GUESS FOR KB AND SC IF NOT GIVEN
NUM = 4
IF (GUESS EQ. 6HYES ) GO TO 1
SUMLNT = 0.0
SMLNT2 = 0.0
SUMSE = 0.0
SSELNT = 0.0
DO 50 1 = 1,NUM
II = N + 1 -1
ALOGT = ALOG(T(II))
SUMLNT = SUMLNT + ALOGT
SMLNT2 = SMLNT2 + ALOGT**2
SUMSE = SUMSE + S(II)
50 SSELNT = SSELNT + S(II)*ALOGT
SLOPE = (NUM*SSELNT-SUMSE*SUMLNT)/(NUM*SMLNT2-SUMLNT**2)
CON = (SLOPE*SUMLNT - SUMSE)/NUM
KB = Q/(SLOPE*4.0*3.14159)
ALNTO = -CON/SLOPE
SC = 4.0*KB/( R**2*EXP(.5772-ALNTO))
IF (UNIT EQ. 3HYES ) SC = SC/7.48
C PRINT THE CALCULATED GUESS FOR KB AND SC
WRITE (6,110) KB, SC
110 FORMAT (' THE CALCULATED GUESS FOR KB AND SC IS ', 2E16.8)
1 ICOUNT = 0
5 TEMPKB = KB
TEMPSC = SC
C USE SENSITIVITY ANALYSIS AND LEAST SQUATES TO FIND A BETTER KB AND SC DO 10 1 = 1, N
10 CALL THEIS(SC)KB,Q,R,T(I),SP(I),DSDT(I),DSDSC(I),UNIT)
SSUS = 0.0
SSUT = 0.0
SUTUS = 0.0
SUSDIF = 0.0
SUTDIF = 0.0
DO 15 1 = 1, N
SSUS = DSDSC(I)**2+SSUS
SSUT = DSDT(I)**2+SSUT
SUTUS = DSDSC(I)*DSDT(I) + SUTUS
SUSDIF = DSDSC(I)*(S(I)-SP(I)) + SUSDIF
15 SUTDIF = DSDT(I)*(S(I)-SP(I)) + SUTDIF
DELTSC = (SSUT*SUSDIF-SUTUS*SUTDIF)/(SSUS*SSUT-SUTUS**2)
IF (DELTSC .LT. -.95*SC ) DELTSC = -.95*SC
IF (DELTSC .GT. 2.0*SC ) DELTSC = 2.0*SC
SC = SC + DELTSC
IF (SC .GT. 1.0 ) SC = 1.0
DELTKB = (SUTDIF-DELTSC*SUTUS)/SSUT
IF (DELTKB .LT. -.95*KB) DELTKB = -.95*KB
IF (DELTKB .GT. 2.0*KB DELTKB = 2.0*KB
KB = KB + DELTKB
WRITE (6,104) KB, SC
104 FORMAT( ' BEST FIT KB AND SC THIS ITERATION ARE ', 2E16.8)
ICOUNT ICOUNT + 1
C IF IT DID NOT CONVERGE PRINT MESSAGE
IF ( ICOUNT .GT. 20) GO TO 40
C DO ANOTHER ITERATION IF IT HASNT CONVERGED AND ITERATION LIMIT NOT
C EXCEEDED.
IF(ABS((TEMPKB-KB)/KB) .GT. .001 .OR. ABS((TEMPSC-SC)/SC) .GT. .001)
& GO TO 5
C PRINT BEST FIT DRAWDOWN-TIME PAIRS IF CONVERGENCE OBTAINED
108 FORMAT (' THE BEST FIT DRAWDOWN-TIME PAIRS ')
WRITE (6,108)
DO 20 1 = 1, N
CALL THEIS(SC,KB,Q,R,T(I),SP(I),DSDT(I),DSDSC(I),UNIT)
WRITE (6,109) SP(I), T(I)
109 FORMAT(2E20.8)
20 CONTINUE
C CALCULATE THE RMS ERROR IN DRAWDOWN AFTER BEST FIT
SUM = 0. 0
DO 30 1 = 1, N
30 SUM = SUM + (SP(I)-S(I))**2
SIGMA = SQRT(SUM/N)
WRITE(6,111) SIGMA
111 FORMAT(' THE RMS ERROR FOR DRAWDOWN IS ', E20.8)
STOP
40 WRITE (6,112)
112 FORMAT(' DID NOT CONVERGE IN 20 ITERATIONS ')
STOP
END
SUBROUTINE THEIS(SC,KB,Q,R,T,S,DSDT,DSDSC,UNIT)
C PURPOSE CALCULATES THE DRAWDOWN FOR A WELL PUMPING AT
C CONSTANT DISCHARGE FROM AN -INFINITE AQUIFER BY THE
C THEIS EQUATION
C
C ARGUMENTS (L IS ARBITRARY LENGTH, T IS ARBITRARY TIME)
C SC STORAGE COEFFICIENT FOR AQUIFER (UNITLESS)
C KB TRANSMISSIVITY OF AQUIFER (GAL/DAY/FT) OR (L**2/T)
C Q CONSTANT PUMPAGE OF WELL (GAL/DAY) OR (L**3/T)
C R OBSERVATION DISTANCE FROM WELL (FT) OR (L)
C T OBSERVATION TIME (DAYS) OR (T)
C S DRAWDOWN (FT) OR (L)
C DSDT SENSITIVITY W.R.T. TRANSMISSIVITY (FT**2*DAY/GAL) OR (T/L)
C DSDSC SENSITIVITY W.R.T. STORAGE (FT) OR (L)
C UNIT SYSTEM OF UNITS; YES FOR GAL, DAY, FT; OTHERWISE
C CONSISTENT UNITS ARE ASSUMED
C -------------------------------------------------------------
C
C CALCULATE THEIS SOLUTION
C THEIS EQUATION FROM 'INTRO TO HYDROLOGY' BY VIESSMAN PAGE 257
C RATIONAL APPROXIMATION FROM 'HANDBOOK OF MATH FUNCTIONS'
C BY ABRAMOWITZ AND STEGUN PAGE 231
REAL KB
CHARACTER UNIT
U = R*R*SC/(4.0*KB*T)
IF (UNIT EQ. 6HYES ) U = 7.48*U
IF (U .LE. 0.0) GO TO 96
IF (U .GT. 1.0) GO TO 86
W= -ALOG(U)-.57721566+.99999193*U-.24991055*U*U+
& 1.05519968*U*U*U-.00976004*U**4+.00107857*U**5
GO TO 87
86 W = (EXP(-U)/U)*(U*U+2.334733*U+.250621)/
& (U*U+3.330657*U + 1.681534)
87 S = (Q/(4.0*3.14159*KB))*W
DSDT = (Q/(4.0*3.14159*KB**2))*(-W + EXP(-U))
DSDSC = -(Q/(4.0*3.14159*KB*SC))*EXP(-U)
RETURN
96 PRINT 13
13 FORMAT ( ' DRAWDOWN UNDEFINED U TOO SMALL ')
99 STOP
END
Kansas Geological Survey, The Theis Equation
Placed on web June 18, 2012; originally published in 1980.
Comments to webadmin@kgs.ku.edu
The URL for this page is http://www.kgs.ku.edu/Publications/Bulletins/GW3/02_append.html