Simulation of Drying‐Rewetting Processes in Numerical Groundwater Models Using a New Picard Iteration‐Based Method

When simulating groundwater flow in unconfined and convertible aquifers using a groundwater model with the block‐centered finite‐difference approach, such as MODFLOW, it frequently encounters drying and rewetting of cells. Although many drying and rewetting simulation methods have been proposed in the past, balancing simulation accuracy and convergence capability all at once is difficult. MODFLOW‐2005, which has second‐order accuracy, employs a trial‐and‐error method, but it suffers from computational instability when large quantities of grid cells are dried. MODFLOW‐NWT adopts the upstream‐weighting approach and Newton iteration method to ensure the stability of the drying and rewetting simulations. However, the upstream‐weighting approach has only first‐order accuracy, and the Newton iteration method is complex to implement because it necessitates the establishment of an additional Jacobian matrix. The methods employed by MODFLOW‐NWT are also available in MODFLOW 6, therefore it inherits both the strengths and weaknesses of MODFLOW‐NWT. In this study, a new method, Picard iteration‐based always active cell (PAAC), is proposed. Similar to MODFLOW‐NWT, the PAAC method also uses dry cells as active cells. The PAAC method, however, does not use the upstream‐weighting approach and has second‐order accuracy. Moreover, it ensures good convergence stability even under the Picard iteration method. In addition to discussing the algorithm, five cases were used to comprehensively compare the simulation effects of the PAAC method with MODFLOW‐2005 and MODFLOW‐NWT, including an analytical solution, repeated drying‐rewetting of multi‐layer grids, pumping well problem, perched aquifer problem and a nearly dry single‐layer grid, which verified the practicability of the PACC method.

In early versions of MODFLOW, such as MODFLOW-96 (Harbaugh & McDonald, 1996), MODFLOW-2000(Harbaugh et al., 2000), and MODFLOW-2005 (Harbaugh, 2005), once a grid cell is dewatered (the saturated thickness is zero), it is declared dry and converted to an inactive cell, and all intercell conductances associated with the dry cell are set to zero, to prevent unsolvable matrix equations during groundwater simulations (McDonald et al., 1992).However, the treatment of dry cells in this manner introduces non-physical artifacts, such as the inadvertent loss of sources associated with the dry cells and the interuption of groundwater flow at these locations (Doherty, 2001;Juckem et al., 2006;Merritt, 1994).These artifacts may lead to unacceptable deviations in water balance and groundwater head simulation results (Bedekar et al., 2012;Feinstein et al., 2012).
On the other hand, solving problems of the natural world often required the ability to reactivate dry cells when the water table rose above the dry cell bottom and early versions of the MODFLOW model employed the trial-anderror method to achieve the goal.This method determines whether the dry cell can be rewetted based on whether the hydraulic heads of adjacent grid cells around the dry cell meet the predetermined rewetting conditions (Harbaugh, 2005).However, repeated drying and rewetting attempts of grid cells in critical regions may lead to the model failing to converge (Hunt & Feinstein, 2012;Painter et al., 2008).
Current methods for preventing or addressing problems with the drying and rewetting of grid cells in groundwater models can be summarized into two categories.The first has no physical mechanism and is more of a modeling technique, such as increasing the tolerance of the model convergence and lowering the bottom of key areas (Painter et al., 2008), adjusting the vertical discretization of the model (Keating & Zyvoloski, 2009;Markstrom et al., 2008), assigning minimal saturation thicknesses to dry cells (Doherty, 2001;Lin et al., 2010), converting the nearly dry unconfined domain into a confined domain to simulate (Fitts, 2018;Fitts et al., 2015), etc., which avoid the drying of the simulation domain.In some cases, this method can produce an acceptable simulation effect, but such nonphysical adjustments also raise concerns about the accuracy and reliability of the simulation results (Hunt & Feinstein, 2012;Painter et al., 2008).The second method has a physical mechanism.It uses the dry cell as the active cell and considers the finite-difference flow equation of the dry cell to solve the drying and rewetting simulation problem based on the groundwater flow dynamics.Currently, this method is mainly applied by the upstream-weighting approach, with representative models being MODFLOW-NWT and the current core version, MODFLOW 6 (Hughes et al., 2017;Niswonger et al., 2011;Painter et al., 2008).
The upstream-weighting approach is an alternative method of calculating intercell horizontal hydraulic conductance (Niswonger et al., 2011) and can be written as follows (Taking two adjacent cells in the row direction as an example, with a similar procedure along the column direction): where ∆C is the column width (L); KR ave is the intercell average hydraulic conductivity (L/T, harmonic, logarithmic, or arithmetic mean); ∆R is the distance (L) between the center points of these two adjacent grid cells; h up is the hydraulic head (L) of the grid cell with the relatively higher hydraulic head (i.e., the upstream cell) among the two grid cells; Top up and Bot up are the top and bottom elevations (L) of the upstream cell, respectively.
Equation 2 means using the saturated thickness of the upstream cell to calculate the intercell horizontal hydraulic conductance, thereby avoiding groundwater flowing out of the dry cells, which would not be realistic and could cause model convergence failure, and allowing groundwater flowing into the dry cells from adjacent wet cells (Niswonger et al., 2011;Painter et al., 2008).Consequently, this method allows dry cells to remain active without losing any sources.To enhance convergence in nonlinear groundwater simulations, the MODFLOW-NWT model incorporates the Newton iteration method, coupled with an enhanced preconditioning and acceleration PCG solver χMD solver, and an advanced under-relaxation scheme (Niswonger et al., 2011).It should be acknowledged that the MODFLOW-NWT model has successfully addressed the convergence challenges associated with simulating the drying-rewetting problems.However, it also faces several issues.First, the early versions of the MODFLOW model provided algorithms for intercell conductance, including the harmonic mean method and other optional methods, all of which had second-order accuracy, but the upstream weighting method has first-order accuracy (Edwards, 2011).Second, compared with the Picard iteration method, the Newton iteration method needs to establish a Jacobian matrix corresponding to the groundwater dynamics equations, and smoothing of hydraulic conductance, water storage terms, and some source/sink terms of the model in order to obtain continuous derivatives, which significantly increases the programming complexity of the model (Lipnikov et al., 2016;Maina & Ackerer, 2017;Niswonger et al., 2011).It should be noted that the methods employed by MODFLOW-NWT are also integrated into the current core version of MODFLOW, known as MODFLOW 6.As a result, MODFLOW 6 inherits both the strengths and weaknesses of MODFLOW-NWT.
In this study, a new physically-based drying and rewetting simulation algorithm, namely the Picard iterationbased always active cell (PAAC) method, is proposed for the groundwater model.Consistent with MODFLOW-NWT, the PAAC uses dry cells as active cells and demonstrates convergence comparable to MODFLOW-NWT, even with the Picard iteration method and a general PCG solver.It also achieves secondorder accuracy, similar to MODFLOW-2005.Furthermore, owing to the adoption of the Picard iteration method, the programming complexity of implementing this approach is also relatively low.In this study, the PAAC method is not coupled with the MODFLOW-2005, MODFLOW-NWT or MODFLOW 6 codes, but integrated into the COMUS code.The full, non-open-source version of COMUS, along with Chinese documentation, has been formally published at http://www.hydrotools.top/hydroModel/SkyCoMuS,and the open-source part of the software, along with the English documentation associated with this paper for calculations, has been published on Zenodo at https://doi.org/10.5281/zenodo.10807403(Lu et al., 2024).COMUS utilizes the blockcentered finite-difference flow equations, the modified incomplete Cholesky preconditioning PCG solver, and stress-period-based data management, which are identical to those of MODFLOW-2005, with the primary differences being its development using the C++ programming language and the incorporation of an original algorithm for simulating lake-groundwater interaction (Lu et al., 2021;Wu et al., 2023).However, aside from the PAAC method, this paper does not involve the differences in simulation principles between COMUS and MODFLOW-2005.Therefore, we do not expect the absence of development of the PAAC method in MODFLOW-2005 to hinder a fair comparison between the PAAC method and the two drying-rewetting simulation methods used by MODFLOW-2005 and MODFLOW-NWT.We also do not expect this absence to affect the conclusions drawn in this paper.
The following sections first introduce the principles of the PAAC method and elaborate on its technical similarities and differences with MODFLOW-NWT in handling drying-rewetting problems.Following this, case #1 will be employed to compare the simulation accuracy of COMUS-PAAC with that of MODFLOW-2005 and MODFLOW-NWT.Subsequently, case #2, 3, 4, 5 will be utilized to assess the convergence and self-consistency of simulation results among the three models in modeling drying-rewetting problems.Finally, the characteristics and advantages of the PAAC method will be discussed.

