Document ID: EPA-HQ-OW-2009-0819-9928
Agency: epa
Document Type: Supporting & Related Material
Title: 
Posted Date: 2023-03-29T04:00Z

# -*- coding: utf-8 -*-
"""
Code for Landfills

This Python script calculates pumping volume estimates for potential landfill leachate pumping systems from data facility specific data and input parameter from the SE10413_Landfill Leachate Analysis_Calculations.xlsx.

SE10413_Landfill Leachate Analysis_Calculations.xlsx must be saved in the same folder as this script for the code to successfully run.

Created on Mon Sep 19 21:17:57 2022

@author: Andrew Rimelman, PG Environmental
"""

import pandas as pd
import numpy as np
from scipy.special import exp1

#Calculates aquifer transmissibility
def transmissivity (k,b):
    trans=k*b
    return trans

#Calculates of the dimensionless parameter u
def ufunct(r,s,trans,time):
    u=1.87*r**2*s/(trans*time)
    return u

#Calculates pumping required to reach a drawdown of at a radius r and time t
def qpump(w,T,s):
    Q=s*T/(114.6*w)
    return Q

#Calculates drawndown using Theis Equation
def theis (q, T, r0, S, t):
    u0=ufunct(r0, S, T, t)
    w0=exp1(u0)
    dd=114.6*q/T*w0
    return dd

#Impots and cleans up data from the Leachate Facilites Spreadsheet (Note: Leachate Facilities Spreadsheet must be located in the same folder as this script.)
indata=pd.read_excel("SE10413_Landfill Leachate Analysis_Calculations.xlsx", "Landfill Data", header=2)
pddata=indata.iloc[:,1:]
data=pddata.to_numpy()
inputs=pd.read_excel("SE10413_Landfill Leachate Analysis_Calculations.xlsx", "Theis Equation Inputs", header=1)
inputs=inputs.iloc[:,1:]
wellnotab=pd.read_excel("SE10413_Landfill Leachate Analysis_Calculations.xlsx", "Well Number Buckets", header=1)
wellnotab=wellnotab.iloc[:,1:]
sumstat=pd.read_excel("SE10413_Landfill Leachate Analysis_Calculations.xlsx", "Summary Statistics", header=1)
sumstat=sumstat.iloc[:,1:]

#Assigns variables to Equation inputs
s=float(inputs.iloc[0,1])
mos=float(inputs.iloc[1,1])
t=float(inputs.iloc[2,1])
dd=float(inputs.iloc[3,1])
r0=float(inputs.iloc[4,1])
ddratio=float(inputs.iloc[5,1])
maxit=float(inputs.iloc[6,1])
minpump=float(inputs.iloc[7,1])
jc=inputs.iloc[14,1]

#For loop variables
count=0
exiter=0
speccount=0
exwellcount=0
quaderrorcount=0
#wellno=5.

#Loops through all landfills and calculates pumping requirements at all unlined facilities that are marked "Y" in the test column.
for i in data[:,0]:
    if data[count,32]=="Y":
        #Facility Specific parameters
        k=float(data[count,34])
        b=float(data[count,22])
        xc=float(data[count,25])
        
        #Applies Jacobs correction to drawndown if enabled
        if jc == "Y":
            ddc1=(-1.+(1.**2.-4./(2.*b)*dd)**(1./2.))/(-1./b)
            ddc2=(-1.-(1.**2.-4./(2.*b)*dd)**(1./2.))/(-1./b)
            if abs(ddc1-dd) < abs(ddc2-dd):
                if ddc1 >= 0.:
                    ddu=ddc1
                else:
                    ddu=dd
                    quaderrorcount += 1
                    
            if abs(ddc1-dd) >= abs(ddc2-dd):
                if ddc2 >= 0.:
                    ddu=ddc2
                else:
                    ddu=dd
                    quaderrorcount += 1
            print (ddu)
        #While loop variables
        count3=0
        test1=True
        if jc !="Y":
            ddu=dd
            
        #Verifies whether drawdown at the well is sufficiently less than the aquifer thickness and adds wells accordingly.
        while test1 == True:
            #Sets number of wells for first iteration
            if count3==0:
                counttwo=0
                for j in wellnotab.iloc[:,0]:
                    if xc > j and xc <= wellnotab.iloc[counttwo,1]:
                        wellno= wellnotab.iloc[counttwo,2]
                        break
                    counttwo += 1
            #Calculates pumping rate for individual wells
            r=xc/wellno/2
            tr=transmissivity(k,b)
            u=ufunct(r,s,tr,t)
            w=exp1(u)
            qp=qpump(w,tr,ddu)
            
            #Tests for adequate aquifer thickness
            dd0=theis(qp, tr, r0, s, t)
            maxdd=ddratio*b
            if dd0<=maxdd:
                test1=False
            
            #Tests for maximum iteration limit
            elif count3>maxit:
                test1=False
                exiter+=1
                #print (dd0)
            
            #Adds five wells and repeats if drawdown exceeds aquifer thickness
            else:
                wellno += 5
                print (dd0)
            count3 += 1
        print (exiter)
        
        #Calculates total flow rate for well system
        qtotal=wellno*qp*(1.+mos/100.)
        
        #Counts facilies that require more wells than the initial tabulated values
        if count3 > 1:
            exwellcount += 1
            
        #Inputs calculated values into dataframe
        pddata.iloc[count,35]=tr
        pddata.iloc[count,36]=wellno
        pddata.iloc[count,37]=r
        pddata.iloc[count,38]=u
        pddata.iloc[count,39]=w
        pddata.iloc[count,40]=qp
        pddata.iloc[count,41]=qtotal
        pddata.iloc[count,42]=dd0
        pddata.iloc[count,43]=count3
    count+=1    
    #print(pddata.iloc[count,40])

