淡水透镜体数学模型的分类与发展探析

淡水透镜体数学模型可以按照空间维数分为一维模型、二维模型、三维模型和零维模型。另外，与机理性数学模型相对应，还有所谓“概念模型”。概念模型是对水文地质单元进行科学概化，根据现场勘测和实验收集的数据，从地下水系统中提取“对模拟目标有用”的概念，按要求组织信息建立的模型，或者为实际地下水系统建立的初始模型。如1942年，Wentworth提出了一个厚透镜体底部贮存（Bottom Storage of Thick Lens）的概念用以解释Honolulu厚透镜体含水层的流动特性。他假设透镜体界面对淡水水头变化的响应有一个“滞后”时间，从而给出了厚透镜体的动力学平衡。然而，他只是叙述，没有用数学公式来表达含水层系统的动力学问题。另一种概念模型是Wentworth 1948年提出的过渡带模型，用以解释淡水和海水间过渡带的产生与扩展。他的解释使用了离散混合的概念。按空间维数划分的一维模型只有一个空间变量，用Ghyben-Herzberg近似关系确定淡水透镜体的厚度。二维模型有两个空间变量，常忽略淡水与海水间的过渡带，而认为淡水-海水间存在一突变界面，以水头即淡水透镜体厚度为目标函数，模型的构建以Ghyben-Herzberg近似为基础；但如将流动和变密度流的传输现象耦合在一起，二维模型也可以模拟界面处由于密度差产生的过渡带，如SUTRA模型就是一个能模拟过渡带的二维模型。三维模型包含三个空间变量，以淡水水头和水中盐浓度为目标函数，可以模拟过渡带的存在，模型的构建以质量守恒和Darcy定律为基础。零维模型则是近年在淡水透镜体实测数据和数值模拟的基础上发展起来的代数模型，根据珊瑚岛的年降雨量、岛屿宽度和水力传导系数等估算淡水透镜体的厚度，用以评估地下淡水资源，使用十分方便。

按照解的形式，还可将淡水透镜体模型分为解析模型和数值模型。解析模型具有封闭形式的解，能清楚地表明变量和重要参数之间的函数关系。这类模型使用简单，易于理解，经过现场测试数据的验证，这类模型可有效地用于地下水的计划管理。如Mink于1980年针对具有突变界面厚透镜体的垂直剖面，开发了一个适用的解析模型。用现场勘测数据进行测试后，该模型可用于估算含水层的持续开采量。随着电子计算机和计算技术的发展，近年使用的淡水透镜体模型大都是数值模型。由于模型中主管流体运动的动力学方程和主管溶质运移的扩散方程都是二阶非线性偏微分方程，在通常的边界条件下难以求得解析解，只能借助计算机和现代计算技术，将求解域（即珊瑚岛）划分为若干计算单元，用有限差分法、有限元法等数值计算方法将主管方程和边界条件离散为各个单元上成立的一组线性代数方程，从而将求解连续区域内满足偏微分方程的函数问题转化为在一系列孤立点上求解满足偏微分方程的函数值的问题。虽然计算工作量大，数值求解也还有一定的困难和局限性，但它容易改变计算条件，具有灵活、经济和限制少、能获得精确可靠的结果等优点，因此，数值模型在淡水透镜体的数学模拟中有广泛的应用。

历史上淡水透镜体的研究从现场勘测和资料的积累开始，包括水文地质测绘、钻探、水井试验等，最初构建的模型主要是定性的概念性模型。虽然也出现过一些简单的计算方程，但并不完善。如1933年，S.T.Hoyt在Hawaii进行了一次现场试验：从一口井中以不同的流量抽水，同时测量相应的水位降深。然后，他构建了一个方程H＝KQn，式中H是水位降深，Q是抽水速率，通过曲线拟合由试验数据求解K和n。S.T.Hoyt的试验本应得到n＝1的结果，这实际上就是Darcy定律。可惜的是，H是在井壁附近测量的，流动为紊流而非层流。因而，Hoyt求得的n值在1与2之间。直到1946年，Wentworth分析和补充了Hoyt的试验数据，剔除了几个较大抽水速率时的测值，剩下的数据基本上是在层流条件下获得的，n值才正好等于1，从而建立了第一个地下水，也是淡水透镜体的数学模型。