Principle of PAAC Method
Except for those initially designated as inactive grid cells (i.e., mesh elements outside the simulation area), all other cells are active regardless of whether they are dewatered or not during the entire simulation period, and the finite-difference flow equations of these cells must be established to simulate their hydraulic heads.The hydraulic head of a grid cell that is higher than its bottom elevation is referred to as a "wet cell" while the hydraulic head that is lower than its bottom elevation is referred to as a "dry cell."

Horizontal Flow Calculation
With flow entering the cell considered positive, the PAAC method calculates the horizontal flow (L 3 /T ) between any two active grid cells along the row direction as follows (similarly along the column direction): where MaxB is the relatively higher bottom elevation (L) in these two adjacent grid cells; h up is the hydraulic head (L) of the grid cell with the relatively higher hydraulic head (i.e., the upstream cell) among the two grid cells; h down is the hydraulic head (L) of the grid cell with the relatively lower hydraulic head (i.e., the downstream cell) among the two grid cells; CR is the intercell horizontal conductance (L 2 /T ).
The intercell horizontal conductance CR is calculated as follows: where ∆C is the column width (L); KR ave is the intercell average hydraulic conductivity (L/T ); ∆R is the distance (L) between the center points of these two adjacent grid cells; HT ave is the average effective saturated thickness (L).
The intercell average hydraulic conductivity KR ave is calculated using the harmonic mean method as follows: KR up and KR down are the hydraulic conductivity (L/T ) of the upstream and downstream units, respectively; ∆R up and ∆R down are the row width (L) of the upstream and downstream units, respectively.
The average effective saturated thickness HT ave is linearly weighted by the upstream and downstream cell effective saturated thicknesses (Figure 1) as follows: where HT up and HT down are the upstream and downstream cell effective saturated thicknesses (L), respectively; Top up and Top down are the top elevations (L) of the upstream and downstream units, respectively; Bot up and Bot down are the bottom elevations (L) of the upstream and downstream cells, respectively.
Equations 3 and 4 specify eight groundwater flow scenarios between two adjacent grid units in the same layer with different bottom elevations and hydraulic head conditions (Figure 2).Table 1 presents the expressions of intercell conductance calculated by COMUS-PAAC for these eight conditions.In the cases depicted in Figures 2a, 2b, and 2e, the intercell conductance provided by MODFLOW-NWT is zero and the computational expression corresponds to the non-zero portion of Equation 2 in all other scenarios.As shown in Figure 2 and Table 1, when use the PAAC method, groundwater cannot flow horizontally out of any dry cell (Figures 2a, 2b,  and 2e), but it can flow into a dry cell from the horizontally adjacent wet cell when the conditions for bottom elevation and hydraulic head are met (Figures 2c and 2g).
There are two significant differences between COMUS-PAAC and MODFLOW-NWT in calculating horizontal flow.First, COMUS-PAAC simultaneously considers the effective saturated thickness of upstream and downstream cells for weighted calculation of intercell conductance, achieving higher-order accuracy than the upstream weighting method.Second, COMUS-PAAC introduces the MaxB condition, leading to differences in their flow calculation modes.This means that in the case of Figure 2f, MODFLOW-NWT calculates intercell conductance as non-zero value, indicating flow between the grid cells, while PAAC considers that there is no flow between the grid cells due to the head of the wet cell not higher than the bottom elevation of the dry cell (here is the MaxB).In the cases depicted in Figures 2c, 2d, and 2g, if h down is not higher than MaxB, PAAC considers that, in those cases, the flow from the upstream cell to the downstream cell depends on MaxB rather than h down , and the original head difference h down h up in flow calculation is modified to MaxB h up (Equation 3).
Furthermore, it should be noted that the flow calculated by Equation 3 and intercell conductance calculated by Equations 4-6, can maintain the continuity relationship with the hydraulic heads of the upstream and downstream cells, which is one of the prerequisites for model-convergence.

