I have a two layer model. The top layer in unconfined and the bottom layer is confined. I can not seem to get my flopy models vertical conductance to match the results of a PMWIN model where it is computed automatically within the program. Here is the relevant bit of code:
MODEL DIMENSIONS
nlay = 2 # number of layers
nrow = 80 # number of rows
ncol = 57 # number of columns
laytyp = [1,0] #unconfined top layer, confined aquifer below
zbot = -300 # bottom of aquifer same for entire model
Elevations for each model layer (top & bottom) is imported from text-file:
Model_Top_Elevation = np.loadtxt('model_top_elev.txt', dtype=np.float32)
Model_Top_Elevation = np.array([Model_Top_Elevation])
Lay1_Bottom_Elevation= np.loadtxt('bottom_lay1.txt', dtype=np.float32)
# bottom elevation of Layer 1 is imported from txt file
Lay1_Bottom_Elevation= np.loadtxt('bottom_lay1.txt', dtype=np.float32)
Lay1_Bottom_Elevation = np.array([Lay1_Bottom_Elevation])
Lay2_Top_Elevation = np.loadtxt('top_lay2.txt', dtype=np.float32)
Lay2_Top_Elevation = np.array([Lay2_Top_Elevation])
DETERMINE AQUIFER THICKNESS of each model layer:
Lay1_Width = -1*(Lay1_Bottom_Elevation-Model_Top_Elevation)
Lay2_Width = -1*(-1*Lay2_Top_Elevation + -300)
Calculate Lay 1 aquifer thickness if it begins from watertable (not top elevation):
Initial head = 211
Lay1_Watertable_Width = -1*(Lay1_Bottom_Elevation-Initial head)
Hydraulic conductivity of 10m/d adopted for entire model
hk = 10
Vertical conductivity
vk = 0.00001*np.ones((nlay, nrow, ncol), dtype=np.float32)
Calculate transmissivity of layer 2
T = hk*(Lay2_Width)
Calculate vertical conductance (vcont)
vcont = 1/(((Lay1_Watertable_Width)/ vk[0])+((Lay2_Width*0.5)/ vk[1]))
Write BCF File:
bcf = flopy.modflow.ModflowBcf(mf,ipakcb=50,laycon=[1,0], hdry=-1e+30, iwdflg=0, wetfct=0.1, iwetit=1, ihdwet=0, trpy = [1,0],hy=hk, vcont=vcont,tran=[0,T])
MODEL OUTPUT DOES NOT MATCH PMWIN CALCULATION FOR VCONT
If anyone has any recommendations as to where I am going wrong please help!
