CN102567634A - 一种基于水循环的地下水数值仿真方法 - Google Patents

Abstract

本发明涉及一种基于水循环模拟的地下水数值仿真方法，所述方法的步骤如下：建立水循环模拟系统的步骤；建立地下水数值仿真系统的步骤；以日尺度为时间步长，通过水循环模拟方法和地下水数值仿真系统在时间步长内的双向信息交互，构成两者之间的同步耦合，从而形成能够适合大尺度地下水数值仿真应用的解决方案。本发明利用基于子流域空间离散的水循环模拟，实现了在大空间尺度流域或区域和长仿真期条件下应用的能力。本发明使用有限差分法，使水循环模拟与地下水数值仿真过程具有双向反馈的能力。本发明还利用网格单元与子流域的面积直接叠加的方式，提高了数据处理运行效率。本发明所述方法技术方法通用，易于推广应用。

Description

一种基于水循环的地下水数值仿真方法

技术领域

本发明涉及一种基于水循环的地下水数值仿真方法，是一种对水文水资源的数据信息进行快速处理方法，是一种通过模拟和建立数学模型对水文资源的数据信息进行处理的方法。

背景技术

关于地下水数值仿真技术及其系统，目前技术主要有以下几种。一种是单纯的地下水数值仿真技术及系统，仅仿真地下水自身循环过程，不考虑地下水循环与外部水分循环条件的信息交互与相互影响。此类仿真技术及系统很多，代表仿真技术及相关系统有Visual Modflow、FEFLOW、Visual Groundwater等等。二是文件交换形式的水循环模拟与地下水数值仿真技术耦合方法。先用水循环模拟计算出地下水数值仿真所需的前期数据信息，再将数据信息处理成符合地下水数值仿真要求的数据文件格式，最后地下水数值仿真系统读入上述数据文件完成仿真过程。该技术属于松散数据耦合式的解决方案。三是网格式交互的水循环模拟与地下水数值仿真技术方法，这类技术方法代表有MIKE-SHE、IGSM、MODHMS等。主要技术关键是将水循环模拟时的网格单元与地下水数值仿真时的网格单元构成严格的一一对应关系，通过每个网格单元内数据的同步交互，可实现水循环模拟与地下水数值仿真的统一。该技术方案是目前本领域内相对较为先进的代表技术。

以上三种现存的地下水数值仿真技术都不同程度存在若干不足。

对于单纯的地下水数值仿真技术及系统，缺陷源于仅从地下水自身循环的观点看待地下水，基本上不考虑地下水与土壤水、地表水之间水分转化的相互影响效应，造成一些地下水仿真所需的关键数据，如降水入渗补给量、灌溉渗漏量补给量、河道渗漏补给量等，只能在地下水数值仿真系统外部显式输入。在较小时空尺度地下水数值仿真分析时，如堤防渗流、矿井疏干排水、地下水水源地评价等，由于这些补给量数据在地下水流动循环中不占主导地位，因此对仿真可靠度的影响较小。但在地下水数值仿真的空间尺度较大（如区域或流域）、仿真期较长（如多年状况）时，这些数据是地下水流动循环的主要通量，此时数据的合理性和精度是影响地下水数值仿真可靠性的重要因素。由于下垫面条件和地表岩性参数的复杂性、气象变化和人类活动等因素的影响，直接确定这些数据十分困难，精度也难以保证。

由于单纯地下水数值仿真技术在大尺度地下水问题应用上的不足，现有技术中发展了文件交换形式的水循环模拟与地下水数值仿真技术耦合方法，试图改善这一应用问题。该方法既有优点也有缺点。优点在于通过水循环模拟将土壤水、地表水对地下水的动态影响纳入到地下水数值仿真过程中，是对原技术方法的改进和拓展；同时由于仅是水循环模拟系统与地下水数值仿真系统之间文件交换形式的松散耦合，该技术方法的灵活性较大，也比较容易实现。可以选取不同复杂程度的水循环模拟系统与地下水数值仿真系统进行配合，只需要将水循环模拟系统的输出按地下水数值仿真系统的输入要求改造成相应的文件格式即可。同时地下水数值仿真系统也可以选用不同的方案，如有限差分形式、有限元形式等。缺点在于若研究的地下水问题较为复杂，时间步长要求较小，则交换文件有可能十分庞大，一方面提高了对计算系统存储容量的要求，另一方面庞大的交换文件也会引起传输之间的不便，同时文件交换还会显著影响地下水数值仿真的工作效率。另外一个本质的不足是，通常这种技术方法只能实现从水循环模拟到地下水数值仿真的单向数据信息传递，地下水数值仿真的数据信息无法同步反馈到水循环模拟过程中实现双向作用过程。水循环模拟系统的优势在于可以模拟大气水-土壤水-地表水-地下水的一体化过程，但地下水循环部分多数以均衡模式（水桶模式）处理，缺乏模拟地下水侧向流动过程的能力，而这正是地下水数值仿真系统的优势。若能将地下水数值仿真系统的输出同步反馈到水循环模拟系统，则可以改善水循环模拟的精度，进而其输出又可提高地下水数值仿真可靠度和效果，真正实现两个系统之间的优势互补。但文件交换形式的水循环模拟与地下水数值仿真技术耦合方法从技术上讲很难实现双向反馈，因此虽然对流域/区域尺度地下水数值仿真有一定改善，但程度有限。

网格式交互的水循环模拟与地下水数值仿真技术方法正是针对以上文件交换形式方法的不足而提出的。这种方法将水循环模拟与地下水数值仿真融为一个系统，无需数据文件的交换，所有数据信息的传递在内存中进行。同时最重要的是在网格单元的基础上实现了水循环模拟和地下水数值仿真过程的统一，具有真实的双向数据信息反馈，因此该技术方法相对而言比较先进。目前的主要不足在于两方面，一是对水循环模拟和地下水数值仿真技术方法的要求比较严格，如两者都必须基于矩形网格，而且必须共享同一网格单元剖分。问题在于网格单元尺度太大时水循环模拟将产生明显的尺度效应，影响模拟精度，应用过程中需要将网格单元控制在较小尺度范围内。对于面积较大的流域或区域，建模时网格单元的规模将十分庞大，导致运行十分耗时，这样对硬件系统的存储能力和计算能力要求都很高。同时由于网格单元是地下水数值仿真系统和水循环模拟系统共享的，庞大的网格单元规模也会严重影响地下水数值仿真系统的运行效率。二是基于网格单元的水循环模拟系统虽然物理机制较强，但一般结构都比较复杂、需要大量参数和数据支撑，专业性很强，不易被一般用户掌握。基于以上原因，网格式交互的水循环模拟与地下水数值仿真技术方法虽然比较先进，但仍然只适合在较小尺度的流域/区域上应用，同时不易推广。