Vertical Hydraulic Conductance Calculation
The interblock vertical hydraulic conductance is involved in multilayer aquifer simulations, with PAAC using the harmonic mean method similar to MODFLOW-NWT to calculate it.The difference lies in the fact that the vertical conductance in MODFLOW-NWT is calculated as the conductance of two one-half cells in a series with continuous saturation between them and remains constant (Niswonger et al., 2011).However, the PAAC method specifies that the aquifer thickness of the lower grid cell is consistently utilized to compute the vertical conductance regardless of whether that grid is dewatered or not, while no such limitation applies to the upper grid cells.

Storage Calculation
Similar to MODFLOW-NWT, the active grid cells under the PAAC method can also be divided into three states: wet-confined, wet-unconfined, and dry according to the relative position between the hydraulic head and the top/ bottom elevation.For transient simulations, the storage coefficient, specific yield, and zero-value storage coefficient would be used to calculate the water storage terms for grid cells under wet-confined, wet-unconfined, and dry states, respectively (Figure 3) and this is the same approach as used in MODFLOW-NWT.Note that storage term conversion is implemented where cells alternate between different states, and more details about this process could be found in the specification of MODFLOW-2005(Harbaugh, 2005).

Vertical Flow Correction Under Perched Groundwater Conditions
Perched groundwater occurs in multilayer aquifer simulations when a cell is unconfined while the cell directly above is fully or partially saturated, and the vertical flow is not dependent on the head of the lower grid cell under this condition, requiring vertical flow correction.PAAC and MODFLOW-NWT employ similar methods for vertical flow correction, and technical details can be referenced from the MODFLOW-2005 specification (Harbaugh, 2005).It should be noted that, as described in Section 2.2, the vertical conductance used by PAAC and MODFLOW-NWT for vertical flow correction is different.

Drying Resistant Cell
In the PAAC method, the active grid cells in the deepest layer and the active grid cells directly above the inactive grid cells cannot be dewatered, which means that the heads of these grid cells cannot be lower than their bottom elevations; otherwise, model-convergence failure may occur, as it does in MODFLOW-NWT (Niswonger et al., 2011;Painter et al., 2008).In MODFLOW-NWT, when the head calculated by the current iteration of a drying resistant cell falls below the cell-bottom altitude, the arithmetic average of the head from the previous iteration and the cell-bottom altitude is used as the updated head value of that cell for this iteration.Nonetheless, when this method is combined with the Picard iteration method, it may result in severe numerical oscillations.Therefore, the PAAC framework includes a  Water Resources Research 10.1029/2022WR034334 new head updating method (Figure 4), which can significantly improve the numerical stability.
where h m,E is the updated head (L) of the drying resistant cell for this iteration; Bot is the bottom elevation (L); h m 1 is the simulated head (L) from the previous iteration; and h m is the originally simulated head (L) for this iteration.
In Equation 7, when h m < Bot and h m ≈ Bot, that is, the originally simulated head is only slightly lower than the bottom elevation, and then the updated head is submitted to h m 1 > h m,E > Bot and h m,E ≈ Bot, which means that the updated head is only slightly higher than the bottom elevation.This is mainly considered because the head would not be corrected if it were slightly higher than the bottom elevation.Therefore, when the head is only slightly lower than the bottom elevation, the updated head should be close to the bottom elevation to maintain the continuity of head solutions to some extent.
In Equation 7, when h m << Bot, means the originally simulated head is much lower than the bottom elevation, the updated head is submitted to h m 1 > h m,E > Bot and h m,E ≈ h m 1 , which would result in slower convergence.To overcome this shortcoming, the head updating is further restricted as follows: which means the updated saturated thickness can only reach maximum 0.6 times of the saturated thickness from the previous iteration.
Furthermore, the PAAC method shares the same characteristics as the MODFLOW-NWT upstream-weighting approach, namely, that when the head of an upstream drying resistant cell approach the cell bottom, the horizontal conductance along the downstream direction is close to zero, thereby preventing flow from leaving the unit.
But the inflow to a nearly dry drying resistant cell from neighboring wet cells would not be affected, therefore, in the absence of other sinks of groundwater, the drying resistant cell can always maintain a head that higher than its bottom elevation, theoretically.
In transient simulations and the majority of the steady-state simulations, the new head updating method can lead to model convergence success.However, in some particularly demanding steady-state simulations, extensive head modifications for a significant number of drying resistant cells may be required within a single iteration, posing challenges for the model to achieve convergence.As a solution, PAAC supports special treatments for the steady-state finite difference flow equation of such cells: where m, m 1 represents the two consecutive iterations; λ is the user-specified retardation factor between 0.0001 and 0.001.The role of λ(h m i,j,k h m 1 i,j,k ) is to add a small storage term to the finite difference flow equation of drying resistant cells during each iteration which approximates the steady-state simulation in a manner similar to transient simulation, thereby effectively improve convergence.When the simulation successfully converge, h m i,j,k ≈ h m 1 i,j,k , and therefore the impact of this storage term on water balance is minimal.The Newton iteration method typically exhibits stronger convergence properties, thereby MODFLOW-NWT does not employ any special treatments for the steady-state finite difference flow equation of such drying resistant cells.

Treatment of the Boundary Conditions
In the PAAC method, dry cells are regarded as active cells that can receive sources and sinks.Generally, sources can be applied to dry cells without any issues, such as recharge.However, model convergence failure when applying sinks to dry cells, as they violate the principle of water mass balance.Therefore, PAAC incorporate specific treatments to handle sinks.
On one hand, sinks related to groundwater head only require a data reasonability check.Taking the River Package as an example, groundwater is discharged into the river when the river stage is below the groundwater head.If there is bottom elevation information in the grid cell, as long as the user-input river stage is reasonable and consistently higher than the bottom elevation of the grid cell, there won't be any computational issues.However, if the user-input river stage is unreasonable, PAAC will generate an error message.Similar sinks include the General-head Package, Drain Package, Stream Package, Lake Package, and Reservoir Package.
On the other hand, sinks unrelated to groundwater head require specialized treatment.These types of sinks primarily involve Well Package.In the PAAC method, the pumping rates are linearly decreased as the head drops below a user-specified saturated thickness threshold: where Ф is the user-specified saturated thickness threshold of the well cell.Generally, it is recommended to set the value of Ф between 0.1 and 0.25 times the thickness of the grid cell for easier convergence.Q Well is the userspecified pumping rate (L 3 /T ) and Q Pump is the applied pumping rate (L 3 /T ).MODFLOW-NWT also employs a similar approach for managing pumping wells.However, to satisfy the requirement of smooth equations with continuous derivatives in the application of the Newton iteration method, the reduction of pumping rates is not accomplished through a linear attenuation formula, but rather through the utilization of a cubic formula.

