New insights into the differences between the dual node approach and the common node approach for coupling surface–subsurface flow

The common node approach and the dual node approach are two widely applied approaches to coupling surface–subsurface flow. In this study both approaches are analyzed for cell-centered as well as vertex-centered finite difference schemes. It is shown that the dual node approach should be conceptualized and implemented as a one-sided first-order finite difference to approximate the vertical subsurface hydraulic gradient at the land surface. This results in a consistent dual node approach in which the coupling length is related to grid topology. In this coupling approach the coupling length is not to be interpreted as a nonphysical model parameter. Although this particular coupling approach is technically not new, the differences between this consistent dual node approach and the common node approach have not been studied in detail. In fact, this coupling scheme is often believed to be similar to the common node approach. In this study it is illustrated that in comparison to the common node approach, the head continuity at the surface–subsurface interface is formulated more correctly in the consistent dual node approach. Numerical experiments indicate that the consistent dual node approach is less sensitive to the vertical discretization when simulating excess infiltration. It is also found that the consistent dual node approach can be advantageous in terms of numerical efficiency.

Abstract.The common node approach and the dual node approach are two widely applied approaches to coupling surface-subsurface flow.In this study both approaches are analyzed for cell-centered as well as vertex-centered finite difference schemes.It is shown that the dual node approach should be conceptualized and implemented as a one-sided first-order finite difference to approximate the vertical subsurface hydraulic gradient at the land surface.This results in a consistent dual node approach in which the coupling length is related to grid topology.In this coupling approach the coupling length is not to be interpreted as a nonphysical model parameter.Although this particular coupling approach is technically not new, the differences between this consistent dual node approach and the common node approach have not been studied in detail.In fact, this coupling scheme is often believed to be similar to the common node approach.In this study it is illustrated that in comparison to the common node approach, the head continuity at the surface-subsurface interface is formulated more correctly in the consistent dual node approach.Numerical experiments indicate that the consistent dual node approach is less sensitive to the vertical discretization when simulating excess infiltration.It is also found that the consistent dual node approach can be advantageous in terms of numerical efficiency.

Introduction
There is a variety of hydrogeological problems, such as the hydrologic response of hillslopes and river catchments, which requires an integrated analysis of surface and sub-surface flows.This has led to the development of physically based, distributed parameter models for simulating coupled surface-subsurface flows.Well-known examples of such models include MODHMS (Panday and Huyakorn, 2004), InHM (Ebel et al., 2009), HydroGeoSphere (Therrien et al., 2010), CATHY (Camporese et al., 2010), WASH123D (Yeh et al., 2011), ParFlow (Kollet and Maxwell, 2006) and OpenGeoSys (Kolditz and Shao, 2010).Typically, subsurface flow is governed by the Richards' equation whereas surface flow is either governed by the kinematic wave or the diffusive wave equation.
The coupling between subsurface and surface flow may be either based on the common node approach (Kollet and Maxwell, 2006) or on the dual node approach (Ebel et al., 2009;Panday and Huyakorn, 2004;VanderKwaak, 1999).In the common node approach, coupling is formulated by a continuity in head between surface and subsurface nodes.The dual node approach is based on formulating an exchange flux between the surface and subsurface nodes.Typically, the dual node approach is conceptualized as a hydraulic separation of the surface and the subsurface by an interface with a given thickness (Liggett et al., 2012).The thickness of this interface defines a coupling length between the dual nodes to formulate the discrete exchange flux between the dual nodes.
It has been argued that the coupling length represents a nonphysical model parameter, because there is often no evidence to support the existence of a distinct interface between the two flow domains (Kollet and Maxwell, 2006).As such, it appears that the common node approach is a more physically based coupling approach (Kollet and Maxwell, 2006;Liggett et al., 2012).It has also been found that accurate simulations based on the dual node approach typically require Published by Copernicus Publications on behalf of the European Geosciences Union.5710 R. de Rooij: Coupling surface-subsurface flow a very small coupling length (Ebel et al., 2009;Liggett et al., 2012Liggett et al., , 2013)).Since it is known that the dual node approach mimics the common node in the limit as the coupling length goes to zero (Ebel et al., 2009), it thus seems that the dual node approach is most accurate if it mimics the common node approach.Nonetheless, it has been argued that the dual node approach remains an attractive alternative coupling approach since it offers more flexibility than the common node approach.Namely, while it can mimic the common node approach, the dual node approach offers the possibility to simulate a less tight coupling of surface-subsurface flow which results in increased computational efficiency (Ebel et al., 2009).
In this study a detailed analysis of both coupling approaches is provided for cell-centered as well as vertexcentered finite difference schemes.This analysis starts with the crucial observation that the topmost subsurface nodal values as computed by the finite difference schemes represent the mean values within the topmost discrete control volumes.Numerical experiments to compare the coupling approaches are carried out with the model code DisCo (de Rooij et al., 2013b).It is shown that the dual node approach should be interpreted and implemented as a one-sided finite difference approximation of the vertical hydraulic gradient at the land surface.This yields a consistent dual node scheme in which the coupling length is defined by the half the thickness of the topmost subsurface cells.The scheme of An and Yu (2014) as well as the scheme of Kumar et al. (2009) are essentially very similar to this consistent dual node scheme.In the work of Panday and Huyakorn (2004), one of the suggestions to define the coupling length is to use half the thickness of the topmost subsurface cells, which yields a consistent dual node scheme.While the idea that the coupling length can be based on the grid topology is not new (Panday and Huyakorn, 2004), the idea that it must be related to grid topology to obtain a consistent approach is a significant new insight.Namely, since the coupling length in the consistent dual node approach is not to be interpreted as the thickness of a layer that separates the subsurface from the surface, the consistent dual node approach is not automatically less physically based than the common node.In fact, as explained in this study, in comparison to the common node approach the implementation of a head continuity at the surface-subsurface interface is formulated more correctly in the consistent dual node approach.
The current consensus about how the dual node approach compares to the common node approach is based on alternative dual node approaches which, as explained in this study, are different from the consistent dual node approach.In this study the consistent dual node approach is compared in detail with the common node approach.It is shown that if the vertical discretization is sufficiently fine, then the common node approach and the consistent dual node approach are equally accurate.However, when simulating excess infiltration the consistent dual node approach is found to be less sensitive to the vertical discretization in comparison to the common node approach.This advantage in accuracy is related to the fact that head continuity is more correctly formulated in the consistent dual node approach.Moreover, it is also shown that the consistent dual node approach can be advantageous in terms of numerical efficiency when simulating runoff due to both excess saturation as well as excess infiltration.The finding of this study show that the consistent dual node approach compares more positively with respect to the common node approach than other dual node approaches.