发明内容

为了克服现有技术的问题，本发明提出了一种新的基于水循环模拟的地下水数值仿真方法，所述的方法具有在大空间尺度流域/区域和长仿真期条件下应用的能力，并具有进行水循环模拟与地下水数值仿真过程的双向反馈能力，在保证精度的基础上具有较高的运行效率。

本发明的目的是这样实现的：一种基于水循环模拟的地下水数值仿真方法，所述方法的步骤如下：

建立水循环模拟系统的步骤：用于采用基于子流域空间离散模型的水循环模拟方法，以确定区域范围内地下水的垂向循环通量的时空分布，从而产生地下水数值仿真需求的数据信息，包括降水入渗补给量、河道/水库/湿地/渠系等地表水体渗漏补给量、井灌回归补给量、地下水基流量、潜水蒸发量、地下水开采量；

建立地下水数值仿真系统的步骤：用于采用基于矩形网格单元空间离散技术的地下水数值仿真系统，模拟地下水在研究区域内的侧向流动，其仿真原理基于以下三维地下水动力学方程：

              （1）

其中

，

和为渗透系数在X, Y和Z方向上的分量，量纲为(LT^-1)；

为水头(L)；

为单位体积流量(T^-1)，用以代表来自源汇处的水量；为孔隙介质的贮水率((L^-1)；

为时间(T)，使用网格单元划分法对被模拟区域进行空间划分后，对被模拟区域行数学模拟；

进行水循环模拟系统与地下水数值仿真系统空间嵌套处理的步骤：用于在水循环模拟系统所划分的子流域与地下水数值仿真系统所划分的网格单元之间建立空间关联。该步骤将子流域与网格单元进行空间叠加，确定子流域和网格单元之间的从属关系，并计算各网格单元在各子流域中所占的面积比例。空间嵌套处理将为水循环模拟系统和地下水数值仿真系统在同步耦合过程中的双向数据信息交互提供基础；

水循环模拟系统与地下水数值仿真系统融合的步骤：用于以日尺度为时间步长，通过水循环模拟方法和地下水数值仿真系统在每个时间步长内的双向信息交互，构成两者之间的同步耦合；所述的“同步耦合”指在同一个时间步长内，两个系统完成信息之间的交互反馈；时间步长是进行模拟或仿真时的时间片段，整个模拟期由多个时间步长组成，对每个时间步长两个系统均进行同步耦合，直至模拟结束；所述的“双向信息交互”是以水循环模拟过程中使用的子流域和地下水数值仿真系统中使用的网格单元之间的空间嵌套处理作为交互基础，双向信息交互内容包括两方面：一是水循环模拟系统得出的垂向循环通量时空分布信息传递给地下水数值仿真系统，提供地下水仿真所需的源汇处的水量。垂向循环通量时空分布信息包括降水入渗补给量、河道/水库/湿地/渠系等地表水体渗漏补给量、井灌回归补给量、地下水基流量、潜水蒸发量、地下水开采量；二是地下水数值仿真系统得出的地下水位和地下水埋深的信息传递给水循环模拟系统，辅助水循环模拟系统计算地下水与地表水、土壤水之间的转化量。

本发明产生的有益效果是：本发明利用基于子流域空间离散的水循环模拟系统方法确定地下水垂向循环通量的时空分布，为地下水数值仿真提供了在大空间尺度流域/区域和长仿真期条件下应用的能力。本发明利用在内存中进行信息同步耦合的步骤，使水循环模拟与地下水数值仿真过程具有双向信息反馈的能力。本发明还利用网格单元与子流域的面积直接叠加的方法，提高了数据处理运行效率。本发明所述方法技术方法通用，易于推广应用。

 

附图说明

下面结合附图和实施例对本发明作进一步说明。

图1是基于网格的空间离散技术示意图；

图2是基于地形单元的空间离散技术示意图；

图3是基于子流域的空间离散技术示意图；

图4是水循环模拟系统描述的流域/区域水循环转化路径(一般情况)的示意图；

图5是局部加密的有限差分网格示意图；

图6是局部加密的有限元网格剖分示意图；

图7是三角有限元与子流域相容性示意图；

图8是单元中心格式示意图；

图9是节点中心格式示意图；

图10 是子流域与网格单元的空间嵌套图；

图11 是水循环模拟系统与地下水数值仿真系统的同步耦合流程图；

图12 是邯郸市子流域划分与河道分布图；

图13 是邯郸市平原区子流域与网格单元空间嵌套图；

图14 是2007年邯郸市平原区地下水数值仿真效果图。

具体实施方式

实施例一：

本实施例是一种基于水循环模拟的大尺度地下水数值仿真方法。本实施例所述的方法可以编制为计算机程序，运行在PC机或其他通用计算机中。本实施例所述方法的步骤如下：

1)   建立水循环模拟系统的步骤：用于采用基于子流域空间离散模型的水循环模拟方法，用以确定大尺度区域范围内地下水的垂向循环通量的时空分布，从而产生地下水数值仿真需求的数据信息；

2)   建立地下水数值仿真系统的步骤：用于采用基于矩形网格单元空间离散技术的地下水数值仿真系统，模拟地下水在研究区域内的侧向流动，其仿真原理基于以下三维地下水动力学方程：

              （1）

其中

，

和

为渗透系数在X, Y和Z方向上的分量，量纲为(LT^-1)；

为水头(L)；    

为单位体积流量(T^-1)，用以代表来自源汇处的水量；

为孔隙介质的贮水率(L^-1)；

为时间(T)，使用网格单元划分法对被模拟区域进行空间划分后，对被模拟区域进行数学模拟；

3)  水循环模拟系统与地下水数值仿真系统融合的步骤：用于以日尺度为时间步长，通过水循环模拟方法和地下水数值仿真系统在每个时间步长内的双向信息交互，构成两者之间的同步耦合，从而形成能够适合大尺度地下水数值仿真应用的解决方案；所述的“双向信息交互”是以水循环模拟过程中使用的子流域和地下水数值仿真系统中使用的网格单元之间的空间嵌套处理作为交互基础，交互过程在内存中进行；双向信息交互内容包括两方面：一是水循环模拟系统得出的垂向循环通量时空分布信息传递给地下水数值仿真系统，提供地下水仿真所需的源汇处的水量，垂向循环通量时空分布信息包括降水入渗补给量、河道/水库/湿地/渠系等地表水体渗漏补给量、井灌回归补给量、地下水基流量、潜水蒸发量、地下水开采量等；二是地下水数值仿真系统得出的地下水位和地下水埋深的信息传递给水循环模拟系统，辅助水循环模拟系统计算地下水与地表水、土壤水之间的转化量。