Under-Relaxation Iterative Schemes
To improve model-convergence for highly nonlinear problems or simulations involving a large number of drying resistant cells whose head values must be determined by repeated modifications, it is necessary to activate underrelaxation.
There are two under-relaxation mechanisms in the PAAC framework: the global relaxation mechanism and the cell-by-cell relaxation mechanism.These two mechanisms can be adopted independently or in combination.Generally, the global relaxation mechanism is used first, followed by the cell-by-cell relaxation mechanism when the global relaxation mechanism fails to achieve convergence.
For all active cells, the new head values for the current iteration after the global relaxation mechanism are computed according to: where h m are the original head values (L) for this iteration, obtained directly by solving the matrix equation, h m 1 are the head values (L) from the previous iteration, and Damp is the user-specified global relaxation factor between 0.0 and 1.0.Damp = 1.0 indicates that there is no global relaxation performed, and the smaller the value, the deeper the global relaxation.

Water Resources Research
10.1029/2022WR034334 LU ET AL.
A simplified and improved adaptive learning algorithm (delta-bar-delta technique) is used for cell-by-cell relaxation (Smith, 1993): where ∆h m,T is the difference (L) between the head value for the current iteration after global relaxation and the head value from the previous iteration of the grid cell, ω m (0 < ω < 1) is the adaptive relaxation factor of the grid cell in the current iteration, and the smaller the value, the deeper the cell-by-cell relaxation; h m,C is the relaxed head value (L) for the current iteration of the grid cell after cell-by-cell relaxation.
The adaptive relaxation factor ω m is calculated as: Non-oscillation lasts for every NIT times where ω m 1 is the adaptive relaxation factor value from the previous iteration.The decreasing and increasing coefficients α and δ of ω are defined by the user, generally with constant values of 0.35-0.95and 1-5, respectively; β is the user-specified single-step increasing factor of ω, typically with constant values range from 0.0 to 0.2.With the same δ value, the smaller the α and β values, the deeper the cell-by-cell relaxation.ω decrease or increase automatically based on whether the solution oscillates over nonlinear iterations.Numerical oscillation occurs when the change direction of the calculated head value in this iteration is opposite to that in the previous iteration: where ∆h m 1,T is head change value (L) from the previous iteration.
Equation 13 indicated that if the head solution oscillates, the adaptive relaxation factor ω would decrease proportionally to the value of α; and whenever the number of consecutive non-oscillations accumulates more than NIT times, ω would increase according to the values of δ and β.Based on current application experience, the parameter combination of α = 0.7, δ = 3, β = 0.01, and NIT = 5 can be used, similar to the "Simple" parameter combination in MODFLOW-NWT, to solve problems with relatively low nonlinearity, while α = 0.7, δ = 3, β = 0.001, and NIT = 8 is similar to the "Complex" parameter combination in MODFLOW-NWT, and can be used to solve problems with high nonlinearity.If parameter tuning is necessary to achieve better convergence, adjusting NIT and β can be attempted to influence the speed of ω increasement and improve convergence.
Although the cell-by-cell relaxation schemes in PAAC and MODFLOW-NWT are both based on the delta-bardelta technique, there are significant differences between them.PAAC introduces two new parameters, δ and NIT.
The role of δ is to enable ω recover quickly at an exponential rate, thereby improve the convergence speed and reduce the water budget error.The role of NIT is to enhance convergence.

Test Simulations
Five distinct simulation cases are presented to demonstrate some of the features and capabilities of the PAAC method.Case 1 is designed to examine the accuracy and error decay rate of the intercell horizontal conductance algorithm in PAAC method.Case 2 demonstrates the effect of the intercell horizontal conductance algorithm in PAAC on maintaining the self-consistency of drying-rewetting simulation, as well as the reasonability of the eight horizontal groundwater flow scenarios derived by such algorithm.elucidates the reasonability of the flow calculation formula proposed by PAAC.Case 5 confirms the effectiveness of the handling approach for drying resistant cells (Equations 7-9) and the modified under-relaxation iterative schemes in PAAC.For the purpose of comparison, simulations were also conducted using both MODFLOW-2005 (Harbaugh, 2005) and MODFLOW-NWT (Niswonger et al., 2011).The modified incomplete Cholesky preconditioning PCG solver was adopted for both MODFLOW-2005 and COMUS-PAAC, while the enhanced preconditioning and acceleration PCG solver, namely the χMD solver was employed by MODFLOW-NWT.MODFLOW-2005 supports a number of methods for calculating intercell conductance, with the most commonly used harmonic mean method being used in this study.Similarly, in MODFLOW-NWT, the intercell conductivity is also computed using the widely employed harmonic mean method.The head change criterion for convergence (HCLOSE) and flow change criterion for convergence (RCLOSE) of the three models are set to 0.0001 m and 0.0001 m 3 /d, respectively.Furthermore, the test simulations were conducted on a Windows 10 64-bit system with an Intel Core i9 8950HK processor and 48 GB of RAM.