#Produces Summary Statistics
pumpdata=pddata.iloc[:,40:42]
statarr=pddata["Total Groundwater Pumping Rate (gpm)"]
vallist=[]
vallistadj=[]

countstat=0

#Removes "nan" values
for i in statarr:
    tst=np.isnan(i)
    if tst==False:
        vallist.append(i)
        if pumpdata.iloc[countstat,0] >= minpump:
            vallistadj.append(i)
    countstat += 1
    
#Summary stats (Smaller flows included)        
sumstat.iloc[0,1]=len(vallist)
sumstat.iloc[1,1]=0
sumstat.iloc[2,1]=np.mean(vallist)
sumstat.iloc[3,1]=np.std(vallist)
sumstat.iloc[4,1]=np.min(vallist)
sumstat.iloc[5,1]=np.percentile(vallist,10)
sumstat.iloc[6,1]=np.percentile(vallist,25)
sumstat.iloc[7,1]=np.percentile(vallist,50)
sumstat.iloc[8,1]=np.percentile(vallist,75)
sumstat.iloc[9,1]=np.percentile(vallist,90)
sumstat.iloc[10,1]=np.max(vallist)
sumstat.iloc[11,1]=exwellcount

#Summary stats (Minimum flow rate per pump appied)
sumstat.iloc[0,2]=len(vallistadj)
sumstat.iloc[1,2]=minpump
sumstat.iloc[2,2]=np.mean(vallistadj)
sumstat.iloc[3,2]=np.std(vallistadj)
sumstat.iloc[4,2]=np.min(vallistadj)
sumstat.iloc[5,2]=np.percentile(vallistadj,10)
sumstat.iloc[6,2]=np.percentile(vallistadj,25)
sumstat.iloc[7,2]=np.percentile(vallistadj,50)
sumstat.iloc[8,2]=np.percentile(vallistadj,75)
sumstat.iloc[9,2]=np.percentile(vallistadj,90)
sumstat.iloc[10,2]=np.max(vallistadj)
sumstat.iloc[11,2]=exwellcount

print(vallist)
#Exports Dataframe to spreadsheet
with pd.ExcelWriter('Leachate Output.xlsx') as writer:  
    pddata.to_excel(writer, sheet_name='Output Data', header=True, index=False)
    sumstat.to_excel(writer, sheet_name='Summary Statistics', header=True, index=False)

print (pddata)
print("Data Exported")

#End of landfill leachate script

# -*- coding: utf-8 -*-
"""
Code for Impoundments

This Python script calculates pumping volume estimates for potential impoundment leachate pumping systems from data facility specific data and input parameter from the SE10412_Impoundment Leachate Analysis_Calculations.xlsx spreadsheet.

The SE10412_Impoundment Leachate Analysis_Calculations.xlsx spreadsheat must be saved in the same folder as this script for the code to successfully run.

Created on Mon Sep 19 21:17:57 2022

@author: Andrew Rimelman, PG Environmental
"""

import pandas as pd
import numpy as np
from scipy.special import exp1