本实施例采用了一种不同于网格式交互的新形式，实现了水循环模拟系统与地下水数值仿真系统的融合并能适用于大尺度流域/区域。在技术内容方面包括以下部分。一是水循环模拟系统的开发，二是地下水数值仿真系统的开发，三是水循环模拟系统和地下水数值仿真系统的融合，以下分别进行详述。其中单独的水循环模拟系统和地下水数值仿真系统目前比较多见，这里只作为本技术方案实施的基础条件，水循环模拟系统和地下水数值仿真系统的融合方法是本发明的重点所在。

水循环模拟系统：当前的水循环模拟系统在空间离散技术上主要分为三类，第一类是基于网格单元的空间离散技术，如图1所示。在图1中，标记11表示坡面流，标记12表示河道汇流，标记13表示支流河道，标记14表示主河道。这类水循环模拟系统包括SHE、WEP等；第二类是基于地形单元的空间离散技术，如图2所示。在图2中，标记21表示流域出口，标记22表示地形单元，标记23表示河道，标记24表示流域，标记25表示由地形单元构成的流带，标记26表示地形单元的汇流方向。这类水循环模拟系统包括TOPMODEL，GBHM等；第三类是基于子流域的空间离散技术，如图3所示。在图3中，标记31表示根据分水岭划分的子流域，标记32表示子流域中的河道。第一类和第二类技术出现比较早，研究历史比较长。对流域/区域空间分布的描述分别基于网格单元和地形单元，需要给出每个网格单元或地形单元上的参数。其特点是方法比较严格，物理性比较强，但因尺度效应问题对单元大小比较敏感，同时数据资料要求高，在流域/区域尺度太大时运行成本比较高，工作效率比较低，适合于在数据资料密度集中且尺度较小的流域/区域上应用。第三类技术是近期在前两种技术的基础上发展起来的，融合了分布式和集总式的特点。其主要思想是先根据DEM（数字高程网格）分析汇流路径，再分离出空间上相对独立但具有河道水力联系的子流域，流域/区域的分布特征由子流域体现。在子流域的内部，则继续根据土地利用方式、土壤性质、水土管理的不同进行汇总并划分为多个最小模拟单元。这些最小模拟单元在子流域内部不考虑空间分布，具有集总式的特征。相对而言，基于子流域的空间离散技术比较先进和灵活。比如由于汇流路径的形成只取决于DEM分辨率，因此在同一DEM的基础上，子流域划分的大小并不影响整体汇流路径格局，因此在子流域划分得较大时尺度效应不明显，所以对于大空间尺度流域/区域的模拟也能适应，突破了之前基于网格的空间离散技术对于网格单元尺度的限制。再如对子流域内部进行集中式处理，不仅适当降低了对数据资料的要求，同时简化了一些对模拟结果影响不大的流程，因此运行效率很高。由于这些优势，基于子流域空间离散技术的水循环模拟系统发展很快，目前在应用推广方面已经远远超过前两种离散技术的系统。这类水循环模拟系统包括SWAT、新安江、HSPF等。

通过以上分析对比，为实现本实施例的第一个目标，即具有在大空间尺度流域/区域和长仿真期条件下应用的能力，本实施例涉及的水循环模拟系统需要基于子流域空间离散技术。除此之外，水循环模拟系统还至少应能在子流域层次上对大气水、土壤水、地表水、地下水之间的转化关系有较为清晰和合理的计算流程，以产生地下水数值仿真需求的数据，如降水入渗补给量、河道/水库/湿地/渠系等地表水体渗漏补给量、井灌回归补给量、地下水基流量、潜水蒸发量、地下水开采量等。虽然目前基于子流域空间离散技术的水循环模拟系统很多，但其内部结构和计算复杂程度不一而论，功能强弱也会有一些区别，不过对于具体的某个水循环模拟系统是否合适与地下水数值仿真系统进行融合，也要视乎应用需求而定。一般来说对以上的数据项考虑得越详细，说明该水循环模拟系统对地下水循环转化机制的表达能力越强，相应在可靠性和精度上越好。为不失一般性，本实施例通过以下框图表达了通常情况下水循环模拟系统所描述的流域/区域水循环转化路径，见图4。

图4以一种示意的方式描述了水分在流域/区域中的循环转化过程，其中的水循环陆面过程主要在子流域内部进行，描述水分从降水开始直到循环转化成为河道地表水的过程；水循环河网系统过程则用于描述河道地表水通过河网系统逐级演进直到流域出口的过程。需要说明的是，该图仅代表一般情况，有些水循环模拟系统可能模拟的过程比该图中更加复杂，也有可能更为简化。图中关于地下水部分的循环转化项以阴影的形式标出，这些数据项是地下水数值仿真系统在进行流域/区域尺度应用时所需的输入数据，同时也是水循环模拟系统和地下水数值仿真系统在进行融合时的交互基础。

地下水数值仿真系统：

无论哪种地下水数值仿真系统，都是基于相同的三维地下水动力学方程的：

              （1）

其中：

   

，

和

为渗透系数在X, Y和Z方向上的分量，量纲为(LT^-1)；

   ：水头(L);

   

：单位体积流量(T^-1)，用以代表来自源汇处的水量;

   