Interpretation of nodal values
As explained later on, a correct interpretation of nodal values is crucial for understanding the dual and common node approach for coupling surface-subsurface flow.Moreover, both coupling approaches depend on the configuration of surface and topmost subsurface nodes near the land surface.This configuration depends on whether cell-centered or vertex-centered schemes are used.In this study both type of schemes will be covered, but for simplicity only finite difference schemes are considered.
In both cell-centered as vertex-centered schemes the flow variables such as the heads and the saturation are computed on nodes.In vertex-centered schemes these nodes coincide with the vertices of the mesh, whereas in cell-centered schemes the nodes coincide with the cell centers.When employing a finite difference scheme, nodal values correspond to the mean value within surrounding discrete control volumes.In cell-centered finite difference schemes these discrete volumes are defined by the primary grid cells.In vertexcentered finite difference schemes these discrete volumes are defined by the dual grid cells.Ideally, the mean values in the discrete control volumes are derived by applying the midpoint rule for numerical integration such that their approximation is second-order accurate.Therefore, the nodal values should ideally represent values at the centroid of the surrounding discrete control volume (Blazek, 2005;Moukalled et al., 2016).In that regard, a cell-centered finite difference scheme is thus more accurate than a vertex-centered finite difference scheme.Namely, in cell-centered finite difference schemes the nodal values always correspond to the centroids of the cell whereas in vertex-centered finite difference schemes nodes and centroids (of the dual cells) do not coincide at model boundaries and in model regions where the primary grid is not uniform.It is well known that this mismatch between nodes and centroids can lead to inaccuracies since the mean values within affected discrete volumes are not computed by a midpoint rule (Blazek, 2005;Moukalled et al., 2016).
Typically, vertex-centered schemes for simulating coupled surface-subsurface flow are based on mass-lumped finite element schemes (Liggett et al., 2012) and not on finite difference schemes.However, with respect to coupling surface-Hydrol.Earth Syst.Sci., 21, 5709-5724, 2017 www.hydrol-earth-syst-sci.net/21/5709/2017/ R. de Rooij: Coupling surface-subsurface flow 5711 subsurface flow there is actually no difference between a mass-lumped finite element scheme and a vertex-centered finite difference scheme.Similar to those in vertex-centered finite difference schemes, the nodal values in mass-lumped finite element schemes define the mean values inside dual grid cells (Zienkiewicz et al., 2005).Moreover, the coupling approaches establish one-to-one relations between surface and topmost subsurface nodes which do not depend on whether a finite difference or a finite element approach is being used.