Case #1: Comparison to One-Dimensional Analytical Solution
This case, originally employed in the MODFLOW-NWT specification to compare the simulation accuracy between the upstream weighting method in MODFLOW-NWT and the harmonic mean method in MODFLOW-2005 (Niswonger et al., 2011), is utilized in this study for a fair comparison of the simulation accuracy among PAAC's approach for calculating intercell horizontal conductance, the harmonic mean method, and the upstream weighting method.The total length of the model domain along rows was 5,000 m and the grid cells along columns were 50 m wide.The bottom and top altitudes of the model were 0 and 50 m, respectively.There are two discretization modes along the row directions: "equal interval" and "unequal interval" (Figure 5), each consisting of sparse grid conditions (20 grid cells) and dense grid conditions (40 grid cells).The row width of grid cells decreasing by a ratio of 0.92 from left to right for the unequal interval sparse grid (Figure 5b) and the dense grid was obtained by halving the grid size of the sparse grid.Hydraulic conductivity of grid cells was 50 m/d.Constant heads of 10 m at the left end of the model and 50 m at the right end of the model and the rest are variable head cells.The steady-state analytical head (m) and flow (m 3 /d) for this case are (Niswonger et al., 2011):   where h L and h R are the constant heads (m) at the left and right ends of the model, respectively.D is the distance between two constant head cells (m); x is the distance from h R at any situation in the model domain (m); K is the hydraulic conductivity (m/d); ∆C is the width of grid cells along the column directions (m).
Figure 6 and Table 2 show the differences between the simulation results of various models.In this case, the simulation principles involved in the three models differ only in the algorithm for intercell horizontal conductance, thus it is the sole reason for the differences in simulation results.For the equal interval grid configuration, whether simulating under sparse grid conditions or dense grid conditions, the simulated heads of MODFLOW-NWT had the lowest accuracy, attributed to the lower differential accuracy of the upstream weighting formulation.Under both sparse and dense grid conditions, the average head errors are 0.365 and 0.187 m, respectively, with corresponding flow errors of 25.687 and 12.42 m³/d.The accuracy of MODFLOW-2005 was better than that of MODFLOW-NWT, with average head errors of 0.072 and 0.017 m, along with flow errors of 2.5 and 0.581 m³/ d under the two grid conditions, respectively.By halving the grid size, MODFLOW-NWT exhibits a decrease of 48.7% and 51.6% in average head error and flow error, respectively, owing to the first-order accuracy of the upstream weighting method.MODFLOW-2005 experiences reductions of 76.3% and 76.7%, respectively, attributed to the second-order accuracy of the harmonic mean method.COMUS-PACC has the highest accuracy, with simulated heads and flows exactly match the analytical solutions, even under sparse grid conditions.For the equal interval grid configuration, the row width of any grid cell is equal to the distance between the center points of any adjacent grid cells, denoted uniformly as ∆R.Let h 1 and h 2 represent the hydraulic heads of any two adjacent cells, with the assumption that h 1 < h 2 .When simulating with COMUS-PAAC, the intercell horizontal flow q R is given by (Equations 3-6, expressed as a positive value): Equations 17 and 16 have identical forms.In other words, the flow between any two adjacent cells calculated by COMUS-PAAC matches the analytical flow.Therefore, under the equal interval grid configuration, regardless of the number of simulated cells, the simulation results of COMUS-PAAC will be essentially consistent with the analytical solution.
For the unequal interval grid configuration, the simulation accuracy of all three models decreases, but COMUS-PACC is still the least affected, owing to the closer alignment between the flow calculated by PAAC and the analytical solution (Table 2).By halving the grid size, the average head error and flow error of MODFLOW-NWT decreased by 46.5% and 49.9%, respectively.For MODFLOW-2005, the reductions were 75.9% and 76.7%, and for COMUS-PAAC, they were 83.4% and 80.9%, respectively.The error decay rate of COMUS-PAAC is slightly faster than that of MODFLOW-2005 but significantly faster than that of MODFLOW-NWT.This is because the intercell horizontal conductance algorithm in PAAC, similar to the harmonic mean method in MODFLOW-2005, has second-order accuracy, while the upstream weighting method in MODFLOW-NWT has only first-order accuracy.

Case #2: Simulating the Drying-Rewetting Problems of Multi-Layer Model Cells
The two-dimensional transient groundwater simulation of earth dam profile was used as case 2. The model domain is 40 m high and 100 m wide, consists of a grid with 40 layers, one row, 100 columns, and the width of each of the 2,193 active grid cells in the row and column directions is 1 m.There are two grid discretization modes: "regular discretization" and "vertically offset discretization" (Figure 7).The cell thickness is 1 m under the regular discretization scenario and the cell thickness of adjacent grid cells on the same layer is three times different with the vertically offset discretization scenario.The simulation period is divided into three 10-day stress periods, each with a constant time step length of one day.General-head boundary was implemented on the left side of the model domain, and the head values were set to 38, 22, 38 m during the three stress periods in order to simulate the complete process of the water level before the earth dam from falling to recovering to its initial value.It should be noted that this case involves multi-layer aquifer simulations, and besides the algorithm for intercell horizontal conductance, variations in the intercell vertical conductance algorithm can also impact the simulation results.Nevertheless, both MODFLOW-NWT and COMUS-PAAC employ a similar harmonic mean method for calculating intercell vertical conductance, with minor variations in computation details.Therefore, the disparities in simulation results between the two models are expected to primarily stem from differences in the intercell horizontal conductance algorithm.
For both regular and vertically offset discretization scenarios, COMUS-PAAC, MODFLOW-2005, and MODFLOW-NWT all provide reasonable solutions for the transient groundwater simulation with repeated drying-rewetting processes of multi-layer model cells.There is no need to activate the under-relaxation options for both COMUS-PAAC and MODFLOW-2005 in order to achieve convergence.Furthermore, the utilization of the "Simple" parameter combination for the Newton solver adequately guarantees convergence in MOD-FLOW-NWT.
For COMUS-PAAC, on the 11th day following a sudden fall of general-head, the regular and vertically offset discretization scenarios introduced 388 and 375 dry cells, respectively, which accounted for 17.7% and 17.1% of the total number of active cells.Likewise, on the 21st day after a sudden rise of general-head, the two scenarios resulted in a reduction of dry cells by 560 and 553, respectively, representing 25.5% and 25.2% of the total number of active cells.The simulated water tables for the 11th day with the lowest general-head is shown in Figure 8.The simulated water tables for each of the three independent models are similar in both regular and vertically offset discretization scenarios.The water tables simulated by COMUS-PAAC and MODFLOW-NWT are very similar because both models employ dry cells as active cells, which benefit from the upstream-weighting method and the intercell horizontal conductance algorithm in PAAC, effectively preventing unrealistic outflows from dry cells while enabling inflows from adjacent wet cells.However, the simulated water table by MODFLOW-2005 is significantly lower than the results obtained from COMUS-PACC and MODFLOW-NWT, due to differences between the intercell horizontal conductance algorithms (harmonic mean method vs. upstream weighting method and the new approach proposed by PAAC).Additionally, it may also be associated with the introduction of nonphysical artifacts in handling dry cells in MODFLOW-2005.
The distribution of heads simulated by the three models is affected by the grid discretization mode to varying degrees (Figure 8).It can be seen that the distribution of head contours under the vertically offset discretization scenario of MODFLOW-NWT and MODFLOW-2005 differs significantly from that under the regular discretization scenario, and the difference becomes more prominent the closer it is to the general-head boundary.In comparison, the head contour distribution simulated by COMUS-PAAC under the two grid discretization scenarios is relatively similar, suggesting that considering the effective saturated thickness of both upstream and downstream cells for intercell conductance calculation helps maintain self-consistency in the results.It also indicates that the eight horizontal groundwater flow conditions designed by PAAC are reasonable.For MODFLOW-NWT, which calculates the intercell conductance solely based on the saturated thickness of the upstream cell, its simulation results are evidently affected by the chosen grid discretization method.The nonphysical artifacts introduced by MODFLOW-2005 when handling dry cells can impact the self-consistency of simulation results.However, it is unclear whether the harmonic mean method itself has issues in maintaining the self-consistency of simulation results.
The number of outer iterations and total run time required for convergence, along with the water balance output from three models are presented in Table 3.Under the vertically offset discretization scenario, the main source/ sink terms (general-head inflow and drain outflow) simulated by MODFLOW-2005 and MODFLOW-NWT are about 25% and 12% underestimated than those under the regular discretization scenario, respectively, while only about 1.5% lower for COMUS-PAAC.The water balance comparison conclusions of each model are consistent with those of groundwater head, further illustrating that COMUS-PAAC can provide relatively stable simulation results under different meshing methods.It can be observed that the COMUS-PAAC requires significantly fewer outer iterations for convergence compared to MODFLOW-2005 and MODFLOW-NWT.The run time for COMUS-PAAC falls between that of MODFLOW-2005 and MODFLOW-NWT.The fewer outer iterations but longer run time of COMUS-PAAC compared to MODFLOW-2005 is due to the fact that COMUS-PAAC simulates the dry cells as active cells, resulting in a relatively higher workload for each iteration.

