FEMWATER地下水数值模拟基础

1.建立水文地质概念模型

GMS中FEMWATER使用有限单元法求解地下水流模型及溶质运移模型,具体建模过程与MODFLOW有许多相似之处,但也有需要注意的地方。

(1)确定模拟区范围。FEMWATER导入底图的方式与MODFLOW基本相同。

(2)边界条件的概化。FEMWATER中除了MODFLOW中使用的三类边界条件外,还提供了一种变边界(Variable Boundary)条件。变边界条件通常用于模拟最顶层单元表面过程(如空气-土的交界面),用于描述水分的蒸发、入渗模拟。多数情况下,入渗、蒸发等源汇项通过给定流量边界条件来处理,但该类处理方式要求给定流量方向平行于边界单元的法线方向。事实上降雨、入渗等方式的补给由于顶部单元厚度的起伏而并非如此,该类处理方式将人为加大补给量,FEMWATER的变边界性质能更精确地模拟地表起伏变化条件下的降雨入渗和蒸发过程(施小清和姜蓓蕾,2008)。

(3)源汇项的概化。FEMWATER对于源汇项的概化与MODFLOW和MT3DMS基本相同。

(4)含水层系统结构的概化。FEMWATER采用压力水头而非通常采用的总水头作为因变量,将饱和-非饱和带作为一个整体进行模拟,这不同于以往将饱和带地下水流和非饱和带的水分单独模拟的状况。尽管目前已经建立了众多三维或者准三维的饱和水流模型,但这些模型大多都没有包括严格意义上的潜水面移动边界条件,因此,在实际应用中一般只局限于承压含水层。例如,MODFLOW以及基于MODFLOW的一系列相关程序都采用二类边界,用井流项来近似处理潜水面,即采用流量边界条件近似表示真实的自由面边界,显然这种方法不能真实反映实际的模型条件(施小清和姜蓓蕾,2008)。

FEMWATER模块能够克服这类问题,不同于MODFLOW模块需要指定含水层类型为承压、非承压或承压-非承压,FEMWATER通过对压力水头和模拟层顶面标高的对比,自动识别含水层的类型。当压力水头为负值时,则表明为非饱和带,压力水头为零时,则为潜水面位置。采用压力水头为因变量的方法能很好地刻画和模拟含水系统中潜水面的位置。

(5)水文地质参数分区。FEMWATER中含水层参数的赋值方式与MODFLOW不同。进行饱和带水流数值模拟时,MODFLOW可直接对模拟含水层的渗透系数和贮水系数(或给水度)赋值,而FEMWATER则间接通过设定液体的属性(如密度、黏滞度液体的压缩系数等)和固体骨架属性(如孔隙度、骨架的压缩系数等)来计算得到模拟含水层的水文地质参数。后者较为繁琐,要求建模者对于水文地质的概念有更深的理解和认识,但该方法能更合理的刻画和模拟含水层的特性,特别对于溶质运移和饱和-非饱和带的数值模拟。值得注意的是,进行含水层水文地质调参时,建模者应清楚每个参数对应的含义。

(6)水力特征概化。根据地下水流状态将研究区地下水系统概化为稳定流或非稳定流,一维流、二维平面流或剖面流,准三维流或三维流等。

(7)初始条件概化。初始条件表示压力水头和浓度的初始状态,在FEMWATER中有3种初始条件:cold starts、hot starts和flow solutions。

coldstarts为模拟初期指定的一系列水头值,给定压力水头有两种方式:一种是输入一个常数,FEMWATER会用该值减去节点标高自动计算出压力水头;另一种是输入水头空间变化,需要存在数据集。

对于在某一点、某一段时间收敛很好,然后不收敛的情况,可以通过减少该处节点的时间步长来尝试解决。对该点进行重新模拟,需要使用“hot starts”型初始条件。hot starts文件中包含压力水头、含水率(可选)、速率(可选)、浓度和hot starts时间。