淡水透镜体的数学建模其实也主要从这个时期开始。这一时期，M.Muskat和M.K.Hubbert先后证实了Darcy定律用于地下水的正确性。开始主要是一些简单的计算公式。直到20世纪70年代才有了机理性的、较能反映透镜体流动力学特性的数学模型，并随着对透镜体动力学特性认识的深入和计算技术、计算软硬件的发展，淡水透镜体的数学模型日趋完善。1974年，Meyer等人提出了一个地下水模拟模型。这是一个瞬态电流模拟模型，用以模拟Honolulu岛地下水流，计算水头值与咸淡水界面，并对一些特殊的物理现象（如倒锥）进行了讨论。1977年，Larson等人提出了废水回灌模型。这是美国地质调查局组织开发的一个地下含水层动态界面二维模型，用以计算回灌水的分布，模型中的渗透系数用抽水实验求得，模型仅用于薄的淡水透镜体。1979年，Weatcart等人将流体和盐分的输运方程作为一个耦合的动力系统联立求解，用沙箱实验求得标定参数，对一个理想的二维玄武岩含水层进行了模拟计算，提供了系统参数及其验证的有关资料。1981年，Contractor基于有限元法用一个一维模型模拟了关岛的地下水系统，研究过度抽水引起的海水入侵，模拟的时间步长为1个星期。1993年，Giggs研究马绍尔群岛中淡水透镜体时，使用了二维变密度流数字化模型。2002年，J.M.U.Jucson在研究关岛淡水透镜体对回补的响应时，假设淡-盐水界面为一突变界面，应用了二维海水入侵模型（SWIG2D）。2004年，G.C.Hocking研究岛屿的淡水透镜体时，也使用了二维模型来求得淡-盐水界面的稳态解。二维模型的数学基础是Darcy定律和Ghyben-Herzberg近似，该近似认为海水与淡水处于静平衡状态，并假设淡-盐水界面为一突变界面。按照这一近似，理论上淡水透镜体中海平面以下的淡水厚度与海平面以上的淡水厚度之比为一定值，称为“Ghyben-Herzberg比率”。二维模型比较简单，通常由一个水流方程和相应边界条件构成，边界是珊瑚岛与海平面相交的一条线。(www.xing528.com)

二维模型忽略了多孔介质中的弥散作用，边界条件也十分简单，求解容易，但模型有很大的局限性。1995年，Aly I.EL-kadi出版了一本关于地下水模型的论文集，收集了几十位学者的研究成果和论著，对20世纪的地下水模型进行了回顾和总结。Aly在论文集的前言中指出：①预测地下水流所面临的困难来自软件和硬件的能力不足，事实上有效的三维模型还不存在；②21世纪的数学模型要大大提高模拟的可靠性与准确性。

如前所述，淡水透镜体是一种特殊类型的地下水体，从水力学观点看，与大陆地下水相比有两个显著特点：一是大陆地下水底部为一不透水层，用井抽水时，仅在潜水面产生降落漏斗，抽水强度越大，降落漏斗越大；而珊瑚岛礁淡水透镜体漂浮于海水之上，用井抽水时，不仅潜水面会产生降落漏斗，透镜体底部淡-盐水界面还会因降落漏斗引起的水头下降而上升，出现倒锥。二是大陆地下水一般为淡水，水质均一，研究其运动，关注点是水头分布；而珊瑚岛礁淡水透镜体底部与海水相通，海水中的盐分通过弥散作用向淡水输运，在两者之间形成一个过渡带；透镜体中的淡水通过地下径流流向海洋，将渗入的盐分带走，两种作用达到动态平衡时，淡水透镜体便处于相对稳定状态。这样，淡水透镜体的数学模型就与大陆地下水数学模型不同，差别不仅表现在边界条件上，还表现在构建淡水透镜体数学模型时必须考虑弥散作用，将水流处理为变密度流，数学模型由反映水头分布规律的水流方程和反映溶质运移规律的浓度方程及相应的定解条件构成，在这种情况下，模型为三维模型。1997年Larabi等人采用固定网格的有限元法，求解具有自由边界的三维地下水流动，该法将系数矩阵的几何部分和水力传导部分分离出来，从而加快了计算速度，但前提条件是假设地下水的流动仅局限于饱和淡水区域。1999年，F.Ghassemi等基于SALFLOW软件包，建立了一个三维模型来模拟Home岛的淡水透镜体，但该模型存在显著的数值弥散，其结果也仅与二维模型结果吻合，作者指出需要对模型进行改进。2003年，Alex G.Lee.构建了TOUGH2水文地质模型对珊瑚岛礁淡水透镜体进行三维数值模拟，计算得到了淡水透镜体厚度和混合带厚度的控制参数，但模拟时采用了不相溶分界面条件。2005年，Thomas Graf等人采用了修正的数值模型对存在裂隙的多孔介质中三维变密度流动和溶质运移进行模拟，模拟结果与实验室三维测试结果吻合较好，但Thomas Graf等人提出的边界条件非常规则，难以反映真实的多孔介质含水层。Annamaria Mazzia在研究多孔介质中变密度流时，则使用了高阶Godunov混合四面体网格法，将二维的Godunov混合法推广到了三维，并引入一种变时间步长法，从而改进了计算效率。2006年，X.Mao等人介绍了一种新三维模型PHWAT，该模型将化学反应模型与变密度地下水流动和溶质运移模型耦合，用来模拟多种反应物在变密度地下水流中的运移，但X.Mao等人只是在实验室模拟的情况下给出了数值结果，没有对真实的珊瑚岛礁地下水系统进行探讨。同年，Annamaria Mazzia提出了求解多孔介质中含溶质的变密度流的三维混合有限元-有限体积法，采用三维数值模型，将混合杂交有限元格式用于流动方程，而将混合杂交有限元-有限体积时间分裂法用于输送方程，研究结果表明，这一方法在处理对流问题和密度变化引起流场失稳等问题时较为有效。上述部分资料的分析表明，三维模型可以避免二维模型假设引起的失真，使计算结果更准确可靠。因此，近年在淡水透镜体动力学研究中，三维模型获得了广泛的应用。但在一定条件下，二维甚至一维模型仍不失为研究模拟计算淡水透镜体的有用工具。