Case #3: Simulation Involving Pumping Wells
Case 3 was derived from case 2, the vertically offset discretization scenario, by incorporating pumping conditions that are particularly prone to convergence failure.The two pumping wells are applied in the grid cells at the 37th layer, 21st column, and 37th layer, 60th column (Figure 7b).thickness of the grid cell, both COMUS-PAAC and MODFLOW-NWT will automatically reduce the specified pumping rate.
The simulation using MODFLOW-2005 failed to converge due to the repeated transition of the well cells between dewatered and wetted states, while convergence was achieved using MODFLOW-NWT by utilizing the "Simple" Newton solver parameter combination.In the case of COMUS-PAAC simulation, convergence was achieved by activating the cell-by-cell relaxation mechanism with the recommended parameter values, Damp = 1, α = 0.7, δ = 3, β = 0.001, NIT = 8.Continuing to use the simulation results from the 11th day as an example, it can be observed that the simulated head of COMUS-PAAC is generally higher than that of MODFLOW-NWT (Figure 9, this pattern holds true for other simulation time steps as well).Moreover, both COMUS-PAAC and MODFLOW-NWT simulated well extraction volumes fall below the specified amount (6000 m 3 ), with COMUS-PAAC showing a slightly higher extraction volume compared to MODFLOW-NWT (Table 4).Notably, COMUS-PAAC achieves convergence with approximately half the number of outer iterations needed by MODFLOW-NWT, while both models exhibit similar run time requirements (Table 4).
The similarity in simulation results between COMUS-PAAC and MODFLOW-NWT indicates the usability of PAAC in managing well boundary conditions.It is worth noting that, owing to the Picard iteration method only requiring equation continuity, not differentiability, PAAC can handle well boundary in a straightforward continuous piecewise Equation 10 compared to MODFLOW-NWT (the corresponding treatment formulas can be found in Niswonger et al., 2011).

Case #4: Perched Groundwater Simulation
The challenging nonlinear two-dimensional perched groundwater simulation was presented by Bedekar et al. (2012).The model domain (Figure 10) is 2,000 m high and 10,000 m wide, consists of a grid with 10 layers, one row, 100 columns, and the widths of each grid cell in the row, column and layer directions are 100, 100, 200 m, respectively.An aquitard is located in the central portion of layer 4, with a hydraulic conductivity of 0.00001 m/d, while the hydraulic conductivity in the remaining part of the model domain is 10 m/d.Constant   2012) conducted a series of attempts using MODFLOW-2000 to address this problem.Their approach involved adjusting the rewetting parameters of the trial-and-error method or employing upstream weighting, smoothing conductance, backtracking, and under-relaxation techniques simultaneously, etc.However, achieving successful convergence is contingent upon the initial heads under the Picard iteration method, and typically, the perched groundwater cannot be successfully modeled.Similarly, in this study, the simulation using MODFLOW-2005 also proved unsuccessful in simulating perched groundwater, hence no simulation results using MODFLOW-2005 were provided.In the case of MODFLOW-NWT, model convergence is achieved by utilizing the "Simple" Newton solver parameter combination.On the other hand, COMUS-PAAC only requires a slight global relaxation (Damp = 0.95) to achieve convergence.
As shown in Figure 10a, the simulated heads in the bottom four layers by MODFLOW-NWT are generally lower than those by COMUS-PAAC.It should be noted that in the case of MODFLOW-NWT, the simulated thickness of perched groundwater in the grid cells at both ends of the perched layer approaches zero, while in other perched layer grid cells, it is significantly lower compared to COMUS-PAAC (Figure 10b).In fact, even with a tenfold increase in the recharge rates, MODFLOW-NWT consistently simulates perched groundwater thickness close to  Water Resources Research 10.1029/2022WR034334 zero in the grid cells at both ends of the perched layer.This is due to the fact that the head of a dry cell approaches the head in the topmost directly underlain cell at the row/column location, which is at least partially saturated (Hunt & Feinstein, 2012), as occurred in the columns outside the perched layer.Consequently, there is an original head difference of approximately 800 m between the upstream grid cells at both ends of the perched layer and the downstream adjacent grid cell in the same layer outside the perched layer, resembling the scenario depicted in Figure 2c.This substantial head difference renders the grid cells at both ends of the perched layer incapable of retaining groundwater.However, in such situations, COMUS-PAAC modifies the head difference to MaxB h up (Equation 3).As a result, the perched groundwater thickness in COMUS-PAAC varies correspondingly with the recharge rate, indicating that the flow calculation method in PAAC may be more reasonable compared to MODFLOW-NWT.
On the other hand, COMUS-PAAC requires fewer outer iterations to converge compared to MODFLOW-NWT, while the total run time of COMUS-PAAC is approximately twice as long as that of MODFLOW-NWT (Table 5).This is because the χMD solver of MODFLOW-NWT has faster execution speed (Ibaraki et al., 2021).Nevertheless, the computational efficiency of COMUS-PAAC is still satisfactory, considering its relatively short actual execution time.Compared to the study by Bedekar et al. (2012), under the same conditions of utilizing the Picard iteration method, as well as the upstream weighting method or the intercell horizontal conductance algorithm in PAAC, both of which consider dry cells as active cells and contribute to avoiding numerical oscillations, the easier convergence of PAAC indicates the significant role of the flow calculation mechanism in PAAC in ensuring convergence.