流量法(flow solutions)仅在溶质运移模拟中使用,利用已计算获得的流量模型来定义三维溶质运移模型。流量法与cold starts结合使用需要输入初始浓度值,与hot starts结合使用需要输入非稳定流浓度和hot starts时间。

2.选择合适的模拟程序包

(1)基础程序包。FEMWATER的主界面下有多个选项,分别对应不同的功能,包括模拟问题的类型(水流模型、溶质运移模型、水流和溶质运移模型)、地下水流的性质(稳定流、非稳定流)、计算方法(节点/节点求积分、节点/高斯求积分、高斯/节点求积分和高斯/高斯求积分)、权重因子(Crank-Nicolson central法和向后差分法)、松弛参数(非线性流量方程及运移方程和线性化流量方程及运移方程)等。

(2)求解程序包。FEMWATER中提供的求解方法有3种:逐点迭代矩阵法、预处理共轭梯度法(多项式法)和预处理共轭梯度法(不完全平方法)。(www.xing528.com)

逐点迭代矩阵法包括高斯-赛德尔迭代法、逐次低松弛迭代法和逐次超松弛迭代法。在矩阵对角占优时,逐点迭代矩阵法的结果收敛。在使用逐点迭代法能满足计算效率的情况下,应尽量采用该方法,只有当收敛速度很慢时才需要考虑用另外两种方法。

预处理共轭梯度法(多项式法)采用共轭梯度法解决矩阵方程,它将多项式作为预处理因子,当矩阵是对称正定矩阵时,结果收敛。理论上,收敛速度比逐点迭代法快。只有当逐点迭代法太慢时才使用该法。

预处理共轭梯度法(不完全平方根法)采用共轭梯度法解决矩阵方程,它将不完全Choleski分解作为预处理因子,当矩阵是对称正定矩阵时,结果收敛。然而,当矩阵些微不对称,也可产生收敛解。理论上,收敛速度比逐点迭代法快,与多项式共轭梯度法相当,只有当逐点迭代法太慢时才使用该法。它通常比多项式共轭梯度法使用广泛。

(3)可选程序包。GMS选项卡中将外部源汇项(边界条件)和内部源汇项放在一起供用户选择。FEMWATER中可定义节点边界条件、面边界条件、点源汇项边界条件(井)等。应当根据模拟区的实际水文地质条件选择相应的程序包。

3.数值模拟模型的建立

利用有限单元法对模型进行时间剖分和空间离散。时间剖分和空间离散的原则与MODFLOW与MT3DMS相同。但FEMWATER在空间离散方面有着不同的特点:

(1)不同于MODFLOW和MT3DMS生成的正交单元网格,FEMWATER支持的单元类型有六面体、四面体、三棱柱。使用者可以很方便地生成各种形状的有限元网格,在建模中使用三棱柱比较多,因为三棱柱有一些好的性质便于建模。

(2)FEMWATER提供4种创建网格方法:TINs法、Boreholes法、Solid法和网格节点法。

(3)MODFLOW和MT3DMS中对某一单元进行加密时,将对单元格所在的行或列全部加密,而FEMWATER可只对特定的单元加密,而不影响相邻单元。

FEMWATER提供3种加密方法:垂向列加密只在X、Y方向分割六面体和楔体;所有单元格加密为三角形,包括粗糙加密和细致加密方式;保留单元格类型加密,使单元格转化为更小的同类型单元格。

从计算效率上看,单元数目越少并不意味着计算速度越快。单元越少,计算时所需要的迭代次数越多,既没达到精确的计算目的,计算速度也不一定快。因此最好尽可能用足够的单元来刻画模拟区。

基于GMS建立水文地质概念模型之后,将概念模型中定义的属性数据转换到所剖分的网格单元和结点上,即将概念模型转化为数值模拟模型,通过GMS来实现对模拟区的计算。

4.模型运行

FEMWATER的模型检查和模型运行方式与MODFLOW和MT3DMS基本相同。