一种地下水环境的建模及数值模拟方法技术

本发明专利技术涉及一种地下水环境的建模及数值模拟方法，其包括以下步骤：１）应用ＩＤＬ软件的空间插值函数ＫＲＩＧ２Ｄ（）或ＴＲＩ＿ＳＵＲＦ（）函数对地表高程、潜水埋深进行空间网格插值，同时划分空间有限差分数值模拟网格，得到两个相同的差分网格矩阵；设置模拟区土壤属性文件和地表高程与地下水水位属性文件；计算地下水位埋深ｈ０和单元边界流速场；２）建立非饱和带一维溶质运移模型；３）建立饱和带溶质迁移二维数值模型：４）将上式中的各项式进行近似表示：５）利用ＩＤＬ矩阵运算对以上各分量进行求解。本发明专利技术基于ＩＤＬ平台，可高效、快速、简便地模拟地下水环境，因此，可广泛用于地下水环境突发污染事件的数值模拟、预测、应急预测等过程中。

【技术实现步骤摘要】

本专利技术涉及一种数值模拟方法，特别是关于一种地下水环境的建模及数值模拟方法。
技术介绍
地下水环境突发污染事故需要快速掌握污染物在地下水环境中的迁移转化规律，确定污染物扩散范围，进而迅速实施有针对性的治理与防范。由此要求地下水环境数值模型具有快速高效的建模与计算能力。目前，地下水环境的数值模拟方法主要有欧拉法、拉格朗日法、混合欧拉-拉格朗日法等。基于这些方法构建了MT3D、HST3D、WORM、SWIFT、USGS2D-MOC等地下水环境数值模拟软件。这些软件各有其优缺点，但它们或因局限于饱和带或非饱和带单一含水层污染物迁移模拟；或因构建模型程序复杂，难以适应快速建模的需要；或因计算单元较多，计算效率低，难以应用于地下水环境突发污染事故应急数值模拟等，因此，不能对地下水环境突发污染事故进行全面、高效、准确的模拟计算。
技术实现思路
针对上述问题，本专利技术的目的是提供一种地下水环境的建模及数值模拟方法，该方法可对地下水环境的突发污染事故进行全面、高效、准确的模拟计算。为实现上述目的，本专利技术采取以下技术方案：一种地下水环境的建模及数值模拟方法，其包括以下步骤：1)前处理：设置两个属性文件，一是模拟区土壤属性文件，包括研究区边界、X与Y方向网格间距、区域土壤渗透系数、土壤有效孔隙度、污染物在此种土壤中的纵向和横向的弥散度、土壤非饱和带扩散系数和溶质吸附系数；二是地表高程与地下水水位属性文件；应用IDL软件的空间插值函数KRIG2D()或TRI_SURF()对地表高程、地下水位进行空间网格插值，同时划分空间有限差分数值模拟网格，得到两个相同的差分网格矩阵：地表高程矩阵matrix1和地下水位矩阵matrix2；应用SIZE()函数确定研究区网格属性，包括单元网格数、单元网格行列数；应用随机数值生成函数FINDGEN()生成随机数，然后将网格计算单元的节点坐标生成序列数组；2)二维水流模型计算二维水流模型计算包括：溶质非饱和带迁移模拟的地下水埋深矩阵计算以及饱和带二维地下水流场计算部分；利用前处理得到的地表高程矩阵与地下水位矩阵相减得非饱和带厚度h0，即地下水位埋深矩阵matrixh0：h0：matrixh0＝matrix1-matrix2该矩阵直接用于非饱和带溶质迁移模拟计算；饱和带二维溶质迁移模型的水动力条件是地下水流场的分布，该流场由二维水流模型计算得到，应用二维水流模型计算地下水流场：①计算单元节点流速：利用达西定律：v=KIne]]>-->其中，v为流速；K为渗透系数；I为水力坡度；ne为有效孔隙度；②计算单元边界面流速：求得节点流速，根据网格中相邻两个节点的流速，取两个节点流速的平均值，计算得到网格单元边界面上流速分量的IDL矩阵；3)建立非饱和带一维溶质运移模型；t1=ΔHvz---(1)]]>其中，t1为溶质在非饱和带中的运移总时间；ΔH为非饱和带厚度，ΔH＝h0；vz＝D/Δz表示假设非饱和带厚度为Δz＝1m时，污染物垂向扩散速度，D为垂向扩散系数，量纲为L2/T；4)建立饱和带溶质迁移二维数值模型：R∂C∂t1=∂∂x(Dxx∂C∂x)+∂∂x(Dxy∂C∂y)+∂∂y(Dyy∂C∂y)+∂∂y(Dyx∂C∂x)----(2)]]>∂∂x(vxC)-∂∂y(vyC)+qsθCs-λ1C-λ2ρbθC‾]]>其中，R为等温吸附系数，C为溶解浓度，单位为ML-3，t1为饱和带中溶质迁移数值模拟的时间，Dij为弥散系数张量，单位为L2T-1，x、y为计算距离L，vx和vy为地下水流速，单位为LT-1，，qs为源/汇处单位体积含水层流量(，单位为M3T-1，θ为含水层孔隙度，Cs为源/汇项浓度，单位为ML-3，λ1为溶解相的反应速率常数，单位为T-1，λ2为吸附相的反应速率常数，单位为T-1，ρb为空隙介质体积密度，单位为ML-3，为吸附浓度，单位为MM-1；5)将公式(2)中的各项式进行近似表示：在规则间距的节点中心网格中，公式(2)右边第一项在单元(i，j)近似为：∂∂x(Dxx∂C∂x)≈1ΔxDxx(i,j+1/2)(Ci,j+1-Ci,j)-Dxx(i,j-1/2)(Ci,j-Ci,j-1)Δx---(3)]]>公式(3)表示x方向上，由浓度梯度引起从x方向上进入单元(i，j)的净弥散通量；公式(2)右边第二项近似为：∂∂x(Dxy∂C∂y)]]>≈1ΔxDxy(i,j+1/2)(Ci+1,j+1/2-Ci-1,j+1/2)-Dxy(i,j-1/2)(Ci+1,j-1/2-Ci-1,j-1/2)2Δy---(4)]]>=14ΔxΔy[Dxy(i,j+1/2)(Ci+1,j+1+Ci+1,j-Ci-1,j+1-Ci-1,j)-]]>Dxy(i,j-1/2)(Ci+1,j+Ci+1,j-1-Ci-1,j-Ci-1,j-1)]]]>公式(4)表示由y方向的浓度梯度引起的从x方向进入单元(i，j)的净弥散通量；同理，得到由y方向的浓度梯度引起的从y方向上进入单元(i，j)的净弥散通量，即公式(2)右边的第三项为：∂∂y(Dyy∂C∂y)≈1ΔyDyy(i+1/2,j)(Ci+1,j-Ci,j)-Dyy(i-1/2,j)(Ci,j-Ci-1,j)Δy---(5)]]>由x方向的浓度梯度引起的从y方向进入单元(i，j)的净弥散通量，即公式(2)右边的第四项，写为：∂∂y(Dyx∂C∂x)]]>-->≈1ΔyDyx(i+1/2,j)(Ci+1/2,j+1-Ci+1/2,j-1)-Dyx(i-1/2,j)(Ci-1/2,j+1-Ci-1/2,j-1)2Δx---(6)]]>=14ΔyΔx[Dyx(i+1/2,j)(Ci+1,j+1+Ci,j+1-Ci+1,j-1-Ci,j-1)-]]>Dyx(i-1/2,j)(Ci,j+1+Ci-1,j+1-Ci,j-1-Ci-1,j-1)]]]>公式(2)右边第五项为从x方向进入单元(i，j)的净对流通量，近似为：∂∂x(vxC)≈1Δx{本文档来自技高网...

【技术保护点】
一种地下水环境的建模及数值模拟方法，其包括以下步骤：１）前处理：设置两个属性文件，一是模拟区土壤属性文件，包括研究区边界、Ｘ与Ｙ方向网格间距、区域土壤渗透系数、土壤有效孔隙度、污染物在此种土壤中的纵向和横向的弥散度、土壤非饱和带扩散系数和溶质吸附系数；二是地表高程与地下水水位属性文件；应用ＩＤＬ软件的空间插值函数ＫＲＩＧ２Ｄ（）或ＴＲＩ＿ＳＵＲＦ（）对地表高程、地下水位进行空间网格插值，同时划分空间有限差分数值模拟网格，得到两个相同的差分网格矩阵：地表高程矩阵ｍａｔｒｉｘ１和＊ｘ）≈１／ΔｘＤ↓［ｘｘ（ｉ，ｊ＋１／２）］（Ｃ↓［ｉ，ｊ＋１］－Ｃ↓［ｉ，ｊ］）－Ｄ↓［ｘｘ（ｉ，ｊ－１／２）］（Ｃ↓［ｉ，ｊ］－Ｃ↓［ｉ，ｊ－１］）／Δｘ（３）公式（３）表示ｘ方向上，由浓度梯度引起从ｘ方向上进入单元（ｉ，ｊ）的净弥散通量；公式（２）右边第二项近似为：＊／＊ｘ（Ｄ↓［ｘｙ］＊Ｃ／＊ｙ）≈１／ΔｘＤ↓［ｘｙ（ｉ，ｊ＋１／２）］（Ｃ↓［ｉ＋１，ｊ＋１／２］－Ｃ↓［ｉ－１，ｊ＋１／２］）－Ｄ↓［ｘｙ（ｉ，ｊ－１／２）］（Ｃ↓［ｉ＋１，ｊ－１／２］－Ｃ↓［ｉ－１，ｊ－１／２］）／２Δｙ＝１／４ΔｘΔｙ［Ｄ↓［ｘｙ（ｉ，ｊ＋１／２）］（Ｃ↓［ｉ＋１，ｊ＋１］＋Ｃ↓［ｉ＋１，ｊ］－Ｃ↓［ｉ－１，ｊ＋１］－Ｃ↓［ｉ－１，ｊ］）－Ｄ↓［ｘｙ（ｉ，ｊ－１／２）］（Ｃ↓［ｉ＋１，ｊ］＋Ｃ↓［ｉ＋１，ｊ－１］－Ｃ↓［ｉ－１，ｊ］－Ｃ↓［ｉ－１，ｊ－１］）］（４）公式（４）表示由ｙ方向的浓度梯度引起的从ｘ方向进入单元（ｉ，ｊ）的净弥散通量；同理，得到由ｙ方向的浓度梯度引起的从ｙ方向上进入单元（ｉ，ｊ）的净弥散通量，即公式（２）右边的第三项为：＊／＊ｙ（Ｄ↓［ｙｙ］＊Ｃ／＊ｙ）≈１、ΔｙＤ↓［ｙｙ（ｉ＋１／２，ｊ）］（Ｃ↓［ｉ＋１，ｊ］－Ｃ↓［ｉ，ｊ］）－Ｄ↓［ｙｙ（ｉ－１／２，ｊ）］（Ｃ↓［ｉ，ｊ］－Ｃ↓［ｉ－１，ｊ］）／Δｙ（５）由ｘ方向的浓度梯度引起的从ｙ方向进入单元（ｉ，ｊ）的净弥散通量，即公式（２）右边的第四项，写为：＊／＊ｙ（Ｄ↓［ｙｘ］＊Ｃ／＊ｘ）≈１／ΔｙＤ↓［ｙｘ（ｉ＋１／２，ｊ）］（Ｃ↓［ｉ＋１／２，ｊ＋１］－Ｃ↓［ｉ＋１／２，ｊ－１］）－Ｄ↓［ｙｘ（ｉ－１／２，ｊ）］（Ｃ↓［ｉ－１／２，ｊ＋１］－Ｃ↓［ｉ－１／２，ｊ－１］）／２Δｘ＝１／４ΔｙΔｘ［Ｄ↓［ｙｘ（ｉ＋１／２，ｊ）］（Ｃ↓［ｉ＋１，ｊ＋１］＋Ｃ↓［ｉ，ｊ＋１］...

【技术特征摘要】
1.一种地下水环境的建模及数值模拟方法，其包括以下步骤：1)前处理：设置两个属性文件，一是模拟区土壤属性文件，包括研究区边界、X与Y方向网格间距、区域土壤渗透系数、土壤有效孔隙度、污染物在此种土壤中的纵向和横向的弥散度、土壤非饱和带扩散系数和溶质吸附系数；二是地表高程与地下水水位属性文件；应用IDL软件的空间插值函数KRIG2D()或TRI_SURF()对地表高程、地下水位进行空间网格插值，同时划分空间有限差分数值模拟网格，得到两个相同的差分网格矩阵：地表高程矩阵matrix1和地下水位矩阵matrix2；应用SIZE()函数确定研究区网格属性，包括单元网格数、单元网格行列数；应用随机数值生成函数FINDGEN()生成随机数，然后将网格计算单元的节点坐标生成序列数组；2)二维水流模型计算二维水流模型计算包括：溶质非饱和带迁移模拟的地下水埋深矩阵计算以及饱和带二维地下水流场计算部分；利用前处理得到的地表高程矩阵与地下水位矩阵相减得非饱和带厚度h0，即地下水位埋深矩阵matrixh0：h0：matrixh0＝matrix1-matrix2该矩阵直接用于非饱和带溶质迁移模拟计算；饱和带二维溶质迁移模型的水动力条件是地下水流场的分布，该流场由二维水流模型计算得到，应用二维水流模型计算地下水流场：①计算单元节点流速利用达西定律：v=KIne]]>其中，v为流速；K为渗透系数；I为水力坡度；ne为有效孔隙度；②计算单元边界面流速求得节点流速，根据网格中相邻两个节点的流速，取两个节点流速的平均值，计算得到网格单元边界面上流速分量的IDL矩阵；3)建立非饱和带一维溶质运移模型；t1=ΔHvz---(1)]]>其中，t1为溶质在非饱和带中的运移总时间；ΔH为非饱和带厚度，ΔH＝h0；vz＝D/Δz表示假设非饱和带厚度为Δz＝1m时，污染物垂向扩散速度，D为垂向扩散系数，量纲为L2/T；4)建立饱和带溶质迁移二维数值模型：R∂C∂t1=∂∂x(Dxx∂C∂x)+∂∂x(Dxy∂C∂y)+∂∂y(Dyy∂C∂y)+∂∂y(Dyx∂C∂x)----(2)]]>∂∂x(vxC)-∂∂y(vyC)+qsθCs-λ1C-λ2ρbθC‾]]>其中，R为等温吸附系数，C为溶解浓度，单位为ML-3，t1为饱和带中溶质迁移数值模拟的时间，Dij为弥散系数张量，单位为L2T-1，x、y为计算距离L，vx和vy为地下水流速，单位为LT-1，，qs为源/汇处单位体积含水层流量(，单位为M3T-1，θ为含水层孔隙度，Cs为源/汇项浓度，单位为ML-1...

【专利技术属性】
技术研发人员：吴文强，陈求稳，黄国鲜，马金锋，李伟峰，
申请(专利权)人：中国科学院生态环境研究中心，
类型：发明
国别省市：11[中国|北京]