实现K值随深度衰减

做水文地质数值模拟的时候,渗透系数是必不可少的参数,当模型模拟范围达到区域尺度,同时考虑地下水深循环时,深部的岩层在上覆地层的作用力下被压实,地层渗透率会降低,渗透率随深度的增加呈指数式下降,这是渗透率的深度衰减作用。MODFLOW没有描述相应变化的程序包,因此我们将渗透系数从h5文件中提取出来在外部处理,再重新输入回模型。

参考2014年Kuang和Jiao的研究结果,用函数描述渗透系数随深度的变化规律:

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.h5ProjectName.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. https://doi.org/10.1002/2014GL061999

已标记关键词 清除标记
相关推荐
课程简介: 历经半个多月的时间,Debug亲自撸的 “企业员工角色权限管理平台” 终于完成了。正如字面意思,本课程讲解的是一个真正意义上的、企业级的项目实战,主要介绍了企业级应用系统中后端应用权限的管理,其中主要涵盖了六大核心业务模块、十几张数据库表。 其中的核心业务模块主要包括用户模块、部门模块、岗位模块、角色模块、菜单模块和系统日志模块;与此同时,Debug还亲自撸了额外的附属模块,包括字典管理模块、商品分类模块以及考勤管理模块等等,主要是为了更好地巩固相应的技术栈以及企业应用系统业务模块的开发流程! 核心技术栈列表: 得介绍的是,本课程在技术栈层面涵盖了前端和后端的大部分常用技术,包括Spring Boot、Spring MVC、Mybatis、Mybatis-Plus、Shiro(身份认证与资源授权跟会话等等)、Spring AOP、防止XSS攻击、防止SQL注入攻击、过滤器Filter、验证码Kaptcha、热部署插件Devtools、POI、Vue、LayUI、ElementUI、JQuery、HTML、Bootstrap、Freemarker、一键打包部署运行工具Wagon等等,如下图所示: 课程内容与收益: 总的来说,本课程是一门具有很强实践性质的“项目实战”课程,即“企业应用员工角色权限管理平台”,主要介绍了当前企业级应用系统中员工、部门、岗位、角色、权限、菜单以及其他实体模块的管理;其中,还重点讲解了如何基于Shiro的资源授权实现员工-角色-操作权限、员工-角色-数据权限的管理;在课程的最后,还介绍了如何实现一键打包上传部署运行项目等等。如下图所示为本权限管理平台的数据库设计图: 以下为项目整体的运行效果截图: 得一提的是,在本课程中,Debug也向各位小伙伴介绍了如何在企业级应用系统业务模块的开发中,前端到后端再到数据库,最后再到服务器的上线部署运行等流程,如下图所示:
©️2020 CSDN 皮肤主题: 游动-白 设计师:白松林 返回首页