Case #5: Unconfined Aquifer That Is Extremely Close to Drying up
The prototype of Case 5 serves as a representative highly nonlinear unconfined simulation employed in the MODFLOW-NWT specification to validate the robust convergence of MODFLOW-NWT (Niswonger et al., 2011).In this study, the row and column widths of the prototype were diminished by a factor of five, resulting in a 25-fold increase in the number of drying-resistant cells, thereby requiring much more additional repeated modifications (Equations 7-9) to determine the head values of such cells, consequently further significantly enhancing the nonlinearity.The inclusion of Case 5 aims to substantiate the convergence of PAAC in extreme, highly nonlinear simulations encompassing a substantial number of drying resistant cells.The model consists of a thin, unconfined aquifer with an area of 64 km 2 , which was divided into 400 rows and 400 columns with 20 m dimensions along both row and column directions, resulting in a large-scale computation involving 160,000 valid grid cells (Figure 11).The bottom altitude varies between 4 and 80 m, and the top altitude is 200 m.The horizontal hydraulic conductivity is 1 m/d.Constant head of 24 m at 75 cells with the lowest bottom altitude (Figure 11), and the rest are variable head cells with a certain amount of recharge.
This case includes four recharge scenarios: high recharge rates (HRR), medium recharge rates (MRR), low recharge rates (LRR), and ultra-low recharge rates (ULRR).For the HRR scenario, grid cell recharge rates range between 4 × 10 7 -7 × 10 6 m/d, with a total recharge of about 233.723 m 3 /d.Recharge rates for MRR, LRR, and ULRR scenarios are 1/10, 1/100, and 1/1,000 of the HRR scenario, amounting to recharge volumes of 23.3723, 2.3372, and 0.2337 m 3 /d, respectively.Recharge for all scenarios is distributed according to the bottom-altitude distribution, where higher recharge rates correspond to higher altitudes.Initial heads in the model are 20 m above the cell bottom-altitudes except at the constant-head cells, and the simulation consists of a single steady-state stress period.
The drastic variation of the aquifer-bottom altitude across the model domain makes it easy for a large number of drying resistant cells to approach the drystate, as indicated and verified by the MODFLOW-NWT specification (Niswonger et al., 2011).Even under the HRR scenario, the solutions provided by MODFLOW-2005 are very sensitive to the initial heads, therefore the comparison of model performance on this case is only carried out around MODFLOW-NWT and COMUS-PACC.
For this case, the head values of a large number of drying resistant cells close to dry-state must be determined by repeatedly modifications, and the probability of convergence failure increases as the recharge rates decrease.MODFLOW-NWT has to deploy the "Complex" Newton solver parameter combination, which is specifically tailored to address highly nonlinear problems, alongside adjusting the THICKFACT parameter to accelerate convergence.Likewise, for COMUS-PAAC to achieve convergence, it requires the activation of the cell-by-cell relaxation mechanism and the implementation of the retardation factor λ (see in Equation 9), along with tuning the NIT and λ parameters.A summary of the control parameters for both MODFLOW-NWT and COMUS-PAAC models can be found in Table 6.
Figure 12 shows the simulated head and saturated thickness (hydraulic head minus bottom elevation) of COMUS-PAAC under four recharge scenarios, as well as the hydraulic head calculated by MODFLOW-NWT for comparison.The spatial distribution of the head simulated by COMUS-PAAC and MODFLOW-NWT is generally similar under HRR and MRR scenarios, but there are still slight differences in areas with relatively slow terrain slope.However, with the further decrease of recharge rates, most portions of the model domain have heads that are extremely close the cell bottom, leading to the simulated heads of the two models is being nearly identical under LRR and ULRR scenarios.
Both COMUS-PAAC and MODFLOW-NWT provided good mass balance in each of the four scenarios (Table 7).Additionally, while COMUS-PAAC requires more outer iterations compared to the MODFLOW-NWT, the relative computational efficiency of the two models varies depending on the specific scenarios and does not exhibit a systematic pattern.Given that MODFLOW-2005, which also employs the Picard iteration method, failed to successfully simulate this problem, it suggests that the adapted approach for handling drying-resistant cells (Equations 7-9) and the modified under-relaxation schemes in PAAC contribute significantly to improving its convergence.

