實現K值隨深度衰減

做水文地質數值模擬的時候，滲透係數是必不可少的引數，當模型模擬範圍達到區域尺度，同時考慮地下水深迴圈時，深部的岩層在上覆地層的作用力下被壓實，地層滲透率會降低，滲透率隨深度的增加呈指數式下降，這是滲透率的深度衰減作用。modflow沒有描述相應變化的程式包，因此我們將滲透係數從h5檔案中提取出來在外部處理，再重新輸入回模型。

def
k_decay
(ks,z)
: 
z=z*
0.001
# z單位是m，將其轉化為km
rho=
0.997074
# 密度 g/cm3
mu=0.0000011413
# 運動粘度(15℃) 相當於 動力粘度 1.1404e-3 pa·s
g=9.8# 重力加速度 m/s2
beta=
0.25
# 滲透率衰減係數
kr=10**(-
25.4
)# 衰減殘餘滲透率 log(kr)=-25.4
convt=1/
24.0
/3600.0
*mu/g/rho # m/d -> m/s; k=k*rho*g/mu
ks=ks*convt
k=kr*
(ks/kr)**(
(1+z)**
(-beta)
) k=k/convt
return k

上述函式實現深度為z處對應k的計算，傳入地表滲透係數k和地層深度z，返回該處的滲透係數。接下來要從h5檔案中獲取地表的 k 和 z 列表資料，傳入上述函式計算。

採用下面**獲取k和z資料，需要專案名稱，模型層數。

def
read_k_ele
(filename_src,nlay)
:# modflow專案名稱，模型層數
with h5py.file(filename_src,
'r')
as f:
top=np.array(f[
'arrays/top1'])
# 地表高程
hk1=np.array(f[
'arrays/hk1'])
# 地表滲透係數，在這裡我用的是地層**，如果hk1是滲透係數，可忽略下面6行**
lyrt=np.array(f[
'arrays/bot1'])
# 當層的頂板高程
df=pd.dataframe(hk1,columns=
['hk1'])
k_code=pd.read_table(
"k_code.csv"
,sep=
',')
# 獲取地層**對應的滲透係數
df=pd.merge(df,k_code,on=
'hk1'
)# 在資料表中加入滲透係數，相當於arcgis的join功能
hk1=df.k.values
hk=[[
0]*len
(hk1)
]*nlay # 初始化滲透係數儲存陣列
hk[0]
=hk1
for k in
range(2
,nlay+1)
: 
lyrb=np.array(f[
'arrays/bot%i'
%k])
# 當層的地板高程
lyrc=
(lyrb+lyrt)
*0.5
# 含水層中部高程
lyrt=lyrb
z=top-lyrc
hk[k-1]
=k_decay(hk1,z)
return hk

最後，寫入h5副本

def
main()
: 
filename=
'這裡替換成檔名'
# 檔名
filename_src=filename+
'.h5'
filename_dst=
'op_'
+filename_src
#獲取grid屬性
# nlay nrow ncol nper timeunits lenunits
with
open
(filename+
'.dis'
,'r'
)as f:
temp=f.readlines(
) dis_name=temp[3]
[2:-
1].split(
' ')
dis=temp[4]
[:-1
].split(
' ')
dis=
list
(map
(int
,dis)
) grid=
dict
(zip
(dis_name,dis)
)print
(grid)
print
(grid[
'ncol'
]*grid[
'nrow'
]*grid[
'nlay'])
#獲取h5檔案資訊
with h5py.file(filename_src,
'r')
as f:
riv=np.array(f[
'river/07. property'])
rivid=np.array(f[
'river/02. cell ids'])
top=np.array(f[
'arrays/top1'])
rch=np.array(f[
'recharge/07. property'])
rch_multi=np.array(f[
'recharge/08. property multiplier'])
spech=np.array(f[
'specified head/07. property'])
spechid=np.array(f[
'specified head/02. cell ids'])
et=np.array(f[
'et/07. property'])
#改變引數 
hk=read_k_ele(filename_src,
0.25
,grid[
'nlay'])
print
(hk)
#建立h5檔案副本，並寫入相關引數 
shutil.copy(filename_src,filename_dst)
f=h5py.file(filename_dst,
'a')
#寫入k
for k in
range
(grid[
'nlay'])
:del f[
'arrays/hk%i'
%(k+1)
] f.create_dataset(
'arrays/hk%i'
%(k+1)
,data=hk[k]
) f.close(
) 
main(
)

以上函式是後者呼叫前者的關係，執行的時候需要在**相同目錄下放入projectname.h5和projectname.dis兩個檔案，執行後產生op_projectname.h5檔案，可替換原來modflow的h5檔案。

對了，別忘了import必要的庫

import numpy as np
import pandas as pd
import h5py
import shtil

參考文獻

kuang, x., jiao, j.j., 2014. an integrated permeability-depth model for earth』s crust: permeability of earth』s crust. geophys. res. lett. 41, 7539–7545.

實現K值隨深度衰減

k值近鄰法python3實現

C 實現單鏈表按k值重新排序的方法

深度改造K3 庫存實現物料等級顏色等多維管理

實現K值隨深度衰減

k值近鄰法python3實現

C 實現單鏈表按k值重新排序的方法

深度改造K3 庫存實現物料等級 顏色等多維管理

相關推薦

深度改造K3 庫存實現物料等級顏色等多維管理