：孔隙介质的贮水率((L^-1);

   

：时间(T)。

单独的三维地下水数值仿真系统包括有限差分、有限元、有限体积等几种方法，其中有限差分和有限元法应用较多。有限体积法为有限差分法和有限元法之间的产物，由于相对前两种方法优势不明显，同时实现起来比较复杂，目前应用很少。从表面上比较来看，有限差分法的优势在于概念比较直观，编程容易实现，但缺点在于剖分方式不够灵活，如进行空间剖分时需要以矩形单元为剖分基础，因此在描述非规则区域/流域时会形成部分无效单元，同时在对模拟范围内的部分重要地点（如开采井附近）进行加密剖分时，同行同列的其他地点也会被加密，如图5所示。在图5中，标记51表示无效单元区，标记52表示有效单元，标记53表示开采井。有限元法的模拟原理稍复杂，编程难度相对较大。但其优势在于剖分方式相对灵活，通常以三角形单元剖分的方式对空间进行离散，因此对各种形状的模拟区域均能充分适应，不会形成无效单元，同时很容易对重要地点进行加密（图6）而不影响其他部位。从计算效率上看，因为有限元法没有无效单元，似乎给人一种效率应该比差分法高的错觉，但实际上差分法由于计算步骤相对简单，特别是由于单元编号规整，使得最终形成的大型稀疏矩阵方程比较利于求解，因此同等剖分密度下计算效率比有限元法要高很多，完全可以弥补无效单元带来的负担。有限元法的计算效率反而对计算节点编号的优化和矩阵方程的求解方法有较大的依赖性。最后虽然有限元法在数学理论的推导上更具吸引力，但在仿真效果和精度上与有限差分进行比较，两者并没有明显意义的区别。

对于本实施例明而言，由于要实现水循环模拟系统与地下水数值仿真系统之间的融合，同时要实现双向同步反馈，对于地下水数值仿真系统主要有两点要求。一是要易于与基于子流域空间离散技术的水循环模拟系统相结合，这里所说的结合主要为无形状规则的子流域与地下水数值仿真系统中的剖分单元之间的相容。相比较而言，由于有限元法的剖分单元也是无规则的，似乎更容易与无形状规则的子流域相容，但这很难做到。以三角单元为例，无论剖分如何细密，在子流域之间的边界处仍会有一定的误差，见图7。在图7中，标记71表示子流域，标记72表示三角单元，标记73表示三角单元与子流域之间的边界误差。同时在多个子流域进行联合剖分时三角单元的均匀性难以保证，可能出现某些子流域剖分得过粗，某些又过细，造成仿真精度在区域空间上的非均匀性，这对于区域尺度的应用不利。因此在相容性方面同比有限差分网格剖分方法并无明显优势。二是要有较高的运行效率。这主要是因为水循环模拟系统的计算尺度较小，一般为日尺度或日内尺度，要实现双向同步反馈，地下水数值仿真过程也将以这样的小时间尺度运行。由于地下水数值仿真过程需要反复迭代计算，小时间尺度条件下模拟期内运行工作量较大。比较而言这方面差分方法有较大的优势。

综上比较，在技术方案上，为适应需求，对于地下水数值仿真系统，本发明采取的是有限差分方法。有限差分在形式上可分为单元中心格式，见图8，和节点中心格式，见图9。在图8中，标记81表示网格单元中心的位置，标记82表示网格单元的边界（即网格线），标记83表示网格单元的面积。在图9中，标记91表示计算节点所在的位置，标记92表示计算节点控制的矩形区边界，标记93表示网格线。两种方法都有优缺点，主要体现在于不同类型的渗流边界的处理上。但由于单元中心格式对边界的处理要相对容易，因此本技术方案采用单元中心的格式。

有限差分方法在效率上无疑是较优的选择，但主要需要解决的问题是如何设计一种方法使得地下水网格单元与无规则的水循环子流域有机相结合，这将在下面进行详述。

水循环模拟系统和地下水数值仿真系统的融合方法（空间嵌套和同步耦合）：

以上通过水循环模拟系统和地下水数值仿真系统的详细类比，分析以上各种系统的不同技术特点和优劣，确定了为达到本发明的效果，水循环模拟系统和地下水数值仿真系统的合理选取问题，即水循环模拟系统为基于子流域离散技术，地下水数值仿真系统为单元中心差分格式。最后需解决的问题是这两个系统之间的融合技术问题，包括以下两个技术要点，一是子流域如何与地下水数值仿真系统的网格单元相容；二是信息如何在两个系统之间如何同步相互反馈。

子流域与地下水数值仿真系统的网格单元相容（空间嵌套）：

如前所述，目前已经存在一种技术，可以通过网格式交互的方法实现水循环模拟与地下水数值仿真技术的融合，其技术基础在于水循环模拟系统和地下水数值仿真系统必须共享相同的空间离散网格，但由于尺度效应、运行成本、效率等问题使得该方法难以在大尺度流域/区域进行应用。本发明中水循环模拟系统和地下水数值仿真系统分别基于不同的空间离散方法，面临的首个要解决的问题是子流域与网格的相容技术。大尺度流域/区域应用情况下一般子流域的个数远比网格单元的数量少，因此存在单个子流域能够容纳多个网格单元的条件，这样通常的思路是通过多个网格单元的组合去近似子流域的形状。但难点一是由于子流域是根据DEM汇流分析得出的，因此在子流域边界形态上是无规则的，近似时难免产生误差；二是子流域的面积也是不固定的，有可能大多数子流域的面积比网格单元大，但也存在个别子流域比网格单元还小的情况。以上原因使得用多个网格单元的组合去近似子流域的思路有可能不通用。为此，本实施例采取另外一种结合策略，不强制让网格单元去近似子流域，只需要两者之间的叠加效果即可，如图10所示。图10中标记101表示子流域的编号，标记102表示与子流域边界相交的网格单元，标记103表示子流域内部的网格单元，标记104表示子流域的边界。

这种策略通过子流域的边界切分网格单元，可将网格单元分为两类，一类是完全位于某个子流域内部的网格单元，第二类是位于两个或多个子流域的边界上的网格单元。切分的主要目的是确定网格单元的从属，即第一类网格单元具有唯一从属的子流域，第二类则从属于多个子流域。无论是那类网格单元，其在不同子流域的面积比例是可以确定的。该策略的优点是使网格单元的空间离散独立于子流域离散之外，即子流域对网格单元的空间离散没有硬性要求，因此比较灵活和通用，这样网格单元可以充分根据精度和运行效率要求自由选择剖分尺度。

另外，山丘区地下水和平原区地下水的性质存在一些显著的差别。山丘区子流域由自然分水岭分隔，地下水运动主要发生在风化岩体裂隙中，无明显意义的地下水含水层，一般无需进行地下水数值仿真模拟。平原区坡度平缓，含水层厚度大，各子流域之间地下水有明显的水力联系，是进行地下水数值仿真模拟的重点区。采用本发明提出的空间嵌套策略，还可以根据山丘区和平原区的分布范围指定需要进行地下水数值模拟的子流域。与这些子流域没有叠加关系的网格单元都是无效的网格单元，不参与仿真计算。

信息在水循环模拟系统和地下水数值仿真系统之间的同步相互反馈：

解决信息在水循环模拟系统和地下水数值仿真系统之间的反馈机制需要考虑时间步长问题。时间步长是进行模拟或仿真时的时间片段，模型逐个时间片段进行模拟，直至到达模拟期结束时刻。时间步长的大小取决于模型本身的限制以及计算精度的要求。在进行水循环模拟时，一些水文过程如土壤水含水量变化、蒸发蒸腾等日内变化比较明显，为获取合理的模拟结果多数水循环模拟系统通常以日作为时间步长。相对而言地下水的运动相对缓慢，其响应外界变化所需的时间较长，一般进行仿真时用旬、月作为时间步长对仿真精度影响也不大。当然小时间步长也可以，只不过需要较多的模型运行时间。因此水循环模拟系统和地下水数值仿真系统进行融合时，时间步长主要由水循环模拟系统决定，即以日作为两个系统的共同时间步长。

同步耦合指在同一个时间步长内，两个系统完成信息之间的交互反馈。对于水循环模拟系统而言，是将该时间步长内模拟得出的地下水源汇项信息反馈给地下水数值仿真系统；对于地下水数值仿真系统而言，是在获得水循环模拟系统所反馈的信息后，进行该时间步长内的地下水数值仿真模拟，并将仿真得出的地下水位和地下水埋深数据反馈给水循环模拟系统。地下水位和地下水埋深是水循环模拟系统所需的重要数据，河道/水库/湿地与地下水之间的交换量、土壤水深层渗漏、潜水蒸发等水循环过程的模拟都将用到地下水位和地下水埋深数据。在获得这些由地下水数值仿真系统反馈的数据后，下一个时间步长水循环模拟系统将可以使用这些数据。两个系统之间同步耦合的步骤可描述如下(图11)：

1.任务开始。在准备好相关数据后，开始基于水循环模拟的地下水数值仿真工作。

2.模拟之前的初始化。水循环模拟系统和地下水数值仿真系统独自进行模拟之前各自系统内部的初始化，包括开辟内存空间，读入各自所需的输入数据等。

3.模拟期开始。根据输入数据的信息，设置模拟的起始日期，结束日期等。

4.日时间步长模拟。从模拟期的第1天开始，以日为时间步长，两个系统将开始逐日进行模拟，并进行同步耦合。

5.水循环模拟系统进行水量转化模拟。在每日模拟的开始阶段，水循环模拟系统首先进行水循环模拟，对模拟区域内的水量转化过程进行计算，得出该当日内模拟区域的各项水循环通量。

6.反馈地下水循环通量数据给地下水数值仿真系统。水循环模拟系统将循环通量信息中的地下水部分传递给地下水数值仿真系统。

 7.地下水数值仿真系统设置源汇项。地下水数值仿真系统获得信息后，根据子流域和网格单元间的嵌套关系设置仿真模拟时各个地下水网格的源汇项。

8.建立三维地下水数值仿真矩阵方程。地下水数值仿真系统根据源汇项情况和地下水三维动力学方程原理构建待求解的矩阵方程。

9.近似求解矩阵方程。用求解大型矩阵方程的算法对方程进行近似求解。这些算法包括强隐式法(SIP)、共轭梯度法(PCG)等。

 10.判断收敛。由于是近似求解，有可能地下水位的求解结果并不能满足事先规定的精度要求，因此需要进行判断。如果求解结果满足精度要求，则进入下一步地下水均衡分析；若不满足精度要求，则返回建立三维地下水数值仿真矩阵方程这一步反复迭代，直至满足精度为止。

 11.地下水均衡分析。地下水数值仿真系统对当日的地下水平衡状况进行分析。

 12.保存地下水位和地下水埋深结果。地下水数值仿真模型保存当日求解得到的地下水位和地下水埋深数据。

 13.地下水数值仿真系统反馈地下水位和地下水埋深数据。地下水数值仿真系统将求解得出的地下水位和地下水埋深数据信息反馈给水循环模拟系统。

 14.水循环模拟系统根据反馈数据更新各子流域地下水位和地下水埋深。水循环模拟系统按子流域和网格单元的嵌套关系对获得的地下水位和地下水埋深数据信息进行处理，更新各个子流域的地下水位和地下水埋深数据，作为下一个时间步长（即下一日）初始时水循环模拟时的参数。

 15.保存结果并输出。水循环模拟系统和地下水数值仿真系统各自保存当日的计算结果以及其他状态变量，并根据需要输出结果信息。

16.模拟期结束。一个时间步长的计算任务结束，判断当日是否为模拟期的最后一日。如果是最后一日，模拟期的任务全部完成，进入下一步退出；如果不是最后一日。将当前的模拟日期向前推进1日，继续以上同步耦合模拟过程。

17.退出。所有模拟任务完成，进行一些善后维护工作，如释放内存等。

本实施例应用实例：

通辽市域全面积为5.9万km2，区域面积尺度范围较大。水循环模拟系统中，全通辽市根据DEM汇流分析，共划分为616个子流域，对应每个子流域都有自己的主河道，见图12。

对通辽市平原区开展基于水循环模拟的大尺度地下水数值仿真应用，圈定了通辽市平原区所属的491个子流域。平原区网格单元以2km为间距进行剖分，共剖分120行139列，共计正方形网格单元16680个，其中有效网格单元10539个，并建立子流域与网格单元间的空间嵌套关系，见图13。

达到的应用效果如下：

1、体现出了本发明在大空间尺度流域/区域和长仿真期条件下应用的能力。在5.9万km2的通辽市域完成了从1998年至2007年共计10年（3652天）的仿真模拟应用；

2、系统运行效率很高。在CPU为Intel T9800（频率为2.93GHZ），4G内存的运行环境下，系统运行时间仅为10分20秒左右完成所有工作，表现出很高的运行效率；

3、充分展示了水循环模拟与地下水数值仿真过程的双向反馈及其精度优势。以下是数值仿真过程中取得的部分成果。表1为程序运行完后获得的每年的地下水动态评价定量数据，图14为2007年地下水数值仿真获得的水位等值线与通辽市当地实测数据的对比。在图14中，标记141表示通辽市当地实测水位等值线，标记142表示地下水数值仿真获得的水位等值线。

 

表1 通辽市1998~2007年平原区地下水动态评价分析(亿m ³ )

 

实施例二：

本实施例是实施例一的改进，是实施例一关于空间嵌套处理的细化。本实施例所述的“空间嵌套处理”包括以下子步骤：

1)     确定子流域与网格单元的空间联系。使用ARCGIS或其他通用地理信息系统工具对水循环模拟系统划分的多个子流域与地下水数值仿真系统划分的多个网格单元进行空间上的叠加，从而得出模拟区域中各子流域与各网格单元的从属关系。即对某具体的子流域，确定与其有空间重叠的网格单元的个数及位置并进行记录。所述的地理信息系统工具是用于分析、管理多种地理空间实体数据及其关系的一种软件工具，目前比较流行的ARCGIS、ARCVIEW、MAPGIS、MAPWINDOW等。（注：这些软件系统十分常见，网上搜搜到处都是相关介绍，不过都没有中文名称）

