Newton迭代法求变密度2D-Richards方程的数值解及其在MARUN程序中的实现

潮间带浅层含水层地下水流和溶质运移过程受众多非线性因素影响，比如海潮周期性浸没造成的土壤孔隙水饱和度变化、低潮时潮间带地表可能产生的渗出面、沉积物水力性质的非均质性等等，导致其地下水流过程极其复杂；继而使得以变密度2D-Richards方程作为基本控制方程的海滩含水层地下水渗流模型的非线性急剧增加，其数值求解速度一直难以满足用户的实际需求。Newton迭代算法是一种高效稳定的线性化方法，通过使用函数的导数来构建一系列线性方程组去逼近原来非线性方程的解。目前它已经在一系列刻画土壤水渗流过程的Richards方程的求解中得到成功运用，可以合理应对不同土壤水力模型导致的水流方程的强非线性。但是这些问题都没有考虑海岸带的众多非线性因素对Newton迭代计算带来的影响。

本文以Newton迭代算法求解变密度2D-Richards方程为目标，具体开展了如下几方面的工作：详细阐述了变密度Richards方程的有限元数值离散过程；推导了Newton、Picard、Newton-Picard (NP) 三种求解方法的具体迭代格式，并用FORTRAN编程语言在海岸带地下水流数值模拟软件MARUN中实现了Newton、NP方法的迭代求解；通过两个传统案例验证了Newton迭代算法的可靠性，进一步比较了三种方法在求解不同海滩地下水流模型（包括三个海水入侵理论模型和一个野外场地实测模型）的计算效果。

计算结果表明，盐分带来的密度效应虽然显著增加了模型的计算时间，但是影响不同解法计算速度的主要因素仍然是变饱和度。模拟发现，如果海水水位波动的频率高振幅大，含水层的坡度大，网格剖分粗，VG模型参数α大n小，海浪作用剧烈，Newton迭代算法相对于Picard迭代算法就有着更高的收敛精度和更快的计算速度，但是会比Picard迭代算法花费更多的CPU和计算机存储。此外，NP求解方法相对于Newton求解方法而言，克服了求解过程对初值估计的敏感性；相对于Picard求解方法而言，则优化了解的收敛性，有更快的计算速度和更好的收敛行为。

The simulation of groundwater flow and solute transport processes in the shallow aquifer of the intertidal zone is extremely complex due to numerous nonlinear factors, such as the variable saturation caused by the periodic inundation of seawater, the time-dependent seepage faces occurring on the beach surface of the intertidal zone, and the heterogeneity of hydraulic properties. These nonlinear factors result in a rapid increase in computational cost in simulating the coastal groundwater model based on the variable-density 2D-Richards equation, making it difficult to meet the actual needs of users in terms of computation speed. The Newton iteration algorithm is an efficient and stable linearization method that constructs a series of linear equation sets to approximate the solution of the original nonlinear equation by using the derivatives of the nonlinear discretization equations. Currently, the Newton iteration method has been successfully applied to solving a series of Richards equations that describe soil water seepage processes and can reasonably deal with the strong nonlinearity of groundwater flow models. However, the impacts of numerous nonlinear factors in the coastal zone on the Newton iteration calculations are not considered.
This paper mainly studies the computational performance of the Newton iteration method in solving the variable-density 2D Richards equation, including (1) a detailed explanation of the finite element numerical discretization process of the variable-density Richards equation; (2) derivation of the specific iteration scheme of the Newton, Picard, and Newton-Picard (NP) solutions, (3) implementation of the iterative solution of the Newton and NP methods in the MARUN numerical simulation software for groundwater flow in the coastal zone using the FORTRAN programming language; and (4) validation of the reliability of the Newton iteration algorithm through two traditional examples, and further analyses and discussions of its computational performance in solving different groundwater flow models in coastal aquifers (including three seawater intrusion theoretical models and one case-study model).
The results show that although the density effect caused by salinity significantly increases the calculation time, the main reason for the difference in calculation speed between different iteration methods is still the variable saturation. Simulation results show that higher frequency, larger amplitude of seawater level fluctuations, larger slope of the beach surface, coarser grid division, larger VG model parameter α and smaller VG model parameter n, the Newton iteration has higher convergence accuracy and faster calculation speed than the Picard iteration, but requires more CPU and computer storage than the Picard format. On the other hand, compared with the Newton solution, the NP solution overcomes the sensitivity of the solution process to initial value estimation, and compared with the Picard solution, the NP solution optimizes the convergence of the solution, with faster calculation speed and better convergence behavior.

