Development and Application of a New Open-Source Integrated Surface–Subsurface Flow Model in Plain Farmland

: Accurately characterizing rainfall runoff processes in plain farmland, especially at the plot scale with significant micro-topographic features, has presented challenges. Integrated surface–subsurface flow models with high-precision surface flow modules are appropriate tools, yet open-source versions are rare. To address this gap, we proposed an open-source integrated surface–subsurface flow model called the FullSWOF-Plain model, in which the one-dimensional subsurface module Hydrus-1D was integrated with a modified two-dimensional surface water flow module (FullSWOF-2D) using the sequential head method. Various experimental scenarios were simulated to validate the model’s performance, including two outflow cases (i.e., 1D and 2D) without infiltration, a classical one-dimensional infiltration case, and two typical rainfall events at the experimental field. The results demonstrate the accuracy of this proposed model, with the Nash–Sutcliffe efficiency (NSE) of the simulated discharge exceeding 0.90 in the experimental field case and the root mean squared error (RMSE) values for soil moisture at five depths consistently below 0.03 cm 3 /cm 3 . However, we observed a lag in the simulated response time of soil moisture due to the neglect of preferential flow. The micro-topography significantly influenced ponding time and ponding areas. Lower local terrain normally experienced earlier surface ponding. Scattered surface ponding water first occurred in the ditch and followed in the relatively low areas in the main field. The concentration process of surface runoff exhibited hierarchical characteristics, with the drainage ditch contributing the most discharge initially, followed by the connection of scattered puddles in the main field, draining excess surface water to the ditch through rills. This quantitative study sheds light on the impact of micro-topography on surface runoff in plain farmland areas.