2)     确定各网格单元在各子流域中所占的面积比例。面积比例指某网格单元占某子流域的面积除以该子流域的面积，即若编号为

的网格单元与编号为的子流域在空间上有叠加关系，且两者之间的叠加面积为，而该子流域的面积为

，则该网格单元占该子流域的面积比例可表达为：

                          （2）

其中

为网格单元的编号；

为子流域的编号； 

为编号为

的某网格单元与编号为

的某子流域之间的叠加面积；

为编号为

的子流域的面积。若某网格单元与某子流域在空间上没有叠加关系，则

。

实施例三：

本实施例是实施例一的改进，是实施例一关于同步耦合的细化。如前所述，同步耦合指在同一个时间步长内，两个系统完成信息之间的交互反馈。时间步长是进行模拟或仿真时的时间片段，整个模拟期由多个时间步长组成，对每个时间步长两个系统均进行同步耦合，直至模拟结束。在一个时间步长内，本实施例所述的“同步耦合”包括以下子过程步骤：

1)   在一个时间步长内，水循环模拟系统先完成该时间步长的水循环模拟，并计算各子流域内地下水循环的各项通量，其中与地下水有关的垂向循环通量包括降水入渗补给量、河道/水库/湿地/渠系等地表水体渗漏补给量、井灌回归补给量、基流排泄量、潜水蒸发量、地下水开采量。即：

                                  （3）

                                  （4）