#Calculates aquifer transmissibility
def transmissivity (k,b):
    trans=k*b
    return trans

#Calculates of the dimensionless parameter u
def ufunct(r,s,trans,time):
    u=1.87*r**2*s/(trans*time)
    return u

#Calculates pumping required to reach a drawdown of at a radius r and time t
def qpump(w,T,s):
    Q=s*T/(114.6*w)
    return Q

#Calculates drawndown using Theis Equation
def theis (q, T, r0, S, t):
    u0=ufunct(r0, S, T, t)
    w0=exp1(u0)
    dd=114.6*q/T*w0
    return dd

#Impots and cleans up data from the Leachate Facilities Spreadsheet. Combines individual pond and pond groupins data
indata=pd.read_excel(SE10412_Impoundment Leachate Analysis_Calculations.xlsx, "Facilities", header=2)
pdpdata=indata.iloc[:508,1:]
infacdata=pd.read_excel(SE10412_Impoundment Leachate Analysis_Calculations.xlsx, "Pond Group Entries", header=2)
pdfacdata=infacdata.iloc[:,1:]
combdata=[pdpdata, pdfacdata]
pddata=pd.concat(combdata, ignore_index=True)
data=pddata.to_numpy()
inputs=pd.read_excel(SE10412_Impoundment Leachate Analysis_Calculations.xlsx, "Theis Equation Inputs", header=1)
inputs=inputs.iloc[:,1:]
wellnotab=pd.read_excel(SE10412_Impoundment Leachate Analysis_Calculations.xlsx, "Well Number Buckets", header=1)
wellnotab=wellnotab.iloc[:,1:]
sumstat=pd.read_excel(SE10412_Impoundment Leachate Analysis_Calculations.xlsx, "Summary Statistics", header=1)
sumstat=sumstat.iloc[:,1:]

#Assigns variables to Equation inputs
s=float(inputs.iloc[0,1])
mos=float(inputs.iloc[1,1])
t=float(inputs.iloc[2,1])
dd=float(inputs.iloc[3,1])
r0=float(inputs.iloc[4,1])
ddratio=float(inputs.iloc[5,1])
maxit=float(inputs.iloc[6,1])
minpump=float(inputs.iloc[7,1])
jc=inputs.iloc[14,1]

#For loop variables
count=0
exiter=0
speccount=0
exwellcount=0
quaderrorcount=0
#wellno=5.

#Loops through all landfills and calculates pumping requirements at all unlined facilities that are marked "Y" in the test column.
for i in data[:,0]:
    if data[count,31]=="Y":
        #Facility Specific parameters
        k=float(data[count,33])
        b=float(data[count,19])
        xc=float(data[count,22])
        
        #Applies Jacobs correction to drawndown if enabled
        if jc == "Y":
            ddc1=(-1.+(1.**2.-4./(2.*b)*dd)**(1./2.))/(-1./b)
            ddc2=(-1.-(1.**2.-4./(2.*b)*dd)**(1./2.))/(-1./b)
            if abs(ddc1-dd) < abs(ddc2-dd):
                if ddc1 >= 0.:
                    ddu=ddc1
                else:
                    ddu=dd
                    quaderrorcount += 1
                    
            if abs(ddc1-dd) >= abs(ddc2-dd):
                if ddc2 >= 0.:
                    ddu=ddc2
                else:
                    ddu=dd
                    quaderrorcount += 1
            print (ddu)
        #While loop variables
        count3=0
        test1=True
        if jc !="Y":
            ddu=dd
            quaderrorcount += 1
        #Verifies whether drawdown at the well is sufficiently less than the aquifer thickness and add wells accordingly.
        while test1 == True:
            #Sets number of wells for first iteration
            if count3==0:
                counttwo=0
                for j in wellnotab.iloc[:,0]:
                    if xc > j and xc <= wellnotab.iloc[counttwo,1]:
                        wellno= wellnotab.iloc[counttwo,2]
                        break
                    counttwo += 1
            #Calculates pumping rate for individual wells
            r=xc/wellno/2
            tr=transmissivity(k,b)
            u=ufunct(r,s,tr,t)
            w=exp1(u)
            qp=qpump(w,tr,ddu)
            
            #Tests for adequate aquifer thickness
            dd0=theis(qp, tr, r0, s, t)
            maxdd=ddratio*b
            if dd0<=maxdd:
                test1=False
            
            #Tests for maximum iteration limit
            elif count3>maxit:
                test1=False
                exiter+=1
                print (dd0)
            
            #Adds five wells and repeats if drawdown exceeds thickness
            else:
                wellno += 5
                print (dd0)
            count3 += 1
        print (exiter)
        
        #Calculates total flow rate for well system
        qtotal=wellno*qp*(1.+mos/100.)
        
        #Counts facilies that require more wells than the initial tabulated values
        if count3 > 1:
            exwellcount += 1
            
        #Inputs calculated values into dataframe
        pddata.iloc[count,34]=tr
        pddata.iloc[count,35]=wellno
        pddata.iloc[count,36]=r
        pddata.iloc[count,37]=u
        pddata.iloc[count,38]=w
        pddata.iloc[count,39]=qp
        pddata.iloc[count,40]=qtotal
        pddata.iloc[count,41]=dd0
        pddata.iloc[count,42]=count3
    count+=1    
    #print(pddata.iloc[count,40])