Common node approach
The common node approach defines a head continuity between the topmost subsurface nodes and the surface nodes.This continuity requires that the topmost subsurface nodes and the surface nodes are colocated at the land surface such that there is a continuity in the elevation head.This requirement is automatically fulfilled in vertex-centered schemes.Figure 1a illustrates the configuration of common nodes in ParFlow, a cell-centered scheme (R. Maxwell, personal communication in relation to previous work of the author, 2011; de Rooij et al., 2013a).Figure 1c illustrates the configuration of common nodes for vertex-centered schemes.This configuration is similar to the configuration used in HydroGeo-Sphere (Therrien et al., 2010).
Considering that nodal values ideally represent the mean values within discrete control volumes, as described in Sect.2, it can be argued that the head continuity as implemented in the common node approach is not in agreement with the physical principle of head continuity at the land surface.Namely, the common node approach enforces a continuity between surface heads at the land surface and the mean subsurface heads within the topmost subsurface discrete control volumes which have a finite thickness.This is different from enforcing a continuity between surface heads and subsurface heads within an infinitesimally thin subsurface layer directly below the land surface.As such, the common node approach is only numerically correct if the topmost subsurface cells are very thin.
4 Dual node approach 4.1 Basics Figure 1b and c illustrate the classical arrangement of surface and subsurface nodes in cell-centered and vertex-centered finite difference schemes, respectively.Commonly, the dual node approach is expressed in terms of an exchange flux q e (L T −1 ) computed as follows (Liggett et al., 2012;Panday and Huyakorn, 2004): where h s and h ss are the hydraulic heads (L) associated with the surface node and the topmost subsurface node, respectively, f p (-) is the fraction of the interface that is ponded and l is the coupling length (L).The ponded fraction of the interface is typically defined by a function that varies smoothly between zero at the land surface elevation and unity at the rill storage height which defines the minimum water depth for initiating lateral overland flow (Panday and Huyakorn, 2004).In Eq. (1) the term f p K z / l is commonly referred to as the first-order exchange parameter, where first-order means that the exchange flux depends linearly on the hydraulic head difference.Typically, Eq. ( 1) is not derived as a numerical approximation of basic flow equations that govern the exchange flux, but is merely presented a numerical technique to couple two different flow domains (Ebel et al., 2009;Liggett et al., 2012).Subsequently, the dual node approach is conceptualized by interpreting Eq. (1) as an expression that describes groundwater flow across a distinct interface separating the two flow domains (Ebel et al., 2009;Liggett et al., 2012Liggett et al., , 2013)).In the following, it is illustrated that the dual node approach can and should be derived from basic equations that describe infiltration into a porous medium.Using Darcy's law, the infiltration rate at the ponded land surface q s→ss (L T −1 ) can be written as a function of the vertical subsurface hydraulic gradient at the land surface: where h the hydraulic head (L), z the elevation head (L), k r the relative hydraulic conductivity (-) K z the saturated vertical hydraulic conductivity (L T −1 ) and z s the elevation head at the land surface.The relative hydraulic conductivity is unity because Eq. ( 2) applies to the ponded land surface which implies fully saturated conditions at the land surface (i.e., ponding means p s > 0, where p s is the pressure head at the surface).Similarly, the infiltrability (L T −1 ), defined as the infiltration rate under the condition of atmospheric pressure (Hillel, 1982), can be written as follows: The relative hydraulic conductivity is again unity because the saturation equals unity under atmospheric conditions (p s = 0).The infiltration rate at nonponded land surface q atm→ss (L T −1 ) can be expressed as follows: where q R is the effective rainfall rate (i.e., the infiltration rate is limited by either the infiltrability or the available effective rainfall rate).The total exchange flux across the surfacesubsurface interface can now be written as follows: To approximate the vertical subsurface hydraulic gradient in Eqs. ( 2) and ( 3), it is crucial to recognize that, according to the principle of head continuity at the land surface, the surface hydraulic head at a surface node must also represent the subsurface head at the land surface at that location.Moreover, since the subsurface hydraulic heads at the topmost subsurface nodes are ideally associated with the centroids of the topmost subsurface discrete control volumes, these head values do not represent values at the land surface but at some depth below the land surface.Because the subsurface hydraulic heads at the dual nodes can be and should be associated with a different elevation, the vertical subsurface head gradient between the dual nodes can be approximated by a standard finite difference approximation.If this approximation is being used to approximate the gradient at the land surface in Eqs. ( 2) and (3), then this approximation is by definition a one-sided first-order finite difference.By defining the coupling length by l = z dn , where z dn is the difference in the mean elevation head associated with the dual nodes, the infiltration rate and infiltrability can thus be computed with the following one-sided finite difference approximation: The above definition of the coupling length l = z dn ensures a proper approximation of the vertical gradient in elevation head at the land surface: The above derivation of the consistent dual node approach from basic flow equations has implications for how the dual node approach is conceptualized and how it should be implemented.The idea that the coupling length must be directly related to the spatial discretization is an important new insight.Namely, as the coupling length is related to grid topology, it does not represent a nonphysical parameter associated with a distinct interface separating the two domains.It is also crucial to observe the difference between the consistent dual node approach and the common node approach regarding how the head continuity at the surface-subsurface interface is formulated.As explained in Sect.2, the formulation in the common node approach is only correct if the topmost subsurface discrete volumes are very thin.In comparison, the formulation in the dual node approach is correct irrespective of the vertical discretization.Namely, irrespective of the vertical discretization the surface hydraulic heads equal the subsurface heads at the interface.
Since nodal values in cell-centered scheme are located at the centroids of the cells, the coupling length is simply given by l = z s − z ss , where z s and z ss are the elevation heads (L) associated with the surface node and the topmost subsurface node, respectively.This value for the coupling length in cell-centered schemes has also been suggested by Panday and Huyakorn (Panday and Huyakorn, 2004).However, in their work, the particular advantage of choosing this value (i.e., maintaining a unit gradient in elevation head) is not recognized.The coupling schemes used by An and Yu (2014) and Kumar et al. (2009) are also in essence consistent dual node schemes.However, these schemes are not recognized as a dual node scheme.Instead, An and Yu (2014) argue that their scheme is similar to the common node approach of Kollet and Maxwell (2006).Kumar et al. (2009) argue that their scheme is similar to the dual node approach if the coupling length goes to zero, which implies that their scheme would be similar to the common node approach.However, contrary to the common node approach the schemes of An and Yu (2014) and Kumar et al. (2009) Kumar et al. (2009) are actually quite different from the common node approach.As already mentioned, the consistent dual node scheme differs from the common node approach with respect to how the head continuity is formulated at the surface-subsurface interface.As discussed later on, this difference has crucial consequences in terms of accuracy as well as numerical efficiency.
In vertex-centered schemes the commonly used nodal configuration near the surface is such that z s = z ss .However, even though the topmost subsurface node is located at the land surface in a vertex-centered scheme, the elevation head at this node should ideally correspond to the mean elevation head within the topmost subsurface discrete volume.This suggests that the topmost subsurface node should be moved to the centroid of the topmost subsurface discrete volume.Although this is a possible solution, the drawback of this solution is that the subsurface model ceases to be a purely vertex-centered scheme.Moreover, such an operation cannot be performed in finite element schemes since the nodal positions define the geometry of the elements.Therefore, an alternative solution is proposed.Namely, in vertex-centered schemes the elevation of the surface nodes are changed according to z s = z ss +l, where l is equal to half the thickness of the topmost subsurface dual cell.The resulting nodal configuration is illustrated in Fig. 1d.When applying this solution, all the topmost subsurface cells must have the same thickness, such that the topography is increased with the same value everywhere.In essence, the motivation behind this solution is that a more accurate approximation of the hydraulic gradient (i.e., enforcing a unit gradient in elevation head) is more important than the actual elevation of the land surface.Similar to the nodal configuration in ParFlow, the resulting nodal configuration may not seem ideal.Namely, the surface elevation does not coincide with the top of the subsurface grid.Nonetheless, as illustrated later on, simulation results obtained with the resulting scheme are reasonable.
To illustrate that the presented dual node approach exhibits consistent behavior, the necessary conditions for ponding due to excess infiltration and exfiltration are considered.In general, ponding starts when q R > I (Hillel, 1982).Observing that Eq. ( 6) defines the computed infiltrability when p s = 0 and that the gradient in elevation head between the dual nodes is unity, the infiltrability can be expressed by Ponding due to excess infiltration occurs if q R /K z > 1 and implies that saturation in the subsurface starts from the top down (Hillel, 1982).Using q R /K z > 1, it follows from Eq. ( 8) that p ss is still negative at the moment of ponding.This is reasonable, because the pressure head value at the topmost subsurface node represents a value at a certain depth below the land surface.Top-down saturation implies that saturation at the topmost subsurface node occurs after pond-ing and thus a negative pressure head value at this node at the moment that ponding starts.It is noted that if the ratio q R /K z is greater than but close to unity or if the coupling length is very small, then this condition becomesp ss ≈ 0. Once ponding starts the total flux rate between the dual nodes equals K z ((p s − p ss ) / l + 1).Top-down saturation requires that this flux exceeds the vertical hydraulic conductivity.Reaching saturation at the topmost node (p ss = 0) therefore requires p s ≥ 0. Thus, while p ss is still negative at the moment that ponding starts, saturation at the topmost subsurface node will occur some time after ponding started.Ponding due to excess saturation occurs if q R /K z < 1 and implies that saturation in the subsurface starts from the bottom up (Hillel, 1982).It follows from Eq. ( 8) that ponding due to excess saturation occurs while p ss > 0. Thus, ponding starts after reaching fully saturated conditions at the topmost subsurface node, which is again reasonable.It is noted that if the ratio q R /K z is smaller than but close to unity or if the coupling length is very small, then ponding occurs when p ss ≈ 0.