其中

为某子流域的地下水垂向补给通量总和(L³T^-1)；

为子流域的编号；

、

、

分别为该子流域的降水入渗量、河道/水库/湿地/渠系等地表水体渗漏量和灌溉渗漏量(L³T^-1)；为该子流域的地下水排泄通量总和(L³T^-1)；、

、

分别为该子流域的潜水蒸发量、基流排泄量和地下水开采量(L³T^-1)；

2)   水循环模拟系统将与地下水有关的垂向循环通量信息传递给地下水数值仿真系统，这通过子流域与网格单元的从属关系和面积比例关系进行。各地下水网格单元所获得的源汇量强度大小可由下式确定，即：

                        （5）

其中

为编号为

的网格单元获得的源汇量强度(L³T^-1)；为模拟空间内子流域的个数；

为子流域的编号；

为编号为

的网格单元占编号为

的子流域的面积比例(-)，其他符号意义同前；

3)   地下水数值仿真系统获得了该时间步长内所需的全部源汇量信息，即公式（1）中 W 项(代表来自源汇处的水量)已经可以确定，因此完成本时间步长的仿真过程并更新所有网格单元的地下水位。

4)   地下水数值仿真系统将地下水位信息传递给水循环模拟系统。在水循环模拟系统获得来自地下水数值仿真系统在该时段内的信息后，将通过面积加权法更新各子流域的平均地下水位和地下水埋深，该步骤也通过子流域-网格单元空间嵌套关系和面积比例关系进行。

                              （6）

                              （7）

其中为编号为

的子流域的平均地下水位(L)，

为编号为

的子流域的平均地下埋深(L)；

为编号为

的网格单元的地下水位；

为模拟空间内网格单元的总个数；其他符号意义同前。

5)   该时间步长内的同步耦合过程结束，系统进入下一时间步长，在下一个时间步长，水循环模拟系统将用上一个时间步长所更新的地下水位和地下水埋深作为基础进行水循环模拟。重复以上步骤直至所有的时间步长完成，从而结束整个模拟期。模拟期可长可短，根据应用要求而定，具体可从几天到数十年。

实施例四：

本实施例是实施例三的改进，是实施例三关于网格单元划分的细化。本实施例所述的网格单元划分方式是有限差分法、有限元法、有限体积法中的一种。

地下水数值仿真模拟的基本思想是把连续的模拟区域用有限个离散点构成的网格来代替，这些离散点称作网格的节点。通过把模拟区域上的连续变量的地下水三维动力学方程用在节点上定义的离散变量函数来近似，从而进行地下水的仿真模拟。在具体应用时，网格单元可有多种形式的划分。

有限差分法将三维含水层系统划分为一个三维的网格系统，整个含水层系统被剖分为若干层，每一层又剖分为若干行和若干列。按照节点位置选取的不同，有限差分法又可分为单元中心格式和节点中心格式两种，通过建立每个单元/节点与相邻单元/节点之间的水量平衡关系并进行积分来获得近似解。

