У меня есть двухслойная модель. Верхний слой не ограничен, а нижний слой ограничен. Кажется, я не могу добиться того, чтобы вертикальная проводимость моих моделей с флоппи соответствовала результатам модели PMWIN, где она автоматически рассчитывается в программе. Вот соответствующий бит кода:
РАЗМЕРЫ МОДЕЛИ
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
Высота для каждого слоя модели (верх и низ) импортируется из текстового файла:
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)
# нижний уровень слоя 1 импортируется из текстового файла
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])
ОПРЕДЕЛЕНИЕ ТОЛЩИНЫ АКВИФРА каждого модельного слоя:
Lay1_Width = -1*(Lay1_Bottom_Elevation-Model_Top_Elevation)
Lay2_Width = -1*(-1*Lay2_Top_Elevation + -300)
Рассчитайте толщину водоносного горизонта слоя 1, если он начинается с водной таблицы (не высота над уровнем моря):
Initial head = 211
Lay1_Watertable_Width = -1*(Lay1_Bottom_Elevation-Initial head)
Гидравлика c Проводимость 10 м / д, принятая для всей модели
hk = 10
Вертикальная проводимость
vk = 0.00001*np.ones((nlay, nrow, ncol), dtype=np.float32)
Рассчитать коэффициент пропускания слоя 2
T = hk*(Lay2_Width)
Расчет вертикальной проводимости (vcont)
vcont = 1/(((Lay1_Watertable_Width)/ vk[0])+((Lay2_Width*0.5)/ vk[1]))
Запись файла BCF:
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])
ВЫХОД МОДЕЛИ НЕ СООТВЕТСТВУЕТ РАСЧЕТУ PMWIN ДЛЯ VCONT
Если у кого-то есть какие-либо рекомендации относительно того, где я иду не так, пожалуйста, помогите!