Discussion and Conclusion
Simulation of the drying-rewetting processes is one of the most challenging numerical problems of a groundwater model (Al-Muqdadi et al., 2020;Hunt & Feinstein, 2012;Painter et al., 2008).Although many simulation methods have been proposed in the past, balancing simulation accuracy and convergence capability all at once is difficult.For example, the trial-and-error method employed by MODFLOW-2005 and other earlier MODFLOW versions required iterative parameter tuning to identify suitable parameters, and it frequently encountered challenges when simulating complex problems, such as non-convergence and discrepancies with actual conditions (Feinstein et al., 2012;Juckem et al., 2006;Reilly, 2001).MODFLOW-NWT and the current core version MODFLOW 6 employed the upstream-weighting approach and Newton iteration method with a strong mathematical and physical mechanism, and it has successfully addressed the convergence challenges associated with simulating the drying-rewetting problems, which can be considered as a landmark achievement.However, the upstream-  weighting approach, after all, only has first-order accuracy (Edwards, 2011).Despite numerous studies demonstrating that upstream-weighting approach applied to groundwater numerical modeling can achieve "acceptable accuracy" which derived from the comparison of various typical examples (Forsyth et al., 1995;Niswonger et al., 2011;Sun & Yeh, 1983;Zaidel, 2013), the majority of users continue to rely on the traditional MODFLOW method to ensure the relative accuracy in majority of the groundwater modeling cases.For special applications, upstream-weighting method applied to groundwater numerical modeling can achieve certain accuracy but with room for improvement.Moreover, the Newton iteration method necessitates that source/sink terms transition smoothly with continuous derivatives, making it difficult to extend the model-function from a code-development perspective.
In this study, the proposed PACC method exhibits several notable characteristics and advantages.First, it considers the effective saturated thickness of the upstream and downstream cells simultaneously for weighted calculation of the intercell horizontal hydraulic conductance, thus ensuring stronger self-consistency when simulating drying-rewetting problems on different kinds of grids and achieving second-order accuracy.Second, it preserves the continuity of each computational term of the finite-difference flow equation based on more sensible physical mechanisms, such as intercell horizontal conductance calculation and flow calculation, improving convergence of drying-rewetting simulations using the Picard iteration method and the general PCG solver.Third, it incorporates the modified under-relaxation schemes, which reducing solver iterations and enabling the simulation of highly nonlinear problems.Fourthly, it modifies the hydraulic head updating method for drying resistant cells and the steady-state finite-difference flow equation for such cells, significantly enhancing the model convergence in particularly demanding simulations.Lastly, it couples with the Picard iteration method, which results in a relatively minor impact on the programming of various source/sink packages.In most cases, straightforward modifications to the source/sink packages are adequate to accommodate this approach.
The major disadvantage of various existing simulation methods for drying-rewetting problems is that it is difficult to ensure convergence under the Picard iteration method, which was also the main reason for combining the Newton iteration method (Niswonger et al., 2011;Painter et al., 2008).However, the PACC method overcomes this shortcoming and even could be applied to highly nonlinear simulations involving a large number of grid cells, such as case 5, that the MODFLOW-NWT also necessitates the utilization of the "Complex" Newton solver parameter combination to enhance convergence.Furthermore, based on the tested cases, the relative computational efficiency between COMUS-PAAC (based on the Picard iteration method) and MODFLOW-NWT (based on the Newton iteration method) is problem specific and remains inconclusive, even in large scale computational simulations.While the code size of the open-source version of COMUS-PAAC used in this study is much smaller than that of MODFLOW-NWT, considering that the majority of total run time is consumed by the solver, the overhead from the model outside of the PAAC method is relatively small.Therefore, we believe that the difference in code size would not significantly impact the comparison results in computational efficiency.It is worth noting the good convergence and satisfactory computational efficiency that can be compared to MODFLOW-

Figure 1 .
Figure 1.Schematic diagram of calculating the average effective saturated thickness.

Figure 2 .
Figure 2. Groundwater flow between two adjacent grid cells in the same layer using the PAAC method.The dotted line indicates no flow, the solid line indicates flow, and the arrow indicates the potential groundwater flow direction.

Figure 3 .
Figure 3. Model cell for which storage factors are used during different states.

Figure 4 .
Figure 4. Schematic diagram of the head updating method of the drying resistant cell.

Figure 5 .
Figure 5. Case #1.Schematic diagram of the model domain and boundary conditions under (a) equal interval grid configuration and (b) unequal interval grid configuration.

Figure 6 .
Figure 6.Analytical head for case #1 and head, head error as calculated by three models under (a) equal interval sparse grid configuration, (b) equal interval dense grid configuration, (c) unequal interval sparse grid configuration, and (d) unequal interval dense grid configuration.
The general-head boundary conductance was set to 10,000 m 2 /d.The drain boundary was defined on the right side, the drain elevation of the two grid discretization scenarios is same (at the grid center with the regular discretization scenario), and the drain conductance is 1 m 2 /d.The abrupt fluctuations in general head and drainage result in repeated drying and rewetting of numerous grid cells, allowing for testing of the eight groundwater flow relationships in PAAC under vertically offset discretization scenario (Figure7b).Hydraulic conductivity, initial heads, specific storage, and specific yield are 4 m/d, 38 m, 0.0001 and 0.08, respectively.In particular, for MODFLOW-2005, it is necessary to activate the "conversion from dry cell to wet cell" option and set the cell wetting parameters of the trial-and-error method, in this case WETFCT = 0.1, IWETIT = 1, IHDWET = 1, WETDRY = 0.01 m.A smaller absolute value of the empirical parameter WETDRY and a lower value of the empirical parameter IWETIT make the conversion of dry cells into wet ones easier, thereby contributing to the reduction of non-physical artifacts introduced in handling dry cells in MODFLOW-2005.On the other hand, with the same WETDRY and IWETIT values, the parameters WETFCT and IHDWET have minimal impact on the simulation results once the model successfully converges.Therefore, the aforementioned parameter settings are considered reasonable for this case.

Figure 7 .
Figure 7. Case #2.Side view of the model domain and boundary conditions under (a) regular discretization scenario and (b) vertically offset discretization scenario.
During the three stress periods, the pumping rate for each well remains constant at 100 m³/d.If the saturated thickness of the well cell falls below 0.1 times the

Figure 8 .
Figure 8. Case #2.Comparisons of water tables and groundwater heads simulated by three models on day 11.

Figure 11 .
Figure 11.Case #5.Schematic diagram of the model domain and boundary conditions.

Figure 12 .
Figure 12.Case #5.Comparisons of groundwater heads simulated by MODFLOW-NWT and COMUS-PAAC, and saturated thickness simulated by COMUS-PAAC under (a) high recharge rates scenario, (b) medium recharge rates scenario, (c) low recharge rates scenario, and (d) ultra-low recharge rates scenario.

Table 1
The Intercell Horizontal Hydraulic Conductance Calculated by COMUS-PAAC Under Eight Groundwater Flow Scenarios Bot up ,Top up Bot up ] ⋅ ΔR up +min[max(0, h down Bot up ), Top down Bot down ] ⋅ ΔR down Bot down ,Top up Bot up ] ⋅ ΔR up +min[h down Bot down ,Top down Bot down ] ⋅ ΔR down ΔR up +ΔR down

Table 2
Comparison of Solutions of Three Models for Case #1

Table 3
Comparison of Water Balance, Outer Iterations and Total Run Time Provided by Three Models for Case #2

Table 4
Comparison of Water Balance, Outer Iterations and Total Run Time Provided by MODFLOW-NWT and COMUS-PAAC for Case #3

Table 5
Comparison of Water Balance, Outer Iterations and Total Run Time Provided by MODFLOW-NWT and COMUS-PAAC for Case #4 LU ET AL.

Table 7
Comparison of Water Balance, Outer Iterations and Total Run Time Provided by MODFLOW-NWT and COMUS-PAAC forCase #5 when using PAAC method for simulating drying-rewetting problems without compromising selfconsistency or increasing model programming complexity.Simultaneously, it achieves second-order accuracy comparable toMODFLOW-2005. NWT