Introduction
It is vital to characterize rainfall runoff processes accurately in plain farmland areas, as these affect the local water balance greatly and cause soil erosion and rapid solute transport into rivers [1,2].Surface runoff in these areas occurs intermittently as overland flow [3].Micro-topography is the fundamental feature of plain farmland areas, because the slope can be ignored [4,5].Micro-topographic characteristics have a remarkable effect on surface connectivity [6], infiltration volume [7], depression storage [4,8,9], runoff generation, and concentration processes [5,[10][11][12].
However, it is a challenge to accurately characterize rainfall runoff processes in plain farmland areas, especially at the plot scale, with significant micro-topographic features.Fortunately, integrated surface-subsurface flow models with two-dimensional surface flow modules are appropriate and powerful tools in such situations, which can obtain the distribution of surface water depths and velocities as well as the outlet discharge.Since Freeze and Harlan [13] presented a blueprint for a distributed physically-based hydrological model, various physically-based intergrated surface-subsurface flow models have occurred in recent years [14][15][16][17][18][19][20], improving the quantitative characterization of nearsurface hydrologic response and runoff generation greatly [21,22].
Furthermore, another alternative integrated surface-subsurface flow modeling framework is to be noted, with a high accuracy of surface flow as well as high computational efficiency.In this framwork, fully dynamic SWEs are used for surface flow and runoff calculations, coupled with a conceptual infiltration model, such as the Green-Ampt model [46][47][48] or its revised version [49,50], the Horton model [51] and the SCS-CN method [52].However, such models cannot simulate the variations of layered soil moisture content, as well as the groundwater table.
Therefore, a distributed integrated surface-subsurface flow model, with 2D fully dynamic SWEs in the surface flow module and a 1D Richards' equation in the subsurface flow module, is still needed for accurately characterizing rainfall runoff processes in plain farmland areas.In recent years, many excellent open-source surface flow, as well as subsurface flow, software packages have been published.For example, FullSWOF-2D is a very famous surface flow program [50,53] that solves fully dynamic SWEs on a structured mesh using the finite volume method [54].The key module of FullSWOF-2D has been widely used in many models [54][55][56][57], demonstrating excellent portability.Similarly, the open-source Hydrus-1D, basen on the 1D Richards' equation, is known for its stability and effectiveness in simulating subsurface flow [58,59] and has been widely used [60].It seems worthwhile and interesting to integrate existing open-source software packages to develop an integrated surface-subsurface flow model for plain farmland.
In this study, we proposed an open-source integrated surface-subsurface flow model, the FullSWOF-Plain model, in which the one-dimensional subsurface module Hydrus-1D was integrated with a modified two-dimensional surface water flow module (FullSWOF-2D).This proposed model could acquire high-accuracy simulated surface water depths, as well as layered soil moisture content, via applying 2D SWEs and the 1D Richards' equation in the surface and subsurface flow module, respectively, balancing computational accuracy and efficiency.This paper aims to (1) propose an integrated surface-subsurface flow model and (2) evaluate its performance in a plain farmland area.We introduce the situation of the study area and the proposed model structure in Section 2. Section 3 presents the modeling results of various experimental scenarios to validate the model's performance.We discuss some important points and provide an outlook in Section 4. Conclusions are outlined in the final section.

Study Area
We selected a typical plain farmland field (N 31 • 43 ′ 30 ′′ , E 119 • 28 ′ 15 ′′ , Figure 1b) in Hongqi polder, located in Jintan City in the Taihu Lake Basin (Figure 1a).The annual average rainfall and evaporation are 723-1836 mm and 1041-1562 mm, respectively.Oilseed rape was planted from autumn to the following spring, and soybeans were planted in summer.The annual average phreatic water table depth was about 0.6 m.The whole area of the field was 1008 m 2 .We collected 4186 scattered elevation sample points, with an average sampling interval of 0.5 m.The spline interpolation method was applied to obtain a 0.1 m grid-resolution DEM (Figure 1c).The elevations varied from 2.45 m to 3.74 m, with a mean of 3.14 m.The core area was surrounded by banks to obstruct water.A 1.5 m-depth waterproof barrier was inserted into the banks to prevent lateral leakage.The area inside the isolated boundary was about 912 m 2 , with an average length of 87 m and average width of 10.5 m.Meanwhile, both an inner ditch and outer ditch were excavated surrounding the main field.The inner ditch collected the outflow from the main field, and the outer ditch collected the lateral leakage flow from the rice fields outside.We found that the soil profile in the field exhibited distinct layered characteristics.Specifically, the top 0-13 cm layer was the plow layer (A), the 13-22 cm layer was the plow pan (A p ), the 22-47 cm layer was the percogenic horizon (P), and below 47 cm was the waterloggogenic horizon (W).A 90 • V-notch weir was installed at the outlet of the inner ditch to monitor discharge.Local rainfall was monitored using a TE525MM tipping bucket rain gauge, and five soil moisture probes (CS616) were buried at five depths, i.e., 10 cm, 20 cm, 40 cm, 60 cm, and 80 cm.A groundwater observation well was located at the center of the field.

Study Area
We selected a typical plain farmland field (N 31°43′30″, E 119°28′15″, Figure 1b) in Hongqi polder, located in Jintan City in the Taihu Lake Basin (Figure 1a).The annual average rainfall and evaporation are 723-1836 mm and 1041-1562 mm, respectively.Oilseed rape was planted from autumn to the following spring, and soybeans were planted in summer.The annual average phreatic water table depth was about 0.6 m.The whole area of the field was 1008 m 2 .We collected 4186 scattered elevation sample points, with an average sampling interval of 0.5 m.The spline interpolation method was applied to obtain a 0.1 m grid-resolution DEM (Figure 1c).The elevations varied from 2.45 m to 3.74 m, with a mean of 3.14 m.The core area was surrounded by banks to obstruct water.A 1.5 m-depth waterproof barrier was inserted into the banks to prevent lateral leakage.The area inside the isolated boundary was about 912 m 2 , with an average length of 87 m and average width of 10.5 m.Meanwhile, both an inner ditch and outer ditch were excavated surrounding the main field.The inner ditch collected the outflow from the main field, and the outer ditch collected the lateral leakage flow from the rice fields outside.We found that the soil profile in the field exhibited distinct layered characteristics.Specifically, the top 0-13 cm layer was the plow layer (A), the 13-22 cm layer was the plow pan (Ap), the 22-47 cm layer was the percogenic horizon (P), and below 47 cm was the waterloggogenic horizon (W).A 90° V-notch weir was installed at the outlet of the inner ditch to monitor discharge.Local rainfall was monitored using a TE525MM tipping bucket rain gauge, and five soil moisture probes (CS616) were buried at five depths, i.e., 10 cm, 20 cm, 40 cm, 60 cm, and 80 cm.A groundwater observation well was located at the center of the field.Some physical and hydraulic parameters of the five layers in three soil profiles were measured (Table 1).The soil texture above a 60 cm depth was silty clay loam, while it was silty clay below a 60 cm depth.The dry bulk densities (ρ b ) below a 20 cm depth were similar, exceeding the values of the surface soil.The average porosity (n) decreased with depth.The saturated hydraulic conductivities (K s ) of the five measured layers were all small, especially that of the 60-80 cm layer.The two-dimensional shallow water equations for FullSWOF-2D are stated as the following equations for each computational cell: where h (m) is the surface water depth, u (m/s) and v (m/s) are the cell depth average velocities in x and y directions, respectively.r (m/s) is the rainfall intensity, i (m/s) is the infiltration rate, g (m/s 2 ) is the gravity acceleration, and t (s) is time.S ox and S oy are functions of space z, given by S ox = −∂ x z(x,y) and S oy = −∂ y z(x,y), respectively.S fx and S fy are the cell friction slopes, obtained from Manning's formulas , respectively; n is Manning's roughness coefficient (s/m 1/3 ).The wall, Neumann, periodic, imposed depth, and imposed flow boundary conditions are all available in this module.A well-balanced scheme is adapted to guarantee the positivity of water depth and the preservation of steady states for specific hydrological features, such as during wet-dry transitions and a tiny water depth [50].
However, the original version of FullSWOF-2D was designed for rectangular areas [61].We modified the code to make it available for irregular regions.Specifically, a data structure in C++ was set up, which recorded the row number or column number, as well as the serial numbers of the first and last valid cell of this row or column (Figure 2).This data structure marked the boundaries in four directions (i.e., left, right, top, and bottom) in the study area, aiding in the setting of boundary conditions.

1D Subsurface Flow Module
One-dimensional vertical flow in Hydrus-1D is described by Richards' equation [62]: A schematic of marking boundary information with a data structure in an irregular study area.

1D Subsurface Flow Module
One-dimensional vertical flow in Hydrus-1D is described by Richards' equation [62]: where t is the time (min), θ is the volumetric water content (cm 3 /cm 3 ), h is the soil water pressure head (cm), z is the spatial coordinate (cm) defined as positive upward, K(h) is the unsaturated hydraulic conductivity function (cm/min), and S is a source or sink term (cm/min).Van Genuchten's θ − h and K − h relationships are used to determine soil hydraulic properties in this model [63]: where m = 1 − 1/n, S e = (θ − θ r )/(θ s − θ r ), and θ r and θ s denote the residual and saturated volumetric water content (cm 3 /cm 3 ), respectively; S e (-) is the relative saturation; α (cm −1 ), n (-), and m (-) are fitting parameters of the soil water characteristic curve.

Integrated Method
The integrated method used between surface water and subsurface water flow is the key part of an intergrated surface-subsurface flow model [64].To coupld the 2D surface flow interactively with the 1D surbsurface flow, we apply the sequential head method proposed by Morita et al. [26,27,65].The infiltration between two phases is used as the common internal boundary condition to be solved interactively, in which the potential infiltration capacity I p is calculated firstly.However, this variable, calculated using the original formula, shows a remarkable fluctuation during the receding process because of wet-dry transitions.In order to improve the continuity of the potential infiltration capacity, the infiltration is set as the Darcy flow between two vertically adjacent grids on the surface when there is surface ponding water, while the original formula is used when there is no ponding water.Therefore, the Darcy's formulas for calculating potential infiltration capacity can be represented as follows: where K s is the saturated hydraulic conductivity; K 0 and H 0 are the hydraulic conductivity and hydraulic head of the top node, respectively; K 1 and H 1 are the hydraulic conductivity and hydraulic head of the node below the top one, respectively; ∆z is the thickness of each discrete cell; D is the ponding water depth at a previous time step; and Y is the available water depth from rainfall and depth of flow at the surface, calculated as: where r is rainfall intensity and ∆t is the time increment.The water supply rate R s is calculated as: At each time step, we compare R s with I p to determine the actual infiltration rate i and then set the corresponding boundary condition for the subsurface flow module.Further details can be found in the original work.Moreover, the surface flow module is not updated until the maximum surface ponding water depth on the field is larger than 0.1 mm [66], which can improve computional efficiency.

Modeling Cases
We first validated the surface flow module.Kirstetter's artificial rainfall runoff experiment [67] and a classical V-sloping slope overland flow case [41] were selected as one-dimensional and two-dimensional cases, respectively.Secondly, the one-dimensional subsurface flow module was validated with a classical infiltration case in an unsaturated zone, as well as a typical rainfall event without surface runoff at the experimental field.Finally, a rainstorm event at the experimental field was simulated using the FullSWOF-Plain model.We evaluated the model performance using the Nash-Sutcliffe efficiency (NSE) [68] and root mean squared error (RMSE).

One-Dimensional Outflow Case
Kirstetter's artificial rainfall-runoff experiment [67] was selected as a one-dimensional outflow case.The test flume had a length of 4.04 m and width of 11.5 cm.The slope of the flume (S o ), as well as the artificial rainfall intensity (I rain ), could be adjusted manually.The parameters of the three experimental cases are shown in Table 2.The cell size was set as 0.05 m.The initial surface water depths and velocities were zero.The HLL2 flux and Manning's formula were applied in this case.The CFL (Courant-Friedrichs-Lewy) factor was 0.4, and the Manning's coefficient n was 0.025 s/m 1/3 .The minimum water depth h min was set as 1 × 10 −12 m.A fictive cell was added adjacent to the outflow cell to improve computational stability [69,70].The whole simulation duration was 1000 s, while the artificial rainfall lasted for the first 600 s.The results of this proposed model were consistent with those in Kirstetter's work [67] (Figure 3).was 1000 s, while the artificial rainfall lasted for the first 600 s.The results of this proposed model were consistent with those in Kirstetter's work [67] (Figure 3).

Two-Dimensional Outflow Case
A classical tilted V-catchment overland flow case was chosen as the two-dimensional outflow case that was first proposed by Giammarco et al. [41] and then widely used in many articles [25,31,33,[71][72][73].A 20 m-wide low-lying channel was sandwiched between two symmetrical 1000 m × 800 m slopes.The slopes of two symmetrical slopes and the channel were 5% and 2%, respectively.We set the elevation of the top grid of the channel as 0 m, and the elevation of the whole area is shown in Figure 4a.The Manning's coefficient n was 0.015 s/m 1/3 for the valley side and 0.15 s/m 1/3 for the channel.A constant spatial discretization of 20 m (∆x = ∆y) was used.The critical depth condition was applied at the outlet.The rainfall intensity above the whole area was 10.8 mm/h, lasting for 90 min, with a total simulation time of 180 min.The CFL factor was 0.  , ParFlow [37] and this model.

A Classical One-Dimensional Infiltration Case
A classical one-dimensional infiltration case was selected [74], which concerned vertical flow in the unsaturated zone.The length and width of the experimental soil column were both 50 cm, and the height was 200 cm.The bottom and top faces corresponded to the groundwater table and the soil surface, respectively (Figure 5a).The soil in this column was homogeneous, the saturated hydraulic conductivity s K was 10 cm/day, and the po- rosity n was 0.45.The constitutive relations of the soil were: The initial pressure head was 0 cm at the bottom, −90 cm at the soil surface, and −97 cm elsewhere.The rainfall intensity at the soil surface was 50 mm/day, lasting for 10 days.The time step was set as 0.1 day.The simulated profiles of pressure heads at different moments during infiltration were plotted in Figure 5b.Infiltration caused a gradual increase in pressure head values, starting from the top of the soil column.Generally, the numerical results obtained from this proposed model fitted reasonably well with UNSAT2 We took Giammarco's work [41] as a reference solution and compared the simulated hydrographs in the V-catchment of the three models (Figure 4b).The discharge predicted by our proposed model agreed well generally with those of Giammarco et al. [41] and ParFlow [37], especially during the receding limb.The outlet flux of our model was underestimated during the first 30 min, while the peak flow fitted better than that of ParFlow.

A Classical One-Dimensional Infiltration Case
A classical one-dimensional infiltration case was selected [74], which concerned vertical flow in the unsaturated zone.The length and width of the experimental soil column were both 50 cm, and the height was 200 cm.The bottom and top faces corresponded to the groundwater table and the soil surface, respectively (Figure 5a).The soil in this column was homogeneous, the saturated hydraulic conductivity K s was 10 cm/day, and the porosity n was 0.45.The constitutive relations of the soil were: where k rw is relative hydraulic conductivity coefficient, S wr is residual water saturation, and h a is air entry value.As S wr was 0.333 and h a was 0, S w = 1 + 0.667h/100 and k rw = (S w − 0.333)/0.667.

One-Dimensional Infiltration Case at the Experimental Field
The soil moisture observation profile (10 cm, 20 cm, 40 cm, 60 cm, and 80 cm) in the experimental field could be conceptualized as a vertical soil column.We selected a typical light rainfall event (20160402) without surface runoff as the one-dimensional infiltration case.The total rainfall was 35.3 mm, and the initial groundwater depth was 1.13 m.The height of this soil column was set as 150 cm in the model, according to the initial groundwater depth.The vertical grid resolution was 1 cm.According to the characteristics of the genetic horizon and the buried depths of the five soil moisture probes, we divided the generalized soil profile into five layers, i.e., 0-13 cm, 13-22 cm, 22-47 cm, 47-75 cm, and 75-150 cm.Based on the measured parameters of the five layers (Table 1) and the hydraulic parameters in the van Genuchten model of a similar soil type [75,76], the calibrated parameters of the five layers are shown as Table 3.We ignored evaporation during the rainfall period.The simulated soil moisture content at depths of 10 cm, 20 cm, 40 cm, 60 cm, and 80 cm were output (Figure 6).The simulated results fitted very well with the measured values.The RMSE values of soil moisture at all five depths were less than 0.03 cm 3 /cm 3 .However, the simulated response time of soil moisture above 40 cm depth lagged behind, which might be caused by abundant macropore pathways in the topsoil.During rainfall, a portion of the rainfall infiltrated rapidly into the subsoil through macropore-type preferential pathways.However, such a process was not considered in this subsurface module.We calculated the increment of soil moisture in each layer.Interestingly, more than The soil column was discretized into 40 cells, and the vertical cell size ∆z was 5 cm.The initial pressure head was 0 cm at the bottom, −90 cm at the soil surface, and −97 cm elsewhere.The rainfall intensity at the soil surface was 50 mm/day, lasting for 10 days.The time step was set as 0.1 day.The simulated profiles of pressure heads at different moments during infiltration were plotted in Figure 5b.Infiltration caused a gradual increase in pressure head values, starting from the top of the soil column.Generally, the numerical results obtained from this proposed model fitted reasonably well with UNSAT2 and Huyakorn's model [74].

One-Dimensional Infiltration Case at the Experimental Field
The soil moisture observation profile (10 cm, 20 cm, 40 cm, 60 cm, and 80 cm) in the experimental field could be conceptualized as a vertical soil column.We selected a typical light rainfall event (20160402) without surface runoff as the one-dimensional infiltration case.The total rainfall was 35.3 mm, and the initial groundwater depth was 1.13 m.The height of this soil column was set as 150 cm in the model, according to the initial groundwater depth.The vertical grid resolution was 1 cm.According to the characteristics of the genetic horizon and the buried depths of the five soil moisture probes, we divided the generalized soil profile into five layers, i.e., 0-13 cm, 13-22 cm, 22-47 cm, 47-75 cm, and 75-150 cm.Based on the measured parameters of the five layers (Table 1) and the hydraulic parameters in the van Genuchten model of a similar soil type [75,76], the calibrated parameters of the five layers are shown as Table 3.We ignored evaporation during the rainfall period.The simulated soil moisture content at depths of 10 cm, 20 cm, 40 cm, 60 cm, and 80 cm were output (Figure 6).The simulated results fitted very well with the measured values.The RMSE values of soil moisture at all five depths were less than 0.03 cm 3 /cm 3 .However, the simulated response time of soil moisture above 40 cm depth lagged behind, which might be caused by abundant macropore pathways in the topsoil.During rainfall, a portion of the rainfall infiltrated rapidly into the subsoil through macropore-type preferential pathways.However, such a process was not considered in this subsurface module.We calculated the increment of soil moisture in each layer.Interestingly, more than 80 percents of the infiltration water was stored in the upper 40 cm soil layer, and almost one half was stored in the upper 10 cm soil layer.It was reasonable to only consider vertical infiltration in the unsaturated zone at the experimental field.

Validation of FullSWOF-Plain Model at Experimental Field
We simulated the rainfall runoff processes of the 20160808 rainfall event, which was known as a typical short-duration rainstorm.The rainfall was 23.7 mm, while the duration of the rainfall was around 30 min.The observed surface runoff was about 4.77 mm.Although fine-resolution DEMs can capture small-scale terrain structural features, they also cause an excessive computational burden.To balance simulation accuracy and efficiency, we chose a 0.3 m-resolution DEM, in which the number of effective discrete surface cells was 10,159 and the elevations varied from 2.73 m to 3.49 m.Hofer et al. [77] selected a similar spatial resolution (0.25 m) in the Chicken Creek basin (440 × 130 m).The elevation of the bottom boundary in the subsurface module was set as 2.29 m.The vertical grid resolution was 1 cm, with 889,052 effective discrete subsurface cells.The surface roughness coefficient was set as 0.02 s/m 1/3 according to the range of bare land (0.018-0.06 s/m 1/3 ) [78].The time step of the subsurface flow module was 0.1 min, and a variable time step was applied in the surface flow module.The CFL factor was set as 0.8.The soil hydraulic parameters of the five layers in Table 3 were also applied in this case.Variables of h, u, v, and q were output every minute.The discharge was calculated from the water depth, using a water depth-discharge relation.
The simulated hydrograph and the soil moisture content at the five depths are shown in Figure 7.As the simulation period was too short to make the groundwater respond, we did not output the simulated groundwater depth.The discharge at the outlet reached a peak value at the 25th minute and then receded gradually.The simulated hydrograph agreed reasonably well with the measured hydrograph, while the simulated peak value was slightly lower than the measured one.The NSE and RMSE of the simulated water depth were 0.958 and 0.876 cm, respectively.The NSE and RMSE of the simulated dis-

Validation of FullSWOF-Plain Model at Experimental Field
We simulated the rainfall runoff processes of the 20160808 rainfall event, which was known as a typical short-duration rainstorm.The rainfall was 23.7 mm, while the duration of the rainfall was around 30 min.The observed surface runoff was about 4.77 mm.Although fine-resolution DEMs can capture small-scale terrain structural features, they also cause an excessive computational burden.To balance simulation accuracy and efficiency, we chose a 0.3 m-resolution DEM, in which the number of effective discrete surface cells was 10,159 and the elevations varied from 2.73 m to 3.49 m.Hofer et al. [77] selected a similar spatial resolution (0.25 m) in the Chicken Creek basin (440 × 130 m).The elevation of the bottom boundary in the subsurface module was set as 2.29 m.The vertical grid resolution was 1 cm, with 889,052 effective discrete subsurface cells.The surface roughness coefficient was set as 0.02 s/m 1/3 according to the range of bare land (0.018-0.06 s/m 1/3 ) [78].The time step of the subsurface flow module was 0.1 min, and a variable time step was applied in the surface flow module.The CFL factor was set as 0.8.The soil hydraulic parameters of the five layers in Table 3 were also applied in this case.Variables of h, u, v, and q were output every minute.The discharge was calculated from the water depth, using a water depth-discharge relation.
The simulated hydrograph and the soil moisture content at the five depths are shown in Figure 7.As the simulation period was too short to make the groundwater respond, we did not output the simulated groundwater depth.The discharge at the outlet reached a peak value at the 25th minute and then receded gradually.The simulated hydrograph agreed reasonably well with the measured hydrograph, while the simulated peak value was slightly lower than the measured one.The NSE and RMSE of the simulated water depth were 0.958 and 0.876 cm, respectively.The NSE and RMSE of the simulated discharge were 0.953 and 0.211 L/s, respectively.The simulated soil moisture content at the five depths fitted very well with the measured values.However, the simulated response time of soil moisture above an 80 cm depth clearly lagged, as expected.The deeper the soil depth was, the smaller the increment in soil moisture was.FullSWOF-Plain model was proven to be accurate and effective for simulating rainfall -runoff processes in the experimental field.
Water 2024, 16, x FOR PEER REVIEW 11 of 16 soil depth was, the smaller the increment in soil moisture was.FullSWOF-Plain model was proven to be accurate and effective for simulating rainfall -runoff processes in the experimental field.Furthermore, simulated ponding water depths at four typical moments, i.e., the 11th, 16th, 25th and 110th minutes, were output to figure out the characteristics of the runoff generation processes (Figure 8a).At the 11th minute, surface ponding water could only be seen in some scattered depressions in the ditch, with a tiny discharge at the outlet.Five minutes later (i.e., the 16th minute), the ponding water in the ditch increased rapidly.Meanwhile, ponding water could also clearly be found in many depressions in the main field.Subsequently, the discharge reached a peak value at the 25th minute.The surface ponding water depths at this moment were much deeper than those of the 16th minute.Remarkably, some adjacent scattered depressions in the main field connected, forming patches of puddles.These puddles also connected with the rills nearby, which connected with the surrounding ditch.Therefore, the excess surface water in the puddles drained into the ditch through the rills (Figure 8b).After the rainfall, the discharge receded gradually.At the 110th minute, the ponding water close to the outlet receded fastest, due to drainage.The ponding water of the puddles in the main field receded gradually due to infiltration, which made the connections of some puddles break (Figure 8c).Furthermore, simulated ponding water depths at four typical moments, i.e., the 11th, 16th, 25th and 110th minutes, were output to figure out the characteristics of the runoff generation processes (Figure 8a).At the 11th minute, surface ponding water could only be seen in some scattered depressions in the ditch, with a tiny discharge at the outlet.Five minutes later (i.e., the 16th minute), the ponding water in the ditch increased rapidly.Meanwhile, ponding water could also clearly be found in many depressions in the main field.Subsequently, the discharge reached a peak value at the 25th minute.The surface ponding water depths at this moment were much deeper than those of the 16th minute.Remarkably, some adjacent scattered depressions in the main field connected, forming patches of puddles.These puddles also connected with the rills nearby, which connected with the surrounding ditch.Therefore, the excess surface water in the puddles drained into the ditch through the rills (Figure 8b).After the rainfall, the discharge receded gradually.At the 110th minute, the ponding water close to the outlet receded fastest, due to drainage.The ponding water of the puddles in the main field receded gradually due to infiltration, which made the connections of some puddles break (Figure 8c).

Discussion
This proposed model can accurately simulate rainfall runoff processes at a typical plot-scale plain farmland area with significant micro-topographic features.The visualization of temporal surface water distribution is crucial for understanding the mechanism of runoff generation.It can be seen from the whole rainfall runoff processes in rainfall event 20160808 that the micro-topography at multiple scales affected the ponding time and ponding areas greatly.Specifically, scattered surface ponding water first occurred in the ditch, which was the lowest part in the field, corresponding to the smallest soil water deficit.Then, surface ponding water could be seen in relatively low areas in the main field, which was higher than the ditch.Generally, lower local terrain experienced earlier surface ponding.
The concentration process of surface runoff presented hierarchical characteristics [10,79].Obviously, the drainage ditch outflowed first, contributing the majority of the discharge during the early stage, which could be called the partial outflow period.As some adjacent scattered puddles in the main field connected, the excess surface water in the puddles drained into the ditch through some rills.Therefore, the initial partial outflow period turned into a whole-area outflow period.The main field contributed a large part of the outflow in the whole-area outflow period [80].Overall, the micro-topography indeed affected the early outflow stage, while it impacted very little on surface runoff once the whole region outflowed [81].
Nevertheless, this proposed model could only be applied for rainfall event simulation recently, ignoring lateral subsurface flow and evaporation.We expect to introduce water balance theory to obtain a stable regional groundwater table after a single rainfall event, to make up for the role of lateral subsurface flow.Moreover, the evaporation condition will also be added into the model, so that this proposed model can be suitable for long-term simulation.Although Hydrus-1D proved reasonable for simulating the infiltration process in the experimental field, the response time of the simulated soil moisture clearly lagged because of the model neglecting the role of preferential flow.Therefore, it is necessary to incorporate the preferential flow process into the subsurface module in the future [36].

Discussion
This proposed model can accurately simulate rainfall runoff processes at a typical plotscale plain farmland area with significant micro-topographic features.The visualization of temporal surface water distribution is crucial for understanding the mechanism of runoff generation.It can be seen from the whole rainfall runoff processes in rainfall event 20160808 that the micro-topography at multiple scales affected the ponding time and ponding areas greatly.Specifically, scattered surface ponding water first occurred in the ditch, which was the lowest part in the field, corresponding to the smallest soil water deficit.Then, surface ponding water could be seen in relatively low areas in the main field, which was higher than the ditch.Generally, lower local terrain experienced earlier surface ponding.
The concentration process of surface runoff presented hierarchical characteristics [10,79].Obviously, the drainage ditch outflowed first, contributing the majority of the discharge during the early stage, which could be called the partial outflow period.As some adjacent scattered puddles in the main field connected, the excess surface water in the puddles drained into the ditch through some rills.Therefore, the initial partial outflow period turned into a whole-area outflow period.The main field contributed a large part of the outflow in the wholearea outflow period [80].Overall, the micro-topography indeed affected the early outflow stage, while it impacted very little on surface runoff once the whole region outflowed [81].
Nevertheless, this proposed model could only be applied for rainfall event simulation recently, ignoring lateral subsurface flow and evaporation.We expect to introduce water balance theory to obtain a stable regional groundwater table after a single rainfall event, to make up for the role of lateral subsurface flow.Moreover, the evaporation condition will also be added into the model, so that this proposed model can be suitable for long-term simulation.Although Hydrus-1D proved reasonable for simulating the infiltration process in the experimental field, the response time of the simulated soil moisture clearly lagged because of the model neglecting the role of preferential flow.Therefore, it is necessary to incorporate the preferential flow process into the subsurface module in the future [36].
Even though our proposed model has been compared with ParFLOW in a twodimensional outflow case and with UNSAT2 in a one-dimensional infiltration case, some more validation work needs to be conducted in the future, such as comparing its performance with well-known integrated surface-subsurface flow models (e.g., PIHM [32] and WASH123 [29]) in the experimental field.Furthermore, this model needs fine-resolution meshes even when it is applied in plot-scale areas, which definitely constrains the computational efficiency.It is necessary to introduce a GPUs (graphics processing units)-accelerated modeling framework into the model [82] in future studies, which will make it possible to achieve long-term simulation under fine-resolution meshes.

Conclusions
A new open-source integrated surface-subsurface flow model, FullSWOF-Plain model, was proposed, in which the one-dimensional subsurface module Hydrus-1D was integrated with a modified two-dimensional surface water flow module FullSWOF-2D via the sequential head method.Various experimental scenarios were simulated to validate the performance of this model, including two outflow cases (i.e., 1D and 2D) without infilration, a classical one-dimensional infiltration case, and two typical rainfall events in the experimental field.This proposed model could acquire high-accuracy simulated surface water depths, as well as layered soil moisture content, via applying 2D SWEs and the 1D Richards' equation in the surface and subsurface flow modules, respectively.The key points are as follows: (1) The simulated results of this proposed model were accurate.The NSE of the simulated discharge exceeded 0.90 in the experimental field case, and the RMSE values for soil moisture at the five depths were consistently below 0.03 cm 3 /cm 3 .However, the simulated response time of soil moisture lagged due to the neglect of preferential flow.(2) Micro-topography at multiple scales greatly affected the ponding time and ponding areas.Lower local terrain normally experienced earlier surface ponding.Scattered surface ponding water initially occurred in the ditch and subsequently in the relatively low areas in the main field.(3) The concentration process of surface runoff exhibited hierarchical characteristics.The drainage ditch outflowed first, contributing the majority of the discharge during the early stage.As some adjacent scattered puddles in the main field connected, the excess surface water in these puddles drained into the ditch through rills.

Figure 1 .
Figure 1.Overview of study area: (a) location of Hongqi polder; (b) location of study area; (c) topography and observation instrument of study area.

Figure 1 .
Figure 1.Overview of study area: (a) location of Hongqi polder; (b) location of study area; (c) topography and observation instrument of study area.

Water 2024 , 16 Figure 2 .
Figure 2. A schematic of marking boundary information with a data structure in an irregular study area.
Water 2024, 16, x FOR PEER REVIEW 7 of 16

Figure 3 .
Figure 3. (a-c) Observed and simulated hydrographs of three experimental cases.

4 .
The simulated two-dimensional variables were output in the VTK (Visualization Tool Kit) data file format and then imported into ParaView software-4.3 (a visualization program for large 3-dimensional data sets) to export graphs.Water 2024, 16, x FOR PEER REVIEW 8

k.
is relative hydraulic conductivity coefficient, wr S is residual water satura- tion, and a h is air entry value.As wr S was 0.333 and a h was 0,The soil column was discretized into 40 cells, and the vertical cell size z Δ was 5 cm.

Water 2024 ,
16, x FOR PEER REVIEW 10 of 1680 percents of the infiltration water was stored in the upper 40 cm soil layer, and almost one half was stored in the upper 10 cm soil layer.It was reasonable to only consider vertical infiltration in the unsaturated zone at the experimental field.

Figure 6 .
Figure 6.Simulated results of layered soil moisture content in rainfall event 20160402 at experimental field.

Figure 6 .
Figure 6.Simulated results of layered soil moisture content in rainfall event 20160402 at experimental field.

Figure 7 .
Figure 7. Simulated results of experimental field in rainstorm event 20160808.

Figure 7 .
Figure 7. Simulated results of experimental field in rainstorm event 20160808.

Figure 8 .
Figure 8. Simulated results of rainstorm event 20160808: (a) simulated ponding water depths on field at four typical moments; (b) image of concentration process on field; (c) image of field at 110th minute.

Figure 8 .
Figure 8. Simulated results rainstorm event 20160808: (a) simulated ponding water depths on field at four typical moments; (b) image of concentration process on field; (c) image of field at 110th minute.

Table 1 .
The measured physical and hydraulic parameters of the five soil layers.

Table 2 .
Main parameters of three experimental cases.

Table 3 .
Calibrated soil hydraulic parameters in van Genuchten model of five layers in subsurface module.

Table 3 .
Calibrated soil hydraulic parameters in van Genuchten model of five layers in subsurface module.