有限元法通常将三维含水层系统划分为多个称为有限元的小的互连子域，对每一互连子域假定一个合适的(较简单的）近似解，然后推导求解模拟区域总的满足条件，从而得到问题的解。有限元法最常用的是三角形有限元（即线性有限元），为由多个三角形组成的网状系统，也有非线性的有限元，但实用的很少。

有限体积法将模拟区域划分为一系列不重复的控制体积，并使每个节点周围有一个控制体积；将待解的微分方程对每一个控制体积积分，便得出一组离散方程，其中的未知数是网格点上的因变量的数值。

    虽然网格单元划分方式有所不同，从而导致地下水三维动力学方程在数学离散处理上有所差异，但最终原理仍为流入和流出某个网格单元的水流之差应等于该单元中贮水量的变化，即地下水三维动力学方程可以简单地表示为以下水量平衡方程：

                                    （8）

式中：

：单位时间内流进或流出某网格单元的水量(L³T^-1)；：含水层的贮水率(L^-1)。它表示当水头变化为一个单位时，该网格单元单位体积中所吸收或释放的水量；

：网格单元的体积(L³)；

：某一时间段内水头的变化(L)；

：时间变化量(T)。

如果对模型所包含的网格单元逐个写出类似于公式(xx)形式的方程，则可得一个线性方程组。这个方程组可用矩阵的形式表示为：

                                    （9）

其中为水头的系数矩阵；

为所求的水头矩阵；而

表示各个方程中所包含的所有常数项和已知项。有时也称为右侧项。

以上矩阵形式的方程组的规模由网格单元的数量决定，有时甚至包含有上百万个未知数。这样通常需要采用迭代的方法进行求解，如强隐式法(Strongly Implicit Procedure，简称SIP)、分层逐次超松弛法(Slice Successive Overrelaxation，简称SSOR)或预调共轭梯度法(Preconditioned Conjugate-gradient，简称PCG)。

地下水数值仿真的目的是预测模拟区域的地下水水头变化，这取决于初始水头分布、边界条件、各种水文地质参数的分布以及各种外界源和汇的分布与强度。模拟时总是从初始水头开始，每一步求出每个时间步长结束时的水头值，并用该值作为下一时间步长的初始值，不断重复这样的过程，直至所要求的时间结束。

目前用有限差分法的地下水数值仿真系统有Visual MODFLOW、Processing Modflow等。Visual Modflow是由加拿大waterloo 水文地质公司开发研制的，1994年8月首次在国际上公开发行，是目前国际上最盛行、且被各国同行一致认可的地下水数值仿真软件。与其他软件平台相比，Processing Modflow则是一款轻量级的地下水模拟平台，由W.~H. Chiang和W.Kinzalbach开发，界面相对简单但是实用，不足之处在于后处理功能不够强大。

用有限元法最著名的地下水数值仿真系统是Feflow，由德国Wasy水资源规划系统研究所开发研制，由于其功能比较全面，最近其关注程度逐渐上升。

实施例五：

本实施例是实施例四的改进，是实施例四关于有限差分法的细化。本实施例所述的有限差分法的差分方式是单元中心格式和节点中心格式中的一种。

有限差分法的单元划分可以分为单元中心格式和节点中心格式。无论哪种格式，整个含水层均由两组正交的平行线剖分为一系列矩形方格。在单元中心法中，这些平行线构成了网格单元的边界。而在节点中心法中，节点位于这两组正交直线的交点上，网格单元的边界则位于两个相邻节点中线的位置上。

这两种划分方法均有优缺点，两者的区别主要在于模拟区域边界的处理上。在单元中心法中，模拟区域边界总是位于网格单元的边线上。而节点中心法中，模拟区域的边界则与网格单元的中心处在同一位置。单元中心法的优点在于处理流量边界时比较容易，可以将边界流量值直接加入到位于边界上的网格单元上，节点中心法的优点在于处理水头边界时比较容易，因为网格单元中心的位置与边界的位置重合，边界处的水头就是网格单元中心的水头。

实施例六：

本实施例是实施例一的改进，是实施例一中关于基于子流域空间离散模型的细化。水循环模拟系统是对自然界的水文现象进行抽象和概化，建立水文过程的数学结构与逻辑结构并编制程序进行近似计算表达的一种软件工具。所述的基于子流域空间离散技术的水循环模拟系统是指按照河道汇流格局，将模拟空间划分为多个相互分离但具有地表水力联系的子流域，从而进行分布式模拟的水循环模拟系统。

基于子流域空间离散技术的水循环模拟系统包括SWAT、新安江、HSPF等。

SWAT水循环模拟系统（Soil and Water Assessment Tool）是由美国农业部（USDA）的农业研究中心Jeff Amonld博士1994年开发的，采用FORTRAN语言编写。SWAT模型采用日为时间尺度连续计算。该水循环模拟系统不仅模拟水量，还附带对水质、以及杀虫剂的输移与转化过程的模拟。该水循环模拟系统在子流域内分各种水文响应单元以细化下垫面、土壤和人类水土管理，并将水文过程模拟分为子流域内陆面过程和子流域之间的河道汇流过程两部分，整个水分循环系统遵循水量平衡规律。SWAT水循环模拟系统具有较强的物理基础，考虑了较为丰富的水文过程，如采用Green&Ampt方法或SCS曲线法计算陆面产流过程；汇流方面采用分段连续演算的Muskingum法或河槽蓄量法；用田间持水度控制土壤水的分层下渗；在蒸散发方面细分为冠层截留蒸发、土壤蒸发、植被蒸腾、积雪升华等多个子过程并用Penman-Menteith等方法计算；在地下水方面能够计算基流量、降水入渗补给量、潜水蒸发量、河道/湿地/水库等地表水渗漏补给量等。因此具有在子流域层次上对大气水、土壤水、地表水、地下水之间的转化关系进行定量计算的能力，能够产生地下水数值仿真需求的详细数据。

新安江水循环模拟系统是1973年我国河海大学提出的，在我国南方新安江流域首先获得成功应用，后又逐步发展成为我国少有的一个在国际上具有影响力的水循环模拟系统，采用FORTRAN语言编写。该水循环模拟系统属于概念性的水循环模拟系统，以模拟产汇流过程为主。其将子流域内的下垫面分为透水区和不透水区分开模拟，按照三层蒸散发模式计算流域蒸散发，按蓄满产流概念计算降雨产生的总径流量，采用流域蓄水曲线考虑下垫面不均匀对产流面积变化的影响。在径流成分划分方面，对三水源情况，按“山坡水文学”产流理论用一个具有有限容积和测孔、底孔的自由水蓄水库把总径流划分成饱和地面径流、壤中水径流和地下水径流。在汇流计算方面，单元面积的地面径流汇流一般采用单位线法，壤中水径流和地下水径流的汇流则采用线性水库法。河网汇流一般采用分段连续演算的Muskingum法或滞时演算法。新安江水循环模拟系统的缺点在于对地下水分项循环过程的计算较少，虽然涉及大气水、土壤水、地表水、地下水之间的转化关系的模拟，但只笼统给出了地下水的总补给量和基流量，产生的地下水数值仿真需求的数据较粗略。

HSPF(Hydrological Simulation Program-FORTlRAN)是美国环境保护署(EPA)提出的一个水循环模拟系统，由Johanson 等人开发。与SWAT相类似，HSPF水循环模拟系统是一个综合的系统，不仅模拟流域的水文过程，还包括土壤流失、污染物传输，水温、泥沙传输、营养物和化学物相互反应等。在水文方面，该系统能对透水地面、不透水地面、河流和完全混合型湖泊水库三种不同性质的地表水文过程进行模拟。该系统自上而下分为树冠层、植被层、各土壤层(包括表层土壤、上土壤层、下土壤层、地下水涵养层)，考虑降雨(降雪)、植被截留、不同土壤层的水分传导、蓄积和蒸散发等十余个水文子过程，其中对降雨径流的描述较为复杂完善，较好地再现了流域的降雨径流过程。在地下水方面其分项过程不如SWAT考虑得那么细致，但比新安江模型略强，产生的地下水数值仿真需求的数据比较一般。

最后应说明的是，以上仅用以说明本发明的技术方案而非限制，尽管参照较佳布置方案对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案（比如各个步骤的顺序关系等）进行修改或者等同替换，而不脱离本发明技术方案的精神和范围。 

Claims (6)

1.一种基于水循环模拟的大尺度地下水数值仿真方法，其特征在于所述方法的步骤如下：

建立水循环模拟系统的步骤：用于采用基于子流域空间离散技术的水循环模拟系统，以确定区域范围内地下水的垂向循环通量的时空分布，从而产生地下水数值仿真需求的数据信息，包括降水入渗补给量、河道/水库/湿地/渠系地表水体渗漏补给量、井灌回归补给量、地下水基流量、潜水蒸发量、地下水开采量；

建立地下水数值仿真系统的步骤：用于采用基于矩形网格单元空间离散技术的地下水数值仿真系统，模拟地下水在研究区域内的侧向流动，其仿真原理基于以下三维地下水动力学方程：

                  （1）

其中

，和

为渗透系数在X, Y和Z方向上的分量；

为水头；

为单位体积流量，用以代表来自源汇处的水量；

为孔隙介质的贮水率；

为时间，使用网格单元划分法对被模拟区域进行空间划分后，对被模拟区域进行数学模拟；

进行水循环模拟系统与地下水数值仿真系统空间嵌套处理的步骤：用于在水循环模拟系统所划分的子流域与地下水数值仿真系统所划分的网格单元之间建立空间关联；

空间嵌套处理将为水循环模拟系统和地下水数值仿真系统在同步耦合过程中的双向数据信息交互提供基础；

水循环模拟系统与地下水数值仿真系统融合的步骤：用于以日尺度为时间步长，通过水循环模拟方法和地下水数值仿真系统在每个时间步长内的双向信息交互，构成两者之间的同步耦合；所述的“同步耦合”指在同一个时间步长内，两个系统完成信息之间的交互反馈；时间步长是进行模拟或仿真时的时间片段，整个模拟期由多个时间步长组成，对每个时间步长两个系统均进行同步耦合，直至模拟结束；所述的“双向信息交互”是以水循环模拟过程中使用的子流域和地下水数值仿真系统中使用的格单元之间的空间嵌套处理作为交互基础；双向信息交互内容包括两方面：一是水循环模拟系统得出的垂向循环通量时空分布信息传递给地下水数值仿真系统，提供地下水仿真所需的源汇处的水量，垂向循环通量时空分布信息包括降水入渗补给量、河道渗漏补给量、水库渗漏补给量、渠系输水过程渗漏补给量、井灌回归补给量、地下水基流量、潜水蒸发量、地下水开采量；二是地下水数值仿真系统得出的地下水位和地下水埋深的信息传递给水循环模拟系统，辅助水循环模拟系统计算地下水与地表水、土壤水之间的转化量。

2.根据权利要求1所述的方法，其特征在于，所述的“空间嵌套处理”包括以下子步骤：

空间嵌套处理过程为：

确定子流域与网格单元的空间联系，使用通用地理信息系统工具对水循环模拟系统划分的多个子流域与地下水数值仿真系统划分的多个网格单元进行空间上的叠加，从而得出模拟区域中各子流域与各网格单元的从属关系；

确定各网格单元在各子流域中所占的面积比例，即若编号为的网格单元与编号为的子流域在空间上有叠加关系，且两者之间的叠加面积为

，而子流域的面积为，则

网格单元占

子流域的面积比例为：

若网格单元与子流域在空间上没有叠加关系，则

。

3.根据权利要求1所述的方法，其特征在于，所述的“同步耦合”包括以下子过程步骤：

在一个时间步长内，水循环模拟系统先完成该时间步长的水循环模拟，并计算各子流域内地下水循环的各项通量，其中与地下水有关的垂向循环通量包括降水入渗补给量、河道/水库/湿地/渠系地表水体渗漏补给量、井灌回归补给量、基流排泄量、潜水蒸发量、地下水开采量；即：

                                  （3）

                                  （4）

其中为某子流域的地下水垂向补给通量总和(L³T^-1)；为子流域的编号；、

、

分别为该子流域的降水入渗量、河道/水库/湿地/渠系地表水体渗漏量和灌溉渗漏量(L³T^-1)；

为该子流域的地下水排泄通量总和(L³T^-1)；

、

、

分别为该子流域的潜水蒸发量、基流排泄量和地下水开采量(L³T^-1)；

水循环模拟系统将与地下水有关的垂向循环通量信息传递给地下水数值仿真系统，这通过子流域与网格单元的从属关系和面积比例关系进行，各地下水网格单元所获得的源汇量强度大小可由下式确定，即：

                        （6）

其中为编号为的网格单元获得的源汇量强度(L³T^-1)；

为模拟空间内子流域的个数；

为子流域的编号；为编号为

的网格单元占编号为

的子流域的面积比例(-)，其他符号意义同前；

地下水数值仿真系统获得了该时间步长内所需的全部源汇量信息，即公式（1）中 W 项(代表来自源汇处的水量)已经可以确定，因此完成本时间步长的仿真过程并更新所有网格单元的地下水位；

地下水数值仿真系统将地下水位信息传递给水循环模拟系统；

在水循环模拟系统获得来自地下水数值仿真系统在该时段内的信息后，将通过面积加权法更新各子流域的平均地下水位和地下水埋深，本步骤也通过子流域-网格单元空间嵌套关系和面积比例关系进行，即：

                              （2）

                              （2）

其中

为编号为

的子流域的平均地下水位(L)，

为编号为的子流域的平均地下埋深(L)；

为编号为

的网格单元的地下水位；

为模拟空间内网格单元的总个数；其他符号意义同前；

本时间步长内的同步耦合过程结束，进入下一时间步长；

在下一个时间步长，水循环模拟系统将用上一个时间步长所更新的地下水位和地下水埋深作为基础进行水循环模拟；

重复以上步骤直至所有的时间步长完成，从而结束整个模拟期。

4.根据权利要求3所述的方法，其特征在于，所述的网格单元划分方式是有限差分法、有限元法、有限体积法中的一种。

5.根据权利要求4所述的方法，其特征在于，所述的有限差分法的差分方式是单元中心格式和节点中心格式中的一种。

6.根据权利要求1所述的方法，其特征在于，所述的基于子流域空间离散技术的水循环模拟系统是SWAT水循环模拟系统、新安江水循环模拟系统、HSPF水循环模拟系统中的一种。

Images