#Produces Summary Statistics
pumpdata=pddata.iloc[:,39:41]
statarr=pddata["Total Groundwater Pumping Rate (gpm)"]
vallist=[]
vallistadj=[]

countstat=0

#Removes "nan" values
for i in statarr:
    tst=np.isnan(i)
    if tst==False:
        vallist.append(i)
        if pumpdata.iloc[countstat,0] >= minpump:
            vallistadj.append(i)
    countstat += 1
    
#Summary stats (Smaller flows included)        
sumstat.iloc[0,1]=len(vallist)
sumstat.iloc[1,1]=0
sumstat.iloc[2,1]=np.mean(vallist)
sumstat.iloc[3,1]=np.std(vallist)
sumstat.iloc[4,1]=np.min(vallist)
sumstat.iloc[5,1]=np.percentile(vallist,10)
sumstat.iloc[6,1]=np.percentile(vallist,25)
sumstat.iloc[7,1]=np.percentile(vallist,50)
sumstat.iloc[8,1]=np.percentile(vallist,75)
sumstat.iloc[9,1]=np.percentile(vallist,90)
sumstat.iloc[10,1]=np.max(vallist)
sumstat.iloc[11,1]=exwellcount

#Summary stats (Minimum flow rate per pump appied)
sumstat.iloc[0,2]=len(vallistadj)
sumstat.iloc[1,2]=minpump
sumstat.iloc[2,2]=np.mean(vallistadj)
sumstat.iloc[3,2]=np.std(vallistadj)
sumstat.iloc[4,2]=np.min(vallistadj)
sumstat.iloc[5,2]=np.percentile(vallistadj,10)
sumstat.iloc[6,2]=np.percentile(vallistadj,25)
sumstat.iloc[7,2]=np.percentile(vallistadj,50)
sumstat.iloc[8,2]=np.percentile(vallistadj,75)
sumstat.iloc[9,2]=np.percentile(vallistadj,90)
sumstat.iloc[10,2]=np.max(vallistadj)
sumstat.iloc[11,2]=exwellcount

print(vallist)

#Produces Summary Dataframe

sumdatahead=["ID", "State","EIA Plant Code", "Plant Name", "Unit Name","Clean Closed or Composite Lined?","Pond Group Number", "Pumping Rate Calculated","Number of Wells","Groundwater Pumping Rate (gpm) per Well","Total Groundwater Pumping Rate (gpm)"]
sumdata=pddata[sumdatahead].to_numpy()
sumdataout=np.array([sumdatahead])    
countsum=0
for i in sumdata[:,6]:
   if sumdata[countsum, 6] == "N":
       newrow=sumdata[countsum,:]
       sumdataout=np.vstack((sumdataout,newrow))
   countsum += 1
   
countsum2=0
for j in sumdata[:,0]:
    if type(j)==str:     
        newrow=sumdata[countsum2,:]
        sumdataout=np.vstack((sumdataout,newrow))
    countsum2 += 1
   
sumout = pd.DataFrame(sumdataout, columns=sumdatahead)

#Exports Dataframe to spreadsheet
with pd.ExcelWriter('Leachate Output_Impoundments.xlsx') as writer:  
    pddata.to_excel(writer, sheet_name='Output Data', header=True, index=False)
    sumstat.to_excel(writer, sheet_name='Summary Statistics', header=True, index=False)
    sumout.to_excel(writer, sheet_name='Summary Sheet', header=True, index=False)

print (pddata)
print("Data Exported")

#End of impoundment leachate script.