Comparison to alternative coupling approaches
To illustrate that it is crucial to account for the meaning of the values at the topmost subsurface nodes, it is instructive to consider what happens if these values are not taken as the mean values within discrete control volumes.As a first example, consider vertex-centered schemes where the dual nodes are defined such that z ss = z s , as illustrated in Fig. 1c.This is inconsistent because it defines a zero gradient in elevation head between the dual nodes.Since the vertical gradient in elevation head between the dual nodes is zero, the total flux rate after ponding now equals K z (p s − p ss ) / l.Topdown saturation requires that this flux exceeds the vertical hydraulic conductivity.Thus, reaching saturation at the topmost subsurface node (p ss = 0) requires p s > l.Therefore, top-down saturation will not occur if runoff occurs and if the surface water depths remain smaller than the chosen coupling length.Indeed, it has been pointed out in other studies that the coupling length should be smaller than the rill storage height (Delfs et al., 2009;Liggett et al., 2012).The zero vertical gradient in elevation head between the dual nodal also means that the ponding occurs when p ss > −lq R /K z .This implies that ponding due to excess saturation occurs while the topmost subsurface node is not yet saturated.This dual node approach has been compared to the common node approach in vertex-centered schemes (Liggett et al., 2012).
A second example is the dual node approach for cellcentered schemes as implemented in MODHMS, which uses an adapted pressure-saturation relationship for the topmost subsurface nodes such that the topmost subsurface node only becomes fully saturated if hydraulic head at the node rises above the land surface (Liggett et al., 2013).Since the topmost subsurface heads are associated with the cell centroid, this dual node scheme defines a unit gradient in elevation head at the land surface.However, the saturation value at www.hydrol-earth-syst-sci.net/21/5709/2017/ Hydrol.Earth Syst.Sci., 21, 5709-5724, 2017 5714 R. de Rooij: Coupling surface-subsurface flow the topmost node is associated with a location at the land surface and not with the centroid of a discrete control volume.This has undesirable consequences.Namely, saturating the topmost subsurface node (p ss = l) due to excess infiltration requires that p s > l.Indeed, when simulating excess infiltration with MODHMS, a very small coupling length is needed to simulate top-down saturation due to excess infiltration (Gaukroger and Werner, 2011;Liggett et al., 2013).It can also be shown that ponding due to excess saturation occurs while p ss > 0. But, because of the adapted pressuresaturation relationship this means that ponding starts while the topmost subsurface node is not yet saturated.This dual node approach has been compared to the common node approach in cell-centered schemes (Liggett et al., 2013).
The two comparison studies of Liggett et al. (2012Liggett et al. ( , 2013) indicate that the dual node approach is typically only competitive with the common node approach in terms of accuracy once the coupling length is very small.However, the requirement for a very small coupling length, is a logical consequence if the topmost subsurface nodal values are not taken as the mean values within discrete volumes.In essence, by choosing a very small coupling length this inconsistency is minimized.This contrasts with the consistent dual approach in which decreasing the coupling length for a given vertical discretization will result in more inaccurate simulation results as this would be numerically incorrect.
CATHY (Camporese et al., 2010) as well as the model of Morita and Yen (2002) are examples of models which are neither based on the common node approach, nor a dual node approach.Both these models are conjunctive models in which the surface and subsurface flow are computed separately in a sequential fashion and in which coupling is established by matching the flow conditions along the surfacesubsurface interface.A complete discussion is outside the scope of this paper, but it is worthwhile to mention that these models share some crucial characteristics with the consistent dual node approach.Although the two models are different, both models switch between appropriate boundary conditions along the surface-subsurface interface, such that infiltration fluxes are limited to the infiltrability.In both models the infiltration fluxes are computed while accounting for the unit vertical gradient in elevation head near the surfacesubsurface interface.In addition, in both models ponding occurs when the infiltrability is exceeded.

Numerical model
To compare the consistent dual node approach with respect to the common node approach in terms of accuracy and computational efficiency numerical experiments are presented.These experiments are carried out with the model code DisCo.This model code can simulate coupled surface-subsurface flow with the dual node approach using a fully implicit or monolithic scheme (de Rooij et al., 2013b).Subsurface flow is governed by the Richards' equation while surface flow is governed by the diffusive wave equation.
Starting from a dual node scheme, the implementation of a common node scheme is relatively straightforward.If the surface nodes are numbered last, a permutation vector can be constructed which gives the corresponding topmost subsurface node for each surface node.Then, the node numbering used in the original dual node scheme can still be used to compute the surface and subsurface flow terms.Subsequently, using the permutation vector, the surface and subsurface flow terms associated with a common node can be combined into the same row of the global matrix system.In addition, when using the common node approach, there is no need to evaluate exchange flow terms between the two flow domains.It is noted that the surface flow and subsurface flow computations are exactly the same irrespective of the coupling approach.As such, the model permits comparison between the two approaches in terms of accuracy as well as numerical efficiency.
An adaptive error-controlled predictor-corrector one-step Newton scheme (Diersch and Perrochet, 1999) is used in which a single user-specified parameter controls the convergence as well the time stepping regime.Although this scheme may not be necessary the most efficient scheme, it ensures that the time discretization error is the same irrespective of the applied coupling approach.For brevity, further details about the model are not discussed here and can be found elsewhere (de Rooij et al., 2013b).

Hillslope scenarios
The model code is applied to a set of three hillslope scenarios.Table 1 lists the abbreviations used in the figures to distinguish between the coupling approaches, and to distinguish between cell-centered and vertex-centered schemes.Each scenario is solved using different but uniform vertical discretizations, and z specifies the discretization of the primary grid.The first two simulation scenarios consider hillslope problems as designed by Sulis et al. (2010).For the purpose of this study, a third scenario is considered in which the initial and boundary conditions are different in order to create a flooding wave across an unsaturated hillslope.The problems consist of a land surface with a slope of 0.05 which is underlain by a porous medium.The domain is 400 m long and 80 m wide.The subsurface is 5 m thick.In the direction of the length and in the direction of the width, the discretization is 80 m.Different vertical discretizations are considered.The van Genuchten parameters are given by s r = 0.2, s s = 1.0, α = 1 m −1 and n = 2.The porosity is 0.4 and the specific storage is 10 −4 m −1 .The Manning roughness coefficients are given by 3.  The first scenario considers excess saturation, and the saturated conductivity equals 6.94 × 10 −4 m min −1 .Figures 2  and 3 illustrate the simulated runoff and the number of Newton steps, respectively.Figures 4 and 5 illustrate the subsurface pressure heads at the topmost subsurface nodes and the water depths on the surface nodes.For the second scenario, which considers excess infiltration, the saturated hydraulic conductivity equals 6.94 × 10 −7 m min −1 .Figures 6 and 7 show the simulated runoff and the number of Newton steps, respectively.Figures 8 and 9 illustrate the subsurface pressure heads at the topmost subsurface nodes and the water depths on the surface nodes for the finest and the coarsest vertical discretization, respectively.In the third scenario a surface water flood wave crossing the hillslope in the downhill direction is simulated by applying a Neumann boundary condition of 1.0 m 3 s −1 for a duration of 200 min to the surface nodes with the highest elevation.The initial water table is located at a depth of 1.5 m.The vertical saturated hydraulic conductivity equals 6.94 × 10 −6 m min −1 .Figure 10 illustrates the differences in simulated runoff and Fig. 11  nodes for the finest and the coarsest vertical discretization, respectively.

Accuracy
As discussed by Ebel et al. (2009) and confirmed by others (Liggett et al., 2012) the dual node approach mimics the common node approach if the coupling length becomes sufficiently small.When comparing the consistent dual node approach and the common node approach a very similar observation applies.If the topmost subsurface cells are very thin, then the coupling length in the consistent dual node approach is very small.Also, if the topmost subsurface cells are sufficiently thin then the formulation of head continuity at the surface-subsurface interface in the common node approach is correct.Thus, the common node approach will mimic the consistent dual node approach.Indeed, the simulation results indicate that a relatively fine vertical discretization yields similar results for the common node approach as well as for the consistent dual node approach (Figs. 2a,4a,6a,8a,10a and 12a).
A relatively fine uniform vertical discretization also enables the simulation of sharp saturation fronts with the Richards equation (Pan and Wierenga, 1995;Ross, 1990).As such, the simulation results based on the finest vertical discretization can be taken as reference solutions that enable comparisons of the coupling approaches when a coarser vertical discretization is used.

Excess saturation
The simulation results of runoff due to excess saturation, as obtained by the common node approach and the consistent dual node approach as depicted in Fig. 2, illustrate that simulating excess saturation runoff is not significantly affected by the vertical discretization.This is because the time needed to reach fully saturated conditions in the subsurface is a simple function of the flow boundary conditions and the initial water content.It is thus expected that the vertical discretization does not significantly affect the simulation of excess saturation.Although the vertical discretization may affect the computed initial water content, this effect is usually negligible.It has been found in other studies that the vertical discretization has little effect on simulated runoff due to excess saturation (Kollet and Maxwell, 2006;Sulis et al., 2010).

Excess infiltration
When simulating excess infiltration the common node approach requires fully saturated conditions at the topmost subsurface node for ponding to occur.However, top-down satu-ration associated with excess infiltration implies that reaching fully saturated conditions in the topmost subsurface discrete volumes should require more time than reaching fully saturated conditions at the land surface, especially if the vertical discretization is relatively coarse.It is thus expected that the common node approach delays runoff and that this delay increases for a coarser vertical discretization.In addition, if the saturation fronts are less sharp due to a relatively coarse vertical discretization, it takes more time to reach saturated conditions at the common node.This will further delay runoff.Indeed, the simulation results indicate clearly that runoff is delayed when using the common node approach, particularly if the vertical discretization is relatively coarse (Figs. 6,9a,10 and 13a).It has also been found in other studies that the common node approach delays runoff due to excess infiltration if the vertical discretization is relatively coarse (Sulis et al., 2010).As explained in Sect.4.2, when using the consistent dual node approach, ponding due to excess infiltration occurs before reaching fully saturated conditions at the topmost subsurface node.More specifically, ponding occurs when the in-  filtrability is exceeded.Compared to the condition for ponding in the common node approach, this is arguably more correct.Namely, if saturation occurs from the top-down then the saturation at a certain depth occurs later than saturation at the land surface.Indeed, simulation results indicate that, when simulating excess infiltration, the consistent dual node approach is less sensitive to the vertical discretization in comparison to the common node approach.This is clearly indi- cated in Figs.6b-d, 9a, 10b-d and 13a.To further explain this difference in accuracy, it is emphasized that the spatial resolution only affects the accuracy of the flow computations when using the consistent dual node approach and that the formulation of head continuity at the interface remains correct.In contrast, when using the common node approach, if the spatial resolution is too coarse then this does not only affect the accuracy of the flow computations but in addition the formulation of head continuity becomes incorrect.It must be emphasized, however, that regardless of the applied coupling approach, the vertical discretization must be relatively fine.As indicated by 9a, the difference between the simulated results and the reference solution increases for a coarser discretization.Eventually such differences will lead to unreasonable results regardless of the coupling approach.
It is interesting to note that An and Yu (2014) also found that their model was less sensitive to the vertical discretization in comparison to ParFlow when simulating runoff due to excess infiltration.Whereas An and Yu (2014) hypothesized that this difference in performance was related to using irregular grids instead of orthogonal grids as in ParFlow, it is argued here that this difference can be explained by the fact that both models use a different coupling approach.
Although the consistent dual node approach is less sensitive to the vertical discretization in comparison to the common node approach, it is useful to explain in detail how the vertical discretization affects the accuracy of the consistent dual node approach to the vertical discretization.A relatively coarse vertical discretization may result in an underestimation of the vertical pressure gradient at the land surface.This is because, in a soil close to hydrostatic conditions, the pressure heads increase with depth.Therefore, the infiltrability during the early stages of infiltration may be underestimated.If the applied flux rate is sufficiently large such that the underestimated infiltrability is exceeded, then runoff during the early stages will be overestimated.Figure 6d illustrates that runoff is indeed overestimated at early times when simulated with the cell-centered scheme, a relatively coarse vertical discretization and a consistent dual node approach.During the later stages of infiltration the pressure head at the topmost subsurface node will be underestimated due to the combined effect of an underestimated infiltration rate and the overly diffused saturation fronts.This results in an overestimation of the infiltration rate in the later stages.Thus, at some time after ponding has started, it is expected that the amount of runoff is underestimated.
If the underestimated infiltrability is not exceeded, then the overly diffused saturation fronts resulting from a relatively coarse vertical discretization will eventually lead to an underestimation of pressure head at the topmost subsurface node, and as such the infiltrability may be overestimated at later times.Consequently, when using the consistent dual node approach, runoff due to excess infiltration may be delayed.However, the delay in runoff as simulated by the consistent dual node approach will only equal the delay in runoff as simulated by the common node approach in the limit when q R /K z goes to unity.Namely, as explained in Sect.4.2, if q R /K z goes to unity, then the consistent dual node approach behaves similarly to a common node approach.However, in general, if the consistent dual node approach delays runoff, this delay will be smaller than the delay in runoff as simulated by the common node approach.
Comparing Figs.12a and 13a, it can be observed that if the vertical discretization is relatively coarse then a common node can act as an artificial barrier for a surface water wave advancing across an initially unsaturated subsurface domain.Namely, as the wave travels downstream the wave can only advance to the next common node once it is fully saturated.The effect of this artificial barrier is that the front of the surface water wave is steepened.In contrast, the consistent dual approach simulates a wave that becomes less steep as it ad-vances downstream for relatively fine as well as relatively coarse vertical discretizations, as depicted in Fig. 13a.
As illustrated in Figs.6b-d and 10b-d, if the coupling approach and the vertical discretization are identical, then the vertex-centered schemes are closer to the reference solution with respect to the cell-centered schemes.This difference results solely from the fact that the primary mesh is the same for both schemes.As such, the vertical extent of the topmost subsurface volumes is twice as small when using the vertexcentered scheme.This difference in vertical grid resolution near the land surface explains the differences between the schemes.

Computational efficiency
The computational efficiency of the schemes is measured in terms of the number of Newton steps.The number of Newton steps equals the number of times that the linearized system of equations is solved, and this number depends on the time step sizes as well as the number of failed Newton steps.It is emphasized that the measured efficiency depends crucially on the applied model code.Nonetheless, as shown in the following, the measured differences in efficiencies can be explained in terms of abrupt changes in how fast pressure heads near the surface-subsurface interface are evolving with time.Regardless of the type of scheme used to solve the nonlinear flow equations, such abrupt changes are difficult to solve.
Once ponding occurs a surface-subsurface flow model will encounter significant numerical difficulties as surface flow terms are activated.In essence, the activation of these terms represents a discontinuity in flow behavior which is challenging to resolve (Osei-Kuffuor et al., 2014).Indeed, the Newton steps as depicted in Figs. 3 and 7 indicate that simulations encounter difficulties at the moment of ponding.These figures also indicate that the consistent dual node approach can be more efficient in comparison to the common node approach.

Excess saturation
Just before the moment of ponding due to excess saturation, the rate of change in pressure heads at the topmost subsurface nodes is relatively high for both coupling approaches.This high rate is related to the shape of the water retention curve.Typically, the derivative of the saturation with respect to the pressure head goes to zero when approaching fully saturated conditions.Once ponding starts, the surface flow terms are activated and therefore the rate of change in pressure heads at the topmost subsurface nodes decreases drastically.Both approaches must handle this drastic change.However, from Figs. 4b and 5b it can be observed that the rate of change decreases more abruptly when using the common node approach.
When using the common node approach the vertical hydraulic gradients in the subsurface are close to zero at the moment of ponding, since additional water volumes can only be accommodated by means of specific storage.This implies that the infiltration rate drops instantaneously at the moment of ponding.In contrast, in the dual node approach, ponding starts when the infiltrability is exceeded.Thus, at the moment of ponding, the infiltration rate is higher in comparison to the common node approach.After ponding this infiltration rate will decrease quickly as the hydraulic heads at the dual nodes equilibrate.This difference in the infiltration rate at the moment of ponding explains why the topmost subsurface hydraulic heads change more smoothly when using the dual node approach.If the vertical discretization is coarser, then the infiltration rate at the moment of ponding, computed with the consistent dual node approach, is even higher and this results in a lower initial rate initial rate of change in water depth, as depicted in Fig. 5a.
The more abrupt changes in pressure heads at the common node in comparison to the changes in pressure heads at the dual nodes mean that solving the activation of ponding with the common node approach is more difficult.It is noted that the differences in the infiltration rates between the two coupling approaches only occur at the moment of ponding and directly thereafter when water depths are relatively small.Namely, quickly after ponding, the hydraulic heads at the dual nodes will equilibrate and after that the two coupling approaches will behave similar.This explains why these differences in infiltration rates do not significantly affect the accuracy of simulated runoff.

Excess infiltration
Figures 8, 9, 12 and 13 illustrate the evolution of pressure heads at dual nodes and common nodes when simulating excess infiltration.When applying the consistent dual approach, the net flux into a topmost subsurface cell will decrease once ponding occurs, because the applied flux rate will be partitioned between dual nodes (i.e., between the surface flow and subsurface flow domain).This occurs while the topmost subsurface node is not yet fully saturated.After ponding the infiltration rate decreases such that if the topmost subsurface node reaches fully saturated conditions the net flux into the topmost subsurface node is relatively small.In contrast, partitioning of the applied flux rate on a common node between the surface flow and subsurface domain starts when the common node reaches fully saturated conditions at this node.This means that just before ponding the rate of change in pressure head is relatively high as the common node is driven towards fully saturated conditions, while the infiltration rate is relatively high.This means that, similar to the excess saturation scenario, the rate of change in pressure head at the common node is high just before ponding.At the moment of ponding, this rate must drop abruptly as surface flow terms are activated.This abrupt change explains why the common node approach is less efficient.
Figures 7 and 11 also indicate that a coarser vertical discretization only provides a significant gain in efficiency in terms of Newton steps when using the consistent dual node approach.When using the common node approach, a coarser discretization does not change the fact that the topmost subsurface node must reach fully saturated conditions for ponding to occur and that the infiltration rate is relatively high just before ponding.When using the consistent dual node approach, a coarser vertical discretization means that the saturation fronts are more diffused such that the flow problem becomes easier to solve.
Figure 8a and 9a illustrate that for the second simulation scenario, ponding occurs almost simultaneously at all the surface nodes.Figure 12a and 13a show that this is different for the third scenario where ponding occurs at different times as the flooding wave travels downstream.When Fig. 11a is compared with Fig. 12a and when Fig. 11d is compared with Fig. 13a, it is clear that the common node approach encounters difficulties around each time ponding starts at a surface node.Figure 11 shows that these difficulties are encountered for all discretizations.In contrast the consistent dual node approach has much fewer difficulties solving these problems.As discussed in Sect.6.1.2,the common node approach may result in steepening the advancing wave.This implies that water depths will be changing more quickly.This presents an additional difficulty for solving this flow problem with the common node approach.

Conclusions
In this study it is shown that the dual node approach should be conceptualized and implemented as a one-sided finite difference approximation of the vertical hydraulic gradient at the land surface.This provides an important new insight into the coupling length.Namely, if the dual node approach is properly implemented then the coupling length is related to the vertical grid resolution.Thus, the coupling length does not represent an additional nonphysical model parameter, and therefore the dual node approach is not automatically a less physically based approach in comparison to the common node approach.Actually, this study shows if the vertical discretization is not sufficiently fine then the head continuity at the surface-subsurface interface is formulated more correctly in the consistent dual node scheme.This difference in formulation has consequences for how both approaches compare in terms of accuracy and efficiency.
Numerical experiments indicate that the consistent dual node approach is equally accurate or more accurate than the common node approach.It has been shown that in comparison to the common node approach the consistent dual node approach is less sensitive to the vertical discretization when simulating excess infiltration.However, the practical advantage of the consistent dual node approach in terms of accuracy is limited.Namely, if the vertical discretization is re-fined, both approaches will converge to more accurate and eventually similar results when simulating excess infiltration.When simulating excess saturation, both approaches yield similar results even if the vertical discretization is relatively coarse.
Nonetheless, even though the advantage of the consistent dual node approach in terms of accuracy is limited, the fact that the consistent dual node approach is equally or more accurate than the common node approach is a significant finding.Namely, this finding is different from the commonly held view that a dual node approach is most accurate if it mimics the common node approach.Moreover, it also illustrates clearly that the consistent dual node approach is not similar to a common node approach.
Numerical experiments indicate that the consistent dual node approach can be more efficient than the common node approach while being equally or more accurate than the common node approach.It has been shown that this difference in efficiency is related to abrupt changes in the evolution of pressure heads around the moment that ponding is initiated.
Based on the findings in this study the models of An and Yu (2014) and Kumar et al. (2009) are expected to have some advantages with respect to models that are based on the common node approach.This is because these models are based on a consistent dual node approach.Moreover, given a model that uses an alternative dual node approach, it is relatively straightforward to implement the numerically more correct consistent dual node approach.

Figure 1 .
Figure 1.(a) Common nodes in cell-centered schemes.(b) Dual nodes in cell-centered schemes.(c) Common nodes and colocated dual nodes in vertex-centered schemes.(d) Dual nodes in vertexcentered schemes (not colocated).The white squares and white circles represent surface and subsurface nodes, respectively.The solid and dashed lines represent the primary mesh and the dual mesh, respectively.The grey-shaded area is a topmost discrete volume associated with a topmost subsurface node.The black dot represents the centroid of this volume.The coupling length l as depicted in this figure applies to the consistent dual node approach.

Figure 2 .Figure 3 .
Figure 2. Outflow response for excess saturation on a hillslope (first scenario) using different vertical discretizations.

Figure 4 .Figure 5 .
Figure 4. Simulated values at the common nodes for excess saturation on a hillslope (first scenario) with a cell-centered scheme and z = 0.0125 m.(a) Water depths.(b) Pressure heads.Nodes are numbered 1-5 in the downslope direction.

Figure 6 .
Figure 6.Outflow response for excess infiltration on a hillslope (second scenario) using different vertical discretizations.

Figure 9 .
Figure 9. Simulated values for excess infiltration on a hillslope with a cell-centered scheme (second scenario) and z = 0.2 m.(a) Water depths at the surface nodes.(b) Pressure heads at the topmost subsurface nodes.Nodes are numbered 1-5 in the downslope direction.

Figure 10 .
Figure 10.Outflow response for flooding an unsaturated hillslope using different vertical discretizations.

Figure 13 .
Figure 13.Simulated values for excess infiltration (third scenario) on a hillslope with a cell-centered scheme and z = 0.2 m.(a) Water depths at the surface nodes.(b) Pressure heads at the topmost subsurface nodes.Nodes are numbered 1-5 in the downslope direction.
compute exchange fluxes between surface and topmost subsurface nodes, and therefore these schemes are technically dual node schemes.As explained in this study, it is crucial to observe that the schemes ofAn and Hydrol.Earth Syst.Sci., 21, 5709-5724, 2017 www.hydrol-earth-syst-sci.net/21/5709/2017/ R. de Rooij: Coupling surface-subsurface flow 5713 Yu (2014) and

Table 1 .
Abbreviations used in the figures.