[1] Richards L A. Capillary conduction of liquids through porous mediums[J]. Physics, 1931, 1(5): 318-333.
[2] Yeh T C J, Xiang J, Suribhatla R M, et al. River stage tomography: A new approach for characterizing groundwater basins[J]. Water Resources Research, 2009, 45(5): W05409.
[3] Boughton W. A review of the USDA SCS curve number method[J]. Soil Research, 1989, 27(3): 511-523.
[4] Harbaugh A W, Banta E R, Hill M C, et al. Modflow-2000, the u. s. geological survey modular ground-water model-user guide to modularization concepts and the ground-water flow process[J]. Open-file report U S Geological Survey, 2000(92): 134.
[5] Geng X, Heiss J W, Michael H A, et al. Geochemical fluxes in sandy beach aquifers: Modulation due to major physical stressors, geologic heterogeneity, and nearshore morphology[J]. Earth-Science Reviews, 2021, 221: 103800.
[6] Li H, Boufadel M C. Long-term persistence of oil from the Exxon Valdez spill in two-layer beaches[J]. Nature Geoscience, 2010, 3(2): 96-99.
[7] Guo Q, Li H, Boufadel M C, et al. Hydrodynamics in a gravel beach and its impact on the Exxon Valdez oil[J]. Journal of Geophysical Research: Oceans, 2010, 115(C12): C12077.
[8] Xia Y, Li H, Boufadel M C, et al. Hydrodynamic factors affecting the persistence of the Exxon Valdez oil in a shallow bedrock beach[J]. Water Resources Research, 2010, 46(10): W10528.
[9] Brutsaert W F. A functional iteration technique for solving the Richards equation applied to two‐dimensional infiltration problems[J]. Water Resources Research, 1971, 7(6): 1583-1596.
[10] Philip J. The theory of infiltration: 1. The infiltration equation and its solution[J]. Soil science, 1957, 83(5): 345-358.
[11] Warrick A W, Lomen D O, Islas A. An analytical solution to Richards' equation for a draining soil profile[J]. Water Resources Research, 1990, 26(2): 253-258.
[12] Freeze R A, Witherspoon P A. Theoretical analysis of regional groundwater flow: 1. Analytical and numerical solutions to the mathematical model[J]. Water Resources Research, 1966, 2(4): 641-656.
[13] Neuman S P. Saturated-unsaturated seepage by finite elements[J]. Journal of the hydraulics division, 1973, 99(12): 2233-2250.
[14] Rushton K. Discrete time steps in digital computer analysis of aquifers containing pumped wells[J]. Journal of Hydrology, 1973, 18(1): 1-19.
[15] Moridis G J, Pruess K. T2SOLV: An enhanced package of solvers for the TOUGH2 family of reservoir simulation codes[J]. Geothermics, 1998, 27(4): 415-444.
[16] Cooley R L. Some new procedures for numerical solution of variably saturated flow problems[J]. Water Resources Research, 1983, 19(5): 1271-1285.
[17] Eymard R, Gutnic M, Hilhorst D. The finite volume method for Richards equation[J]. Computational Geosciences, 1999, 3(3): 259-294.
[18] Pour M A, Shoshtari M M, Adib A. Numerical solution of Richards equation by using of finite volume method[J]. World Applied Sciences Journal, 2011, 14(12): 1838-1842.
[19] List F, Radu F A. A study on iterative methods for solving Richards’ equation[J]. Computational Geosciences, 2016, 20(2): 341-353.
[20] Van Genuchten M T. A comparison of numerical solutions of the one-dimensional unsaturated—saturated flow and mass transport equations[J]. Advances in Water Resources, 1982, 5(1): 47-55.
[21] Celia M A, Bouloutas E T, Zarba R L. A general mass‐conservative numerical solution for the unsaturated flow equation[J]. Water Resources Research, 1990, 26(7): 1483-1496.
[22] Celia M A, Ahuja L R, Pinder G F. Orthogonal collocation and alternating-direction procedures for unsaturated flow problems[J]. Advances in Water Resources, 1987, 10(4): 178-187.
[23] Hao X, Zhang R, Kravchenko A. A mass-conservative switching method for simulating saturated–unsaturated flow[J]. Journal of Hydrology, 2005, 311(1-4): 254-265.
[24] Paniconi C, Aldama A A, Wood E F. Numerical evaluation of iterative and noniterative methods for the solution of the nonlinear Richards equation[J]. Water Resources Research, 1991, 27(6): 1147-1163.
[25] Paniconi C, Putti M. A comparison of Picard and Newton iteration in the numerical solution of multidimensional variably saturated flow problems[J]. Water Resources Research, 1994, 30(12): 3357-3374.
[26] Miller C T, Abhishek C, Farthing M W. A spatially and temporally adaptive solution of Richards’ equation[J]. Advances in Water Resources, 2006, 29(4): 525-545.
[27] Herbst M, Gottschalk S, Reißel M, et al. On preconditioning for a parallel solution of the Richards equation[J]. Computers & geosciences, 2008, 34(12): 1958-1963.
[28] Juncu G, Nicola A, Popa C, et al. Nonlinear multigrid methods for solving Richards' equation in two space dimensions[J]. Carpathian Journal of Mathematics, 2009: 82-91.
[29] 陈启生, 戚隆溪. 有植被覆盖条件下土壤水盐运动规律研究[J]. 水利学报, 1996(01): 38-46.
[30] 陈家军, 彭胜, 王金生等. 非饱和带水气二相流动参数确定实验研究[J]. 水科学进展, 2001(04): 467-472.
[31] 张华, 陈善雄, 陈守义. 非饱和土入渗的数值模拟[J]. 岩土力学, 2003(05): 715-718.
[32] 李毅, 王全九, 王文焰等. 入渗、再分布和蒸发条件下一维土壤水运动的数值模拟[J]. 灌溉排水学报, 2007(01): 5-8.
[33] 吴梦喜. 饱和-非饱和土中渗流Richards方程有限元算法[J]. 水利学报, 2009, 40(10): 1274-1279.
[34] 陈曦, 于玉贞, 程勇刚. 非饱和渗流Richards方程数值求解的欠松弛方法[J]. 岩土力学, 2012, 33(S1): 237-243.
[35] Zha Y, Shi L, Ye M, et al. A generalized Ross method for two-and three-dimensional variably saturated flow[J]. Advances in Water Resources, 2013, 54: 67-77.
[36] 黄文竹. 应用间断有限元方法模拟一维非饱和土壤水流问题[D]. 北京: 中国农业大学, 2014.
[37] 曾季才, 查元源, 杨金忠. 基于Richards方程切换的土壤水流及溶质运移数值模拟[J]. 水利学报, 2018, 49(07): 840-848.
[38] 陈远强, 杨永涛, 郑宏等. 饱和–非饱和渗流的数值流形法研究与应用[J]. 岩土工程学报, 2019, 41(02): 338-347.
[39] Homeier H H. A modified Newton method with cubic convergence: the multivariate case[J]. Journal of Computational and Applied Mathematics, 2004, 169(1): 161-169.
[40] Darvishi M, Barati A. A fourth-order method from quadrature formulae to solve systems of nonlinear equations[J]. Applied Mathematics and Computation, 2007, 188(1): 257-261.
[41] Cordero A, Martínez E, Torregrosa J R. Iterative methods of order four and five for systems of nonlinear equations[J]. Journal of Computational and Applied Mathematics, 2009, 231(2): 541-551.
[42] Chun C, Lee M Y. A new optimal eighth-order family of iterative methods for the solution of nonlinear equations[J]. Applied Mathematics and Computation, 2013, 223: 506-519.
[43] Petković M, Ilić S, Džunić J. Derivative free two-point methods with and without memory for solving nonlinear equations[J]. Applied Mathematics and Computation, 2010, 217(5): 1887-1895.
[44] Salimi M, Long N N, Sharifi S, et al. A multi-point iterative method for solving nonlinear equations with optimal order of convergence[J]. Japan Journal of Industrial and Applied Mathematics, 2018, 35(2): 497-509.
[45] Hills R, Porro I, Hudson D, et al. Modeling one‐dimensional infiltration into very dry soils: 1. Model development and evaluation[J]. Water Resources Research, 1989, 25(6): 1259-1269.
[46] Matthews C, Braddock R, Sander G. Modeling flow through a one-dimensional multi-layered soil profile using the Method of Lines[J]. Environmental Modeling & Assessment, 2004, 9(2): 103-113.
[47] Kirkland M R, Hills R, Wierenga P. Algorithms for solving Richards' equation for variably saturated soils[J]. Water Resources Research, 1992, 28(8): 2049-2058.
[48] Forsyth P A, Wu Y, Pruess K. Robust numerical methods for saturated-unsaturated flow with dry initial conditions in heterogeneous media[J]. Advances in Water Resources, 1995, 18(1): 25-38.
[49] Forsyth P, Wu Y, Pruess K. Robust numerical methods for saturated-unsaturated flow with dry initial conditions in heterogeneous media[J]. International Journal of Multiphase Flow, 1996, 22(S1): 92-92.
[50] Lehmann F, Ackerer P. Comparison of iterative methods for improved solutions of the fluid flow equation in partially saturated porous media[J]. Transport in porous media, 1998, 31(3): 275-292.
[51] Huang K, Mohanty B, Van Genuchten M T. A new convergence criterion for the modified Picard iteration method to solve the variably saturated flow equation[J]. Journal of Hydrology, 1996, 178(1-4): 69-91.
[52] Crevoisier D, Chanzy A, Voltz M. Evaluation of the Ross fast solution of Richards’ equation in unfavourable conditions for standard finite element methods[J]. Advances in Water Resources, 2009, 32(6): 936-947.
[53] Simunek J, Van Genuchten M T, Sejna M. The HYDRUS-1D software package for simulating the one-dimensional movement of water, heat, and multiple solutes in variably-saturated media[J]. University of California-Riverside Research Reports, 2005, 3: 1-240.
[54] Simunek J, van Genuchten M T, Sejna M. Development and applications of the HYDRUS and STANMOD software packages and related codes[J]. Vadose Zone Journal, 2008, 7(2): 587-600.
[55] Sinha B P, Mukherjee A. Parallel Sorting Algorithm Using Multiway Merge and Its Implementation on a Multi-Mesh Network[J]. Journal of Parallel and Distributed Computing, 2000, 60(7): 891-907.
[56] Li H, Farthing M, Dawson C, et al. Local discontinuous Galerkin approximations to Richards’ equation[J]. Advances in Water Resources, 2007, 30(3): 555-575.
[57] Zheng C, C. Hill M, Cao G, et al. MT3DMS: Model Use, Calibration, and Validation[J]. Transactions of the ASABE, 2012, 55(4): 1549-1559.
[58] Appelo C A J, Rolle M. PHT3D: A Reactive Multicomponent Transport Model for Saturated Porous Media[J]. Groundwater, 2010, 48(5): 627-632.
[59] Langevin C D. SEAWAT: A computer program for simulation of variable-density groundwater flow and multi-species solute and heat transport[R]. Reston: U.S. Geological Survey, 2009.
[60] Xiao-bin Z. The software (GMS) of Groundwater Modeling System[J]. Hydrogeology and Engineering Geology, 2003, 5: 53-55.
[61] Trefry M G, Muffels C. FEFLOW: A finite-element ground water flow and transport modeling tool[J]. Groundwater, 2007, 45(5): 525-528.
[62] Provost A M, Voss C I. SUTRA, a model for saturated-unsaturated, variable-density groundwater flow with solute or energy transport—Documentation of generalized boundary conditions, a modified implementation of specified pressures and concentrations or temperatures, and the lake capability[R]. Reston: U.S. Geological Survey, 2019.
[63] Croucher A E, O’Sullivan M J. Application of the computer code TOUGH2 to the simulation of supercritical conditions in geothermal systems[J]. Geothermics, 2008, 37(6): 622-634.
[64] Wei Y, Cao X. A COMSOL-PHREEQC Coupled Python Framework for Reactive Transport Modeling in Soil and Groundwater[J]. Groundwater, 2022, 60(2): 284-294.
[65] Van Genuchten M T. A closed‐form equation for predicting the hydraulic conductivity of unsaturated soils[J]. Soil science society of America journal, 1980, 44(5): 892-898.
[66] Park E-J. Mixed finite element methods for nonlinear second-order elliptic problems[J]. SIAM Journal on Numerical Analysis, 1995, 32(3): 865-885.
[67] Bergamaschi L, Putti M. Mixed finite elements and Newton‐type linearizations for the solution of Richards' equation[J]. International journal for numerical methods in engineering, 1999, 45(8): 1025-1046.
[68] Radu F, Pop I, Knabner P. On the convergence of the Newton method for the mixed finite element discretization of a class of degenerate parabolic equation[J]. Numerical Mathematics and Advanced Applications, 2006, 42: 1194-1200.
[69] Putti M, Paniconi C. Picard and Newton linearization for the coupled model for saltwater intrusion in aquifers[J]. Advances in Water Resources, 1995, 18(3): 159-170.
[70] Gear C W. Numerical initial value problems in ordinary differential equations[J]. Prentice-Hall series in automatic computation, 1971, 1: 1.
[71] Galucio A C, Deü J F, Mengué S, et al. An adaptation of the Gear scheme for fractional derivatives[J]. Computer Methods in Applied Mechanics and Engineering, 2006, 195(44): 6073-6085.
[72] Henry H R. Effects of dispersion on salt encroachment in coastal aquifers, in" Seawater in Coastal Aquifers"[J]. US Geological Survey, Water Supply Paper, 1964, 1613: C70-C80.
[73] Robinson C, Gibbes B, Li L. Driving mechanisms for groundwater flow and salt transport in a subterranean estuary[J]. Geophysical Research Letters, 2006, 33(3): L03402.
[74] Xiao K, Li H, Wilson A M, et al. Tidal groundwater flow and its ecological effects in a brackish marsh at the mouth of a large sub-tropical river[J]. Journal of Hydrology, 2017, 555: 198-212.
[75] Milly P. A mass-conservative procedure for time-stepping in models of unsaturated flow[M]. Burlington: Springer, 1984.
[76] Maurer D, Wieners C. A parallel block LU decomposition method for distributed finite element matrices[J]. Parallel Computing, 2011, 37(12): 742-758.
[77] Istok J. Groundwater Modeling by the Finite Element Method[M]. Florida: American Geophysical Union, 1989.
[78] Koniges A E, Anderson D V. ILUBCG2: A preconditioned biconjugate gradient routine for the solution of linear asymmetric matrix equations arising from 9-point discretizations[J]. Computer Physics Communications, 1987, 43(2): 297-302.
[79] 赵颖旺. 基于有限元的地下水流速场高精度计算方法与流线可视化[D]. 北京: 中国矿业大学(北京), 2015.
[80] 吴吉春, 曾献奎, 祝晓彬. 地下水数值模拟基础[M]. 北京: 中国水利水电出版社, 2018.
[81] 牛庠均. 现代变分原理[M]. 北京: 北京工业大学出版社, 1992.
[82] 徐长发, 李红. 实用偏微分方程数值解法[M]. 武昌: 华中理工大学出版社, 1990.
[83] Boufadel M C. Nutrient transport in beaches: Effect of tides, waves, and buoyancy[M]. Ohio: University of Cincinnati, 1998.
[84] Boufadel M C, Suidan M T, Venosa A D. A numerical model for density-and-viscosity-dependent flows in two-dimensional variably saturated porous media[J]. Journal of Contaminant Hydrology, 1999, 37(1-2): 1-20.
[85] Boufadel M, Suidan M, Venosa A. Numerical modeling of water flow below dry salt lakes: effect of capillarity and viscosity[J]. Journal of Hydrology, 1999, 221(1-2): 55-74.
[86] Boufadel M C, Suidan M T, Venosa A D, et al. 2D variably saturated flows: Physical scaling and Bayesian estimation[J]. Journal of Hydrologic Engineering, 1998, 3(4): 223-231.
[87] Xiao K, Wilson A M, Li H, et al. Crab burrows as preferential flow conduits for groundwater flow and transport in salt marshes: A modeling study[J]. Advances in Water Resources, 2019, 132: 103408.
[88] Van Der Kamp G, Gale J E. Theory of earth tide and barometric effects in porous formations with compressible grains[J]. Water Resources Research, 1983, 19(2): 538-544.
[89] Li H, Li L, Lockington D. Aeration for plant root respiration in a tidal marsh[J]. Water Resources Research, 2005, 41(6): W06023.
[90] Spiteri C, Cappellen P V, Regnier P. Surface complexation effects on phosphate adsorption to ferric iron oxyhydroxides along pH and salinity gradients in estuaries and coastal aquifers[J]. Geochimica et Cosmochimica Acta, 2008, 72(14): 3431-3445.
[91] Spiteri C, Slomp C P, Charette M A, et al. Flow and nutrient dynamics in a subterranean estuary (Waquoit Bay, MA, USA): Field data and reactive transport modeling[J]. Geochimica et Cosmochimica Acta, 2008, 72(14): 3398-3412.
[92] Spiteri C, Slomp C P, Tuncay K, et al. Modeling biogeochemical processes in subterranean estuaries: Effect of flow dynamics and redox conditions on submarine groundwater discharge of nutrients[J]. Water Resources Research, 2008, 44(2): W02430.
[93] Yeh G-T. On the computation of Darcian velocity and mass balance in the finite element modeling of groundwater flow[J]. Water Resources Research, 1981, 17(5): 1529-1534.
[94] Anderson D, Tannehill J C, Pletcher R H. Computational fluid mechanics and heat transfer[M]. NewYork: Taylor & Francis, 1984.
[95] Wang W, Dai Z, Li J, et al. A hybrid Laplace transform finite analytic method for solving transport problems with large Peclet and Courant numbers[J]. Computers & geosciences, 2012, 49: 182-189.
[96] Herzberg A. Die Wasserversorgung einiger Nordseebader (The water supply on parts of the North Sea coast in Germany)[J]. Journal Gabeleucht ung und Wasserversorg ung, 1901, 44(815-819): 824-844.
[97] Badon-Ghyben W. Nota in Verband met de Voorgenomen Putboering Nabij Amsterdam(Notes on the probable results of well drilling near Amsterdam)[J]. Tijdschr Kon Inst Ing, 1889, 1888/9: 8-22.
[98] Carlston C W. An early American statement of the Badon Ghyben-Herzberg principle of static fresh-water-salt-water balance[J]. American Journal of Science, 1963, 261(1): 88.
[99] Boufadel M C, Suidan M T, Venosa A D, et al. Steady seepage in trenches and dams: effect of capillary flow[J]. Journal of Hydraulic Engineering, 1999, 125(3): 286-294.
[100]Xiao K, Li H, Xia Y, et al. Effects of tidally varying salinity on groundwater flow and solute transport: Insights from modelling an idealized creek marsh aquifer[J]. Water Resources Research, 2019, 55(11): 9656-9672.
[101]Yu S, Wang C, Li H, et al. Field and Numerical Investigations of Wave Effects on Groundwater Flow and Salt Transport in a Sandy Beach[J]. Water Resources Research, 2022, 58(11): e2022WR032077.
[102]Geng X, Boufadel M C, Lee K, et al. Biodegradation of subsurface oil in a tidally influenced sand beach: Impact of hydraulics and interaction with pore water chemistry[J]. Water Resources Research, 2015, 51(5): 3193-3218.
[103]Geng X, Pan Z, Boufadel M C, et al. Simulation of oil bioremediation in a tidally influenced beach: Spatiotemporal evolution of nutrient and dissolved oxygen[J]. Journal of Geophysical Research: Oceans, 2016, 121(4): 2385-2404.
[104]Li H, Boufadel M C, Weaver J W. Tide-induced seawater–groundwater circulation in shallow beach aquifers[J]. Journal of Hydrology, 2008, 352(1-2): 211-224.
[105]Liu S, Li H, Boufadel M C, et al. Numerical simulation of the effect of the sloping submarine outlet-capping on tidal groundwater head fluctuation in confined coastal aquifers[J]. Journal of Hydrology, 2008, 361(3): 339-348.
[106]Walther M, Graf T, Kolditz O, et al. How significant is the slope of the sea-side boundary for modelling seawater intrusion in coastal aquifers?[J]. Journal of Hydrology, 2017, 551: 648-659.
[107]王超月. 潮汐和海浪引起的海岸带含水层地下水动态研究[D]. 北京: 中国地质大学(北京), 2016.
[108]庄振业, 李建华, 仇士华等. 莱州湾东岸的全新世海侵和地层[J]. 海洋湖沼通报, 1987(02): 31-39.