数值法水量计算的步骤

数值法水量计算的正式进行，是在水文地质勘察外业工作完成之后，计算所需的各种原始数据按照数值法水量计算的要求已全部整理、分析、计算完毕，只等调用；计算所需的任何一种数据缺失，都会导致水量计算无法进行。

一般而言，数值法水量计算需要经过如下几个重要的步骤。

1.概化水文地质模型

包括对计算区域的圈定、计算域边界形状及性质的确定、含水层的划分、上下层水力联系的确定、参数分区等。

其中，计算区域的圈定，主要考虑计算任务的需要，兼顾边界条件处理的方便。对于数百、数千、数万甚至数十万平方千米的大区域，多在勘察区边界附近分段选定出方便概化处理的计算区边界即可，计算区边界与勘察区边界两者一般不会相差非常悬殊。然而，对于面积仅有数平方千米、或者最多数十平方千米的供水水源地，计算区边界圈定与水源地范围常常悬殊很大，并且必须比水源地范围要大。根据情况不同，常会作这样的处理：

（1）当近距离内没有已建成的相邻水源地时，计算区域的圈定主要考虑本水源地开采方案布设的方便、边界条件概化的方便。

（2）当近距离内有已建成的相邻水源地、新建水源地开采有可能影响到相邻水源地时，计算区域的圈定除了考虑本水源地开采方案布设的方便、边界条件概化的方便之外，还应当把已建成的相邻水源地也纳入计算域内，计算本水源地布设的开采方案时，把相邻水源地现状的开采方案也纳入计算中，以论证计算本水源地布设的不同开采方案对相邻水源地的影响。

（3）当计算区的某一侧是无限延伸时，也就是说开采方案引起的流场变化远远达不到这一侧的自然边界时，如果也不具备一条已知水头变化规律的“已知变水头的边界”时，可以把这一侧的模型边界向外推至开采方案引不起流场变化的数公里之外，作为已知的定水头边界处理。

2.选用合适的数学模型来描述

是潜水还是承压水还是承压转无压水，或者是二者、三者都有；单层还是多层含水层；各向同性还是各向异性；二维流还是准三维流还是三维流；稳定流还是非稳定流；以及初始条件、边界条件的数学描述。

3.单元剖分

在参数分区的基础上进行单元剖分，单元剖分时应当注意：①剖分的单元不能一个单元骑跨在两个参数区上；②用作拟合孔、预测孔的观测孔要剖分在节点上。

4.编程

自己编制计算程序，或者购买适合的计算软件。如果是自己编制计算程序，则要经过调试程序、检验程序，确认程序无误后，再投入正式使用。如果是购买适合的计算软件，则要经过学习、使用软件的培训。

5.计算时段剖分

先确定模型识别、模型检验所涉及起始日期、时间，结束日期、时间，然后，按照一定的递推关系进行计算时段的剖分。

6.初始流场、末刻流场的模拟

在统测水位绘制各含水层初始流场、末刻流场的基础上，模拟计算出剖分图上各含水层全部节点的初始水位值、末刻水位值，作为数值法水量计算程序所用的初始流场、及对照末刻计算流场的实测流场。

7.模型识别

把所确定的模型识别起始日期到结束日期之间的源汇项数据按时段代入计算程序、把各区的参数初值作为模型的参数初值代入计算程序，运行程序，通过不断调试模型参数得到的拟合孔的计算历时曲线与实测的历时曲线进行拟合，调试出一组仿真度最高的模型参数。

8.模型检验

把所确定的模型检验起始日期到结束日期之间的源汇项数据按时段代入计算程序、把模型识别调试出的模型参数代入计算程序，主要看模型识别结束日期之后拟合孔的计算历时曲线，与实测的历时曲线的吻合程度，来检验模型识别阶段调试出的模型参数是否经得起“外推”的检验；如果“外推”段的拟合曲线误差太大，则应推倒识别阶段的模型参数，重新进行模型识别、模型检验，直到满足误差要求为止。

9.模型预报

经过模型识别、通过模型检验后确定了模型参数，即最终确定了数值模型，即可用其来进行模型预报。进行模型预报时，一般根据供水需要设计多个可供选择的开采方案，分别代入数值模型，来预测未来数十年（城市供水水源地一般不少于30a）后的流场，然后，通过对比、分析，选择其中较合理的开采方案作为水源地建设、施工设计的水文地质依据。

本章以下几节分别介绍不同水文地质条件的几种二维流有限元水量计算问题。

将水盐运移参数加到上述数值方程上即可进行数值模拟。数值模拟包括模型验证和预测预报两个方面。在已知初始条件和边界条件的前提下，模型验证通过以下步骤进行：

（1）根据实测土壤剖面负压h和土壤溶液浓度c，用三次样条插值方法给出剖面上各节点负压和溶液浓度的初始值。

（2）根据气象资料和地表土壤含水率计算蒸发量E。

（3）由根系吸水模型计算根系吸水率Sr。

（4）求解水分方程，给出时段末各节点的土壤负压分布。

（5）由时段末的负压分布，计算土壤孔隙水流速ν。

（6）求解盐分方程，给出时段末各节点的土壤溶液浓度分布。

（7）用实测数据对模型进行检验。

本次模型验证所用资料为1998年4月30日至1999年9月30日，共计518天。时间步长1h，空间步长1cm。计算中所需要的大量数据，如节点初始负压、初始浓度，各时段降雨量、水面蒸发量、地下水位埋深、地下水矿化度等信息，均以数据文件的形式提供。由于三个监测断面的负压由真空表型张力计监测，以kPa表示，所以计算时先将其换算为cm；土壤溶液浓度由盐分传感器监测，以电导率（mS/cm）表示，同样须将其换算为溶质浓度（g/L）。计算的下边界取动边界，随地下水位埋深的变化而变化。与不动水体有关参数的取值，根据文献资料并结合模型调试确定，寅阳1＃、大兴2＃：f=0.975，α=0.005，兴隆沙1＃：f=0.6，α=0.005。

数值计算程序用VB5.0编写，在奔腾机上进行计算。整个计算程序由四个程序模块组成：第一个模块为数据输入模块，第二为求解水分方程模块，第三为求解盐分方程模块，第四为数据输出模块。其中求解盐分方程模块又分为求解可动水体子模块和求解不动水体子模块。

根据描述土壤水盐运移的定解问题，通过数值模拟可以得到土壤盐分运移的动态过程，如果数学模型能够描述实际的物理过程，数值方法可靠，模拟得到的土壤盐分动态过程（模拟值）与实际观测得到的土壤盐分动态过程（实测值）应该完全吻合。

图2.5.3为寅阳1＃模拟值与实测值对比图，由图可见实测值与模拟值拟合相对较好。说明本文所建立的数学模型和提出的数值方法是可行的。这次模型验证，模拟时间较长518天，纵观整个模拟过程，从宏观上来看，模拟值与实测值的动态变化趋势是一致的，并且在模拟过程中没有出现明显的误差累积叠加和扩大的趋势。因此，可以运用所建模型进行有关土壤盐分动态方面的预测预报。

图2.5.3 寅阳1＃模拟值与实测值对比图

---