本发明涉及地下水环境，尤其涉及一种污染物淋溶通量的数值计算方法及装置。
背景技术：
1、在地表污染源与水流入渗的双重驱动作用下，包气带中污染物的淋溶过程直接决定了污染源是否对地下水环境产生影响，而地下水位的波动又影响污染物在包气带中的迁移转化过程。因此，在计算分析包气带中污染物的淋溶通量时，需重点考虑包气带与饱水带间的交互作用。同时，受地表污染源类型、降雨类型、包气带岩性、地下水位波动等因素影响，包气带中污染物的淋溶通量并不是一个稳定值，而是随时间呈现动态波动。
2、关于污染物由地表进入地下水的淋溶过程，现有研究多集中于地表浅层或根系区，对深层包气带与饱水带间的相互作用研究较少；或直接采用包气带底部的最大或稳定的污染物浓度作为饱水带的补给项，未考虑包气带与饱水带间的实时交互作用，存在一定的理论偏差。因此，建立地下环境中包气带-饱水带间互作用耦合模型，以两带交界面间污染物通量作为包气带淋溶通量更为合理。
技术实现思路
1、为解决上述技术问题，本发明提出了一种污染物淋溶通量的数值计算方法及装置。
2、根据本发明的第一方面，提供一种污染物淋溶通量的数值计算方法，所述方法包括以下步骤：
3、步骤s1：构建水文地质模型，由所述水文地质模型向包气带饱水带互作用耦合模型提供信息；
4、所述包气带饱水带互作用耦合模型包括包气带模型及位于包气带模型下方的饱水带模型；将所述包气带模型划分为多个包气带子模型，对每个包气带子模型，所述饱水带模型均具有对应的饱水带模型区域与所述包气带子模型对应；每个包气带子模型及与所述包气带子模型对应的饱水带模型区域组成一个包气带饱水带对；
5、确定所述包气带模型及所述饱水带模型的初始条件、边界条件、源汇项、水文地质参数及污染物运移参数；
6、步骤s2：设置所述包气带饱水带互作用耦合模型进行仿真的应力期数num1、每个应力期内的流动时间步数num2；初始化当前仿真的应力期数knum1为0、当前流动时间步数knum2为0；
7、步骤s3：对每个包气带饱水带对，均并行执行以下操作：
8、步骤s31：设置用于对所述饱水带模型区域中的水流进行仿真的饱水带模型水流模块、以及所述饱水带模型水流模块进行仿真的第一时间步长step1；设置用于对所述饱水带模型区域中的溶质运移进行仿真的饱水带模型溶质模块、以及所述饱水带模型溶质模块进行仿真的第二时间步长step2；设置用于对所述包气带子模型进行仿真的包气带模块，所述包气带模块进行仿真的第三时间步长step3；其中，step1>step2>step3；时间步长是进行一次完整的仿真的时长；
9、步骤s32：若knum1等于num1,则所述包气带饱水带对处理完毕，获取此时所述饱水带模型溶质模块仿真得到的饱水带模型区域的地下水污染物浓度，所述饱水带模型水流模块仿真得到的饱水带模型区域的地下水水位，包气带模块仿真得到的包气带底部的水流通量和污染物通量，作为所述包气带饱水带对的仿真结果；否则，进入步骤s33，将knum2赋值为0；
10、步骤s33：若knum2等于num2，则所述包气带饱水带对的当前应力期处理完毕，将knum1赋值为knum1加1，进入步骤s32；否则，进入步骤s34；
11、步骤s34：所述饱水带模型水流模块将所述饱水带模型的地下水水位传递给所述包气带模块，所述饱水带模型溶质模块将所述饱水带模型的地下水污染物浓度传递给所述包气带模块，所述饱水带模型的地下水水位和所述饱水带模型的地下水污染物浓度作为所述包气带子模型的下边界条件；
12、其中，所述饱水带模型的地下水水位及所述饱水带模型的地下水污染物浓度均是基于全部饱水带模型水流模块、全部饱水带模型溶质模块的实时仿真结果得到的；
13、步骤s35：所述包气带模块、所述饱水带模型溶质模块、所述饱水带模型水流模块进行并行仿真；且在所述第一时间步长step1内，所述饱水带模型溶质模块及所述包气带模块均进行多次的仿真；在所述第一时间步长step1到达时，获取所述饱水带模型的地下水水位、所述饱水带模型的地下水污染物浓度、所述饱水带模型溶质模块仿真得到饱水带模型区域的污染物浓度、以及所述饱水带模型水流模块得到饱水带模型区域的地下水水位；保存所述饱水带模型溶质模块仿真得到的饱水带模型区域的污染物浓度及所述饱水带模型水流模块仿真得到饱水带模型区域的地下水水位；将knum2赋值为knum2加1，进入步骤s33；
14、其中，所述饱水带模型溶质模块基于包气带子模型底部的污染物通量进行仿真，所述包气带子模型底部的污染物通量是所述包气带模块进行仿真时，每达到第三时间步长时仿真得到的；所述饱水带模型水流模块基于包气带子模型底部的水流通量进行仿真，所述包气带子模型底部的水流通量是所述包气带模块进行仿真时，每达到第三时间步长时仿真得到的；
15、步骤s4：基于每个包气带饱水带对的仿真结果，确定污染物淋溶通量。
16、优选地，所述饱水带模型的地下水水位、地下水污染物浓度是基于全部饱水带模型区域各自的地下水水位、地下水污染物浓度求均值得到的。
17、优选地，在所述第一时间步长step1到达时，获取所述饱水带模型的地下水水位、所述饱水带模型的地下水污染物浓度、所述饱水带模型溶质模块仿真得到饱水带模型区域的污染物浓度、以及所述饱水带模型水流模块得到饱水带模型区域的地下水水位，其中：
18、所述包气带模块进行了多次仿真，每次仿真的时长为所述第三时间步长step3，每次仿真均对水流及污染物第一运移过程公式进行多次迭代求解，以生成包气带子模型底部的污染物通量、以及包气带子模型底部的水流通量；每次仿真结束时，所述包气带模块均将本次仿真得到的包气带子模型底部的污染物通量、以及包气带子模型底部的水流通量分别发送给所述饱水带模型溶质模块、所述饱水带模型水流模块；所述饱水带模型水流模块进行了一次仿真，基于包气带子模型底部的水流通量，对水流第二运移过程公式进行多次迭代求解生成饱水带模型区域的地下水水位以及与地下水水位相关的信息；所述溶质模块进行了多次仿真，每次仿真均基于地下水流信息以及包气带子模型底部的污染物通量，对污染物第二运移过程公式进行多次迭代求解，以生成饱水带模型区域的地下水污染物浓度；其中，所述地下水流信息包括水头、饱和带厚度、单元截面各方向水流通量及各源/汇位置和流速。
19、优选地，所述水流及污染物第一运移过程公式为：
20、
21、式中，θ为体积含水量，t为时间，z为空间直角坐标系上竖轴方向上的坐标值，ks为饱和渗透系数，kr为相对渗透系数，h为压力水头，ss为水流源汇项；c、s、g分别为包气带液相、吸附相及气相污染物浓度，ρ为土壤容重，av为体积含气量，dw为液相水动力弥散系数，dg为气相扩散系数，q为达西流速，kw、ks、kg分别为液相、吸附相、气相污染物反应速率常数，sp为污染物源汇项。
22、优选地，所述水流第二运移过程公式为：
23、
24、所述污染物第二运移过程公式为：
25、
26、其中，h为地下水水位；kx、ky、kz分别为空间直角坐标系中横轴、纵轴、竖轴方向上的渗透系数；w为地下水流源汇项，ss为释水系数；r为迟滞因子，ρb为含水层介质密度，i、j分别代表空间直角坐标系中横轴、纵轴的方向，xi、xj分别对应于空间直角坐标系中横轴、纵轴方向的数值，dij为水动力弥散系数张量，λ1为液相反应速率常数，λ2为吸附相反应速率常数，cs为源汇项中污染物浓度。
27、根据本发明第二方面，提供一种污染物淋溶通量的数值计算装置，所述装置包括：
28、包气带饱水带互作用耦合模型构建模块：配置为构建水文地质模型，由所述水文地质模型向包气带饱水带互作用耦合模型提供信息；
29、所述包气带饱水带互作用耦合模型包括包气带模型及位于包气带模型下方的饱水带模型；将所述包气带模型划分为多个包气带子模型，对每个包气带子模型，所述饱水带模型均具有对应的饱水带模型区域与所述包气带子模型对应；每个包气带子模型及与所述包气带子模型对应的饱水带模型区域组成一个包气带饱水带对；
30、确定所述包气带模型及所述饱水带模型的初始条件、边界条件、源汇项、水文地质参数及污染物运移参数；
31、初始化模块：配置为设置所述包气带饱水带互作用耦合模型进行仿真的应力期数num1、每个应力期内的流动时间步数num2；初始化当前仿真的应力期数knum1为0、当前流动时间步数knum2为0；
32、迭代仿真模块：配置为对每个包气带饱水带对，均并行执行以下操作：
33、第一仿真子模块：配置为设置用于对所述饱水带模型区域中的水流进行仿真的饱水带模型水流模块、以及所述饱水带模型水流模块进行仿真的第一时间步长step1；设置用于对所述饱水带模型区域中的溶质运移进行仿真的饱水带模型溶质模块、以及所述饱水带模型溶质模块进行仿真的第二时间步长step2；设置用于对所述包气带子模型进行仿真的包气带模块，所述包气带模块进行仿真的第三时间步长step3；其中，step1>step2>step3；时间步长是进行一次完整的仿真的时长；
34、第二仿真子模块：配置为若knum1等于num1,则所述包气带饱水带对处理完毕，获取此时所述饱水带模型溶质模块仿真得到的饱水带模型区域的地下水污染物浓度，所述饱水带模型水流模块仿真得到的饱水带模型区域的地下水水位，包气带模块仿真得到的包气带底部的水流通量和污染物通量，作为所述包气带饱水带对的仿真结果；否则，触发第三仿真子模块，将knum2赋值为0；
35、第三仿真子模块：配置为若knum2等于num2，则所述包气带饱水带对的当前应力期处理完毕，将knum1赋值为knum1加1，触发第二仿真子模块；否则，触发第四仿真子模块；
36、第四仿真子模块：配置为所述饱水带模型水流模块将所述饱水带模型的地下水水位传递给所述包气带模块，所述饱水带模型溶质模块将所述饱水带模型的地下水污染物浓度传递给所述包气带模块，所述饱水带模型的地下水水位和所述饱水带模型的地下水污染物浓度作为所述包气带子模型的下边界条件；
37、其中，所述饱水带模型的地下水水位及所述饱水带模型的地下水污染物浓度均是基于全部饱水带模型水流模块、全部饱水带模型溶质模块的实时仿真结果得到的；
38、第五仿真子模块：配置为所述包气带模块、所述饱水带模型溶质模块、所述饱水带模型水流模块进行并行仿真；且在所述第一时间步长step1内，所述饱水带模型溶质模块及所述包气带模块均进行多次的仿真；在所述第一时间步长step1到达时，获取所述饱水带模型的地下水水位、所述饱水带模型的地下水污染物浓度、所述饱水带模型溶质模块仿真得到饱水带模型区域的污染物浓度、以及所述饱水带模型水流模块得到饱水带模型区域的地下水水位；保存所述饱水带模型溶质模块仿真得到的饱水带模型区域的污染物浓度及所述饱水带模型水流模块仿真得到饱水带模型区域的地下水水位；将knum2赋值为knum2加1，触发第三仿真子模块；
39、其中，所述饱水带模型溶质模块基于包气带子模型底部的污染物通量进行仿真，所述包气带子模型底部的污染物通量是所述包气带模块进行仿真时，每达到第三时间步长时仿真得到的；所述饱水带模型水流模块基于包气带子模型底部的水流通量进行仿真，所述包气带子模型底部的水流通量是所述包气带模块进行仿真时，每达到第三时间步长时仿真得到的；
40、计算模块：配置为基于每个包气带饱水带对的仿真结果，确定污染物淋溶通量。
41、根据本发明第三方面，提供一种电子设备，包括：
42、处理器，用于执行多条指令；
43、存储器，用于存储多条指令；
44、其中，所述多条指令，用于由所述存储器存储，并由所述处理器加载并执行如前所述方法。
45、根据本发明第四方面，提供一种计算机可读存储介质，所述存储介质中存储有多条指令；所述多条指令，用于由处理器加载并执行如前所述的方法。
46、本发明提出了一种污染物淋溶通量的数值计算方法及装置，基于地下环境中包气带与饱水带之间的实时交互作用，结合包气带底部水流通量与污染物通量、饱水带中地下水位与污染物浓度4个变量的双向传递，构建包气带-饱水带交互作用的污染物迁移转化耦合模型，解析包气带底部污染物淋溶通量的时间变化规律。
47、本发明具有以下技术效果：
48、(1)通过构建包气带-饱水带互作用耦合模型，获取了包气带污染物淋溶通量的动态变化特征，适用于地表污染源分布不均、地下环境复杂多变的场景，避免了高成本、长工期的野外监测工作；
49、(2)提出的包气带污染物淋溶通量的计算方法，符合实际的物理过程，可为评估区域地下水污染风险及保护地下水环境安全提供技术支撑。
50、上述说明仅是本发明技术方案的概述，为了能够更清楚了解本发明的技术手段，并可依照说明书的内容予以实施，以下以本发明的较佳实施例并配合附图详细说明如后。
1.一种污染物淋溶通量的数值计算方法，其特征在于，所述方法包括以下步骤：
2.如权利要求1所述的方法，其特征在于，所述饱水带模型的地下水水位、地下水污染物浓度是基于全部饱水带模型区域各自的地下水水位、地下水污染物浓度求均值得到的。
3.如权利要求2所述的方法，其特征在于，在所述第一时间步长step1到达时，获取所述饱水带模型的地下水水位、所述饱水带模型的地下水污染物浓度、所述饱水带模型溶质模块仿真得到饱水带模型区域的污染物浓度、以及所述饱水带模型水流模块得到饱水带模型区域的地下水水位，其中：
4.如权利要求3所述的方法，其特征在于，所述水流及污染物第一运移过程公式为：
5.如权利要求4所述的方法，其特征在于，所述水流第二运移过程公式为：
6.一种污染物淋溶通量的数值计算装置，其特征在于，所述装置包括：
7.一种电子设备，包括：
8.一种计算机可读存储介质，所述存储介质中存储有多条指令；所述多条指令，用于由处理器加载并执行如权利要求1-5中任一项所述方法。
