A Dimension-Reduced Line Element Method for 3D Transient Free Surface Flow in Porous Media

: In order to reduce the numerical difﬁculty of the 3D transient free surface ﬂow problems in porous media, a line element method is proposed by dimension reduction. Different from the classical continuum-based methods, homogeneous permeable pores in the control volume are conceptualized by a 3D orthogonal network of tubes. To obtain the same hydraulic solution with the continuum model, the equivalent formulas of ﬂow velocity, continuity equation and transient free surface boundary are derivable from the principle of ﬂow balance. In the solution space of transient free surface ﬂow, the 3D problem is transformed into 1D condition, and then a ﬁnite element algorithm is simply deduced. The greatest advantage of the line element method is line integration instead of volume/surface integration, which has dramatically decreased the integration difﬁculty across the jump free surface. Through the analysis of transient free surface ﬂow in the unconﬁned aquifer, trapezoidal dam, sand ﬂume and wells, the transient free surface locations predicted from the proposed line element method generally agree well with the analytical, experimental and other numerical data in the available literatures, the numerical efﬁciency can also be well guaranteed. Furthermore, the hydraulic anisotropy has signiﬁcant effect on the evolution of free surface locations and the shape of depression cones in spatial. The line element method can be expanded to model the 3D unsaturated seepage ﬂow, two-phase ﬂow and thermos problems in porous media because of the similarity between the similarity of Darcy’s law, Buckingham Law and Fourier’s law.


Introduction
Transient free surface seepage problems are prevalent in the reservoir-dam system and have a significant effect on the landslide risk in the reservoir and dam stability [1][2][3][4][5].Due to the indeterminacy of the free surface and seepage boundary over time, the analytical solutions are limited to the homogeneous, isotropic media with simple geometry [6][7][8].Alternatively, the continuum-based numerical methods have been extensively used to solve the transient free surface seepage problems.
In the early stage, Neuman and Witherspoon [9] proposed a finite element method for unsteady free surface flow in porous media, but the element updates near the moving free surface were required in each timestep.In order to avoid the adaptive mesh updates, a residual flow method with fixed mesh was developed by Desai et al. [10] and verified by the consistency with the available closed-form solution and measured data in a trapezoidal dam.To solve the transient seepage problems with the complicated drainage system, Chen et al. [11] established a parabolic variational inequality method by transferring the moving free surface and seepage boundaries into natural boundaries.Its feasibility was Water 2023, 15, 3072 2 of 20 illustrated by the agreements with the in situ monitoring data of an underground powerhouse.To reduce the high cost of mesh generation in the finite element methods, Rafiezadeh and Ataie-Ashtiani [12] presented a boundary element method to model transient seepage flow in both isotropic and anisotropic media with only boundary mesh, whereas the freesurface nodes still needed to be repositioned in each timestep.Moreover, there are some other alternative methods for numerical solutions of transient free surface flow in porous media, such as the finite volume method [13,14], numerical manifold method [15] and so forth.In the above methods, both the pores and grains in the porous media were treated as a homogeneous continuum.In nature, the pores are the occurrence space and flow paths of groundwater rather than solid grains.However, it is quite difficult to capture the flow process of numerous pores.Thus, irregular pore flow is conceptualized as a 3D homogeneous flow through the whole domain of the porous media based on the representative elementary volume (REV, [16]).Alternatively, in the equivalent pipe network model [17][18][19][20] or line element model [21,22] the homogeneous flow through the whole porous media was replaced by 1D pipe/line flow, in which the high-dimensional seepage flow problems were reduced to a lower-dimensional form in solution space, and the computational cost was substantially decreased for memory storage and numerical iteration.Moreover, these approaches provided a method of dimensional reduction to deal with the seepage flow problems through the fractured rock [23,24].The equivalent pipe network model (PNM) or line element model (LEM) was powerful in analyzing the confined and unconfined, steady and transient seepage flow in 2D isotropic and anisotropic media.
To the best of authors' knowledge, previous LEMs are limited in the 2D porous media, and the 3D transient free surface flow in porous media has not been addressed yet.In reality, the 3D transient free surface flow is more practical than the 2D condition.Nevertheless, the surface integration across the jump free surface in the continuum-based methods is very complicated.Firstly the polygon shape and spatial plane of the free surface intersected with the tetrahedron/hexahedron elements need to be confirmed.Then, the contribution of flow exchange in the matrix with the specific yield is calculated by triangulating the polygon, which highly depends on the locations of free surface and Gaussian points of triangular subdivision [11].Therefore, the previous numerical methods including the continuumbased model and PNM/LEM are mainly concentrated in 2D flow conditions and generally mesh-dependent, which may lead to great mathematical difficulty and computational cost in dealing with dynamic free surface.In order to reduce the numerical difficulty and enhance the numerical stability in transient free surface flow in porous media, a line element method is established in this study through dimensional reduction, and the 3D transient free surface flow problem is transformed into a 1D numerical solution space, in which the mesh preparation and line integral with dynamic free surface become quite simple.Firstly, the governing equations and boundary conditions of the 3D transient free surface flow in the porous media are described in mathematics.Secondly, the line element model is proposed with dimension reduction, and a simple finite line element algorithm is developed.Thirdly, the proposed numerical model is verified through comparisons with the unconfined aquifer, trapezoidal dam, sand flume and well.

Governing Equations
When the water level of upstream reservoir suddenly descends from φ 1 to φ 2 , as shown in Figure 1, the free surface in the 3D soil dam will gradually decrease from S f 1 to S f 2 over time.In general, the averaged flow velocity through the control volume below the free surface is given by Darcy's law [16] as Water 2023, 15, 3072 3 of 20 where v is the velocity of seepage flow (LT −1 ), k is the hydraulic conductivity (LT −1 ), and φ is the potential head (L).Note that x, y and z denote the three principal axes of hydraulic conductivity.
where v is the velocity of seepage flow (LT −1 ), k is the hydraulic conductivity (LT −1 ), and φ is the potential head (L).Note that x, y and z denote the three principal axes of hydrau- lic conductivity.During the discharge process, the relationship between the net flow through the surface areas and the rate of change of water amount in the control volume can be given by the following continuity equation [10][11]: where Ss is the specific storage (L −1 ) as the water volume released from a unit soil volume per unit declines in potential head [16], which is determined by the compressibility of solid skeleton and water.Generally, the water release as a result of the specific storage can be ignored in comparison to that of the gravitational specific yield.Thus, Equation ( 2) is reduced to During the discharge process, the relationship between the net flow through the surface areas and the rate of change of water amount in the control volume can be given by the following continuity equation [10,11]: where S s is the specific storage (L −1 ) as the water volume released from a unit soil volume per unit declines in potential head [16], which is determined by the compressibility of solid skeleton and water.Generally, the water release as a result of the specific storage can be ignored in comparison to that of the gravitational specific yield.Thus, Equation ( 2) is reduced to ∂v x ∂x + ∂v y ∂y Although the above equation is identical to the continuity equation of the steady free surface flow [25], the distribution of the potential head φ subject to the transient free surface is a function of space and time, which is different from the solution of the steady flow as it is only dependent on space.Meanwhile, the transient flow process should also satisfy the initial condition and boundary conditions below: (1) The initial condition [11] φ(x, y, z, t)| t=0 = φ 0 (x, y, z) (2) The potential head boundary S 1 [11] where the overline represents prescribed value, and S 1 denotes the flow boundaries below the upstream and downstream water levels and are valued by the given water heads φ 1 (t) and φ 2 (t), respectively.(3) The flux boundary S 2 [11] where n denotes the unit normal vector.As shown in Figure 1a, the base of the soil dam is impervious, and its flux is equal to zero.(4) The Signorini-type seepage boundary S 3 Excluding the potential head and flux boundaries, the remaining potential seepage boundary can be divided into two sections: the upper part above the transient free surface S 3u with φ(t) ≤ z, q n (t) = 0 and the lower part below the transient free surface S 3l with φ(t) = z, q n (t) ≤ 0. Thus, the boundary condition on the potential seepage face can be unified as [11] φ(t) ≤ z, q n (t (5) The transient free surface boundary S f Besides the zero-pressure condition φ = z, the flow exchange through the transient free surface can be calculated as [11] q n (t) = −µ ∂φ ∂t cos θ (8) where µ is the specific yield, and θ is the angle between the outward normal line of the free surface and the upwards vertical line.

Equivalent Flow Velocity and Hydraulic Conductivity
The porous media in micro scale is composed of countless solid grains and pores, and the calculative scale is very great if the Navier-Stokes equations are applied on a huge number of pores.Instead, the porous media can be averaged as a continuum using the concept of REV, where the macro hydraulic property performs stably with various sample sizes.However, this assumption contradicts to the fact that the grains cannot conduct water, with the exception of the pores.Conversely to the traditional continuum model, the homogeneous porous medium can be seen as the combination of regular arrayed spherical grains, as shown in Figure 1b.Due to the similarity of flow property between the pore and tube [26][27][28], the permeable pores between the solid grains in the x, y and z directions are conceptualized by a 3D orthogonal network of tubes, as shown in Figure 1c.For simplicity, the tube is represented by the line element in this study.
To ensure solution equivalence between the line element model and continuum model, these two models should have identical distributions of potential head and hydraulic gradient in the 3D space.Therefore, the flow velocity through the line elements can be expressed in the form of Darcy's law: where v l and k l denote the flow velocity (LT −1 ) and hydraulic conductivity (LT −1 ) of the line elements.
Subsequently, the tube flow rates Q l (L 3 T −1 ) of the line elements yield where A l denotes the cross-section area of the line elements (L 2 ).Combing Equations ( 9) and ( 10), the total flow rates through the control volume, as shown in Figure 1c, yield: where N x , N y and N z are the numbers of line elements, which is obtained by: where x, y and z are the length, width and height of the control volume; and B x , B y and B z are the spacings.
Based on the continuum model, the total flow rates are obtained by combing Equation (1): Through the comparison of Equations ( 11) and ( 13), the equivalent hydraulic conductivities through the line elements yield If A lx, A ly and A lz are specified with a unit of area, Equation ( 14) can be simplified as

Equivalent Continuity Equations
As shown in Figure 1c, there is only horizontal velocity in the x-directional line elements, while the other two vertical velocities are vanished.Inserting Equation (9) into Equation (3) gives Similarly, the continuity equations in the y-directional and z-directional line elements can be derived as ∂v ly ∂y

Equivalent Transient Free Surface Boundary
As shown in Figure 2, the free surface declines continuously from S f 1 to S f 2 in the time interval ∆t, and the total amount of water draining out through the horizontal section is equal to u ∂φ ∂t dt • ∆x∆y.Based on Equation ( 8) and the concept of specific yield, only the vertical line elements can contribute to water release, while the horizontal line elements have cos 90 • = 0, which yields where µL is the equivalent specific yield.Submitting Equation (12) into Equation (19) gives Thus, the two-dimensional planar flow across the free surface has been transferred to the 1D vertical flow in the line elements as

Unified Formulations and Numerical Procedure
By substituting the 3D Cartesian coordinate, a 1D local coordinate system l is applied in the line element model.When the two endpoints of any element are marked with i and j to establish the 1D coordinate system l, i is treated as the origin, with the positive direction points to j.Hence, the flow velocities in Equation ( 1) can be generalized as: Equation ( 22) is generally reasonable for the completely saturated flow below the free surface, the transient flow zone is not fixed with the evolution of free surface, and the finite element mesh needs to be regenerated during the iterative process, which would lead to divergent solutions [29].In order to obtain the transient flow solutions in the whole domain, including the unsaturated domain above the free surface, a continuous penalized Submitting Equation (12) into Equation (19) gives Thus, the two-dimensional planar flow across the free surface has been transferred to the 1D vertical flow in the line elements as

Unified Formulations and Numerical Procedure
By substituting the 3D Cartesian coordinate, a 1D local coordinate system l is applied in the line element model.When the two endpoints of any element are marked with i and j to establish the 1D coordinate system l, i is treated as the origin, with the positive direction points to j.Hence, the flow velocities in Equation ( 1) can be generalized as: Water 2023, 15, 3072 7 of 20 Equation ( 22) is generally reasonable for the completely saturated flow below the free surface, the transient flow zone is not fixed with the evolution of free surface, and the finite element mesh needs to be regenerated during the iterative process, which would lead to divergent solutions [29].In order to obtain the transient flow solutions in the whole domain, including the unsaturated domain above the free surface, a continuous penalized Heaviside function H λ (φ − z) [11] is employed to unify the saturated and unsaturated flow in the entire domain, which can be expressed as where λ is a penalty parameter.Subsequently, the continuity Equations ( 16)-( 18) can also be unified as The 1D forms of the initial condition and boundary condition can be updated as The solution of potential head within a 1D form in the whole domain should also satisfy the continuity Equation (25), initial condition and boundary conditions, i.e., Equations ( 26)- (30).From the variation principle, the functional I(φ) of Equation ( 25) is given as Minimizing Equation (31) yields By applying the 1D linear shape function, the discrete form for Equation (32) can be expressed as where where η is the time step, and ∆t is the time interval.The convergence criterion during each time step is defined as: where m is the iteration number, and ε is the error tolerance, which is specified as 0.001 in the next example analysis.

An Unconfined Aquifer
The proposed method is firstly applied to predict the evolution of transient free surfaces in a semi-infinite isotropic unconfined aquifer, as shown in Figure 3.The water level at the left boundary suddenly descends from 10 m to 5 m, this transient seepage problem through the semi-infinite unconfined aquifer is generally governed by the 1D nonlinear Boussinesq equation [30], and the analytical solution in the x direction is obtained by Zhang [30]  Water 2023, 15, x FOR PEER REVIEW 9 of 22 where η is the time step, and t Δ is the time interval.
The convergence criterion during each time step is defined as: where m is the iteration number, and ε is the error tolerance, which is specified as 0.001 in the next example analysis.

An Unconfined Aquifer
The proposed method is firstly applied to predict the evolution of transient free surfaces in a semi-infinite isotropic unconfined aquifer, as shown in Figure 3.The water level at the left boundary suddenly descends from 10 m to 5 m, this transient seepage problem through the semi-infinite unconfined aquifer is generally governed by the 1D nonlinear Boussinesq equation [30], and the analytical solution in the x direction is obtained by Zhang [30]   Figure 4a presents the evolution of free surface locations over time from the proposed line element method, and good agreements with the analytical solutions can be commonly obtained for different mesh sizes.The variation of iteration steps over time plotted in Figure 4b shows that the numerical solutions in each timestep can be obtained with less than Figure 4a presents the evolution of free surface locations over time from the proposed line element method, and good agreements with the analytical solutions can be commonly obtained for different mesh sizes.The variation of iteration steps over time plotted in Figure 4b shows that the numerical solutions in each timestep can be obtained with less than three iteration steps, which indicates that the convergence of the proposed line element method is fast for different mesh sizes.Figure 5a,b shows the evolutions of free surface locations and iteration steps over time for different time steps, respectively.The mesh size is fixed as Bx = By = Bz = 1 m.The numerical results are almost consistent with the theoretical solutions and converge within three steps, which demonstrates that the proposed method is robust for different time steps.Free Surface(m) x(m)  Figure 5a,b shows the evolutions of free surface locations and iteration steps over time for different time steps, respectively.The mesh size is fixed as B x = B y = B z = 1 m.The numerical results are almost consistent with the theoretical solutions and converge within three steps, which demonstrates that the proposed method is robust for different time steps.

A Trapezoidal Dam
The laboratory test through a trapezoidal dam consists of the glass beads perf by Desai et al. [10] is shown in Figure 6.This isotropic dam is discretized by a 3D o onal network of line elements.The hydraulic parameters are given as kx = ky = kz = 0. and µ = 0.4.The water level at the left boundary gradually increases from 0 cm to 1

A Trapezoidal Dam
The laboratory test through a trapezoidal dam consists of the glass beads performed by Desai et al. [10] is shown in Figure 6.This isotropic dam is discretized by a 3D orthogonal network of line elements.The hydraulic parameters are given as k x = k y = k z = 0.1 cm/s and µ = 0.4.The water level at the left boundary gradually increases from 0 cm to 15 cm.Similar to the first example, the effect of mesh size and time step are shown in Figures 7 and 8.For different mesh sizes, the time step is fixed as 1 min.For different time step, the mesh size is set as 0.6 cm.The model results after 7 min and 9 min agree well with the experimental data conducted by Desai et al. [10] and the computational efficiency can be well guaranteed.For comparison, the simulated data attained from Desai et al. [10] are also given, and the deviation of free surface location at 9 min is much larger than the proposed line element method.Similar to the first example, the effect of mesh size and time step are shown in Figures 7 and 8.For different mesh sizes, the time step is fixed as 1 min.For different time step, the mesh size is set as 0.6 cm.The model results after 7 min and 9 min agree well with the experimental data conducted by Desai et al. [10] and the computational efficiency can be well guaranteed.For comparison, the simulated data attained from Desai et al. [10] are also given, and the deviation of free surface location at 9 min is much larger than the proposed line element method.Similar to the first example, the effect of mesh size and time step are shown in Figures 7 and 8.For different mesh sizes, the time step is fixed as 1 min.For different time step, the mesh size is set as 0.6 cm.The model results after 7 min and 9 min agree well with the experimental data conducted by Desai et al. [10] and the computational efficiency can be well guaranteed.For comparison, the simulated data attained from Desai et al. [10] are also given, and the deviation of free surface location at 9 min is much larger than the proposed line element method. (a)

A Sand Flume
The laboratory experiment through an isotropic sand flume reported by Hu et [31]is illustrated in Figure 9a, which contains five rectangular drainage tunnels.The draulic parameters of the fine sand soil are given as kx = ky = kz = 2.4 × 10 −3 cm/s and µ = 0.The water level at the left boundary gradually decreases from 1.0 m to 0.2 m after 48 m while the right water level remains at 0.2 m.During numerical analysis, the drainage tu nel is specified as a Signorini-type seepage boundary.The comparisons of free surf locations after 18 min, 30 min and 48 min between the proposed line element and obs vations from Hu et al. [31] are presented in Figure 9b.The calculated free surface locatio are very close to the measurements and are generally achieved within 6 iteration steps

A Sand Flume
The laboratory experiment through an isotropic sand flume reported by Hu et al. [31] is illustrated in Figure 9a, which contains five rectangular drainage tunnels.The hydraulic parameters of the fine sand soil are given as k x = k y = k z = 2.4 × 10 −3 cm/s and µ = 0.05.The water level at the left boundary gradually decreases from 1.0 m to 0.2 m after 48 min, while the right water level remains at 0.2 m.During numerical analysis, the drainage tunnel is specified as a Signorini-type seepage boundary.The comparisons of free surface locations after 18 min, 30 min and 48 min between the proposed line element and observations from Hu et al. [31] are presented in Figure 9b.The calculated free locations are very close to the measurements and are generally achieved within 6 iteration steps.

A 3D Well
Different from the previous examples with 2D condition, a 3D cylindrical well model, as shown in Figure 10a, is selected to assess the proposed method and investigate the transient free surface flow behaviors in the anisotropic porous media.The hydraulic parameters of the isotropic sand are given as k x = k y = k z = 9.57 inch/min and µ = 0.3.The water level in the well decreases gradually from 48 inch to 12 inch at the initial time, while the water level of the peripheral boundary of the tank remains at 48 inch.Different from the uniform-sized element mesh in the above examples, the well model is divided into three domains meshed with non-uniformed sizes (coarse mesh 4 m, medium mesh 2 m and fine mesh 1 m) as shown in Figure 10a, in which the well with dense mesh is specified as a Signorini-type seepage boundary.domains meshed with non-uniformed sizes (coarse mesh 4 m, medium mesh 2 m and fine mesh 1 m) as shown in Figure 10a, in which the well with dense mesh is specified as a Signorini-type seepage boundary.
Figure 10b shows the steady free surface location after 30 min.For comparison, the experimental data from Hall [32], numerical predictions from Cooley [33] and Clement et al. [34] are also plotted.The predicted free surface is very close to the measured and other numerical results.Besides the above isotropic case, two anisotropic cases are additionally applied to illustrate the influence of hydraulic anisotropy on the free surface variation [35][36].The anisotropic Ⅰ and Ⅱ cases are specified as kx = 10ky = kz = 9.57 in/min and kx = ky = 10kz = 9.57 in/min, respectively.The comparisons of free surface variation with time from the isotropic and anisotropic Ⅰ and Ⅱ cases are shown in Figure 11.The corresponding 3D depression cones of the three cases after 5 min are presented in Figure 12.
In Figure 11b, the free surfaces almost keep horizontal around the external boundary but decline rapidly near the well with time.In Figure 11c, the variation of free surfaces are distributed evenly.Compared to the isotropic case (Figure 11a), the free surface locations at the profile x = 0 (Figure 11b) of anisotropic Ⅰ turn to be much higher, while those at the profile y = 0 (Figure 11c) become much lower.This is because the decrease in ky can lead to the increase in flow resistance and the slight hydraulic gradient in the y direction.Conversely, the flow paths in the x direction are more unimpeded in terms of seepage flow, and the hydraulic gradients are more substantial.Therefore, the corresponding 3D de- Figure 10b shows the steady free surface location after 30 min.For comparison, the experimental data from Hall [32], numerical predictions from Cooley [33] and Clement et al. [34] are also plotted.The predicted free surface is very close to the measured and other numerical results.
Besides the above isotropic case, two anisotropic cases are additionally applied to illustrate the influence of hydraulic anisotropy on the free surface variation [35,36]

Conclusions
A line element method stimulated by dimension reduction is developed to model the 3D transient free surface flow in porous media, where the traditional homogeneous continuum is displaced by a 3D orthogonal network of line elements in the scale of representative elementary volume.Based on the principle of flow balance, the equivalent equations of flow velocity, continuity equation and boundary conditions are deduced to 12. Three-dimensional depression cone for isotropic and anisotropic cases.
In Figure 11b, the free surfaces almost keep horizontal around the external boundary but decline rapidly near the well with time.In Figure 11c, the variation of free surfaces are distributed evenly.Compared to the isotropic case (Figure 11a), the free surface locations at the profile x = 0 (Figure 11b) of anisotropic I turn to be much higher, while those at the profile y = 0 (Figure 11c) become much lower.This is because the decrease in k y can lead to the increase in flow resistance and the slight hydraulic gradient in the y direction.Conversely, the flow paths in the x direction are more unimpeded in terms of seepage flow, and the hydraulic gradients are more substantial.Therefore, the corresponding 3D depression cone of anisotropic I is elliptical with different curvature at the x and y directions.
In Figure 11d, the evolution of free surfaces are not evident with time, except for the seepage points near the well.The elevations of the free surfaces are much higher than those of the isotropic case, which is in accord with the observations by Wei et al. [24].This is because the decrease in k z leads to the increase in flow resistance in the z direction, while the seepage flow predominantly occurs in the horizontal plane, and the hydraulic gradient becomes mild.Similar to the isotropic case, the depression cone in the anisotropic II case is also circular because of the permeability symmetry k x = k y but with smaller curvature.Note that this special effect cannot be characterized by the two-dimensional seepage analysis [37].

Conclusions
A line element method stimulated by dimension reduction is developed to model the 3D transient free surface flow in porous media, where the traditional homogeneous continuum is displaced by a 3D orthogonal network of line elements in the scale of representative elementary volume.Based on the principle of flow balance, the equivalent equations of flow velocity, continuity equation and boundary conditions are deduced to guarantee identical solution with the continuum model.Then, the finite line-element algorithm is established, in which a local 1D coordinate system is applied instead of a 3D Cartesian coordinate, and the continuous penalized Heaviside function is used to describe the whole-domain flow, including the unsaturated flow above the free surface.It is found that there is no additional calculation of hydraulic parameter and mesh reconstruction during numerical analysis.Compared to the homogeneous continuum model, the 3D transient seepage problem is simplified into a 1D system in solution space to reduce the computational cost and numerical difficulty.Through the transient free surface flow analysis in four typical examples, the numerical predictions can be achieved in few iteration steps and agree well with the analytical, experimental and other numerical results in the current literature.In addition, the hydraulic anisotropy has significant effect on the evolution of free surface locations and the shape of depression cones in spatial.
The significant advantages of the dimensional reduction based on the line element model for the transient free surface flow can be concluded as the following: (1) Easy mesh preparation.The 3D flow domain can be easily meshed into three orthogonal sets of line elements rather than the tetrahedron and hexahedron elements, and each node is connected to a maximum of six adjacent nodes (twenty-six adjacent nodes for hexahedron elements).Thus, the computational cost of mesh generation and memory storage of the nonzero elements in the matrix are greatly reduced by the dimensional reduction in the line element model.(2) Simple line integral.The line integral is much simpler than the area or volume integral through the tetrahedron and hexahedron elements, with the latter depending on the Gaussian integral points.Especially for the surface integral term ∑ S f N T u f Ndxdy across the jump free surface in the continuum model, the polygon shape of the free surface across the tetrahedron/hexahedron elements and the spatial plane equation based on the isoparametric element needs to be confirmed firstly.Then, the triangular subdivision is employed to calculate the contribution of flow exchange in the matrix with the specific yield.However, in the matrix F of the line element model, there is no complicated area integral and no confirmation of the surface shape; only the position of the intersection between the line element and free surface is required to calculate the value of N T u f N. In addition, the element mesh around the free surface is fixed without updating the timestep and iteration number.The proposed method can be extended to model the 3D unsaturated seepage flow, two-phase flow and thermos problems in porous media because of the similarity between the similarity of Darcy's law, Buckingham Law and Fourier's law.

1 SFigure 1 .
Figure 1.(a) Transient free surface flow in a 3D soil dam.(b) The homogeneous media with arrayed spherical grains.(c) The 3D orthogonal network of tubes.

Figure 1 .
Figure 1.(a) Transient free surface flow in a 3D soil dam.(b) The homogeneous media with arrayed spherical grains.(c) The 3D orthogonal network of tubes.

Figure 2 .
Figure 2. The equivalent specific yield of the jump free surface.

Figure 2 .
Figure 2. The equivalent specific yield of the jump free surface.
based on the Laplace transform.The hydraulic parameters are specified as k x = k y = k z = 5 m/d and µ = 0.1.To evaluate the influence of mesh size and time step on the numerical accuracy and convergence, this aquifer is meshed by a 3D orthogonal line element with different mesh sizes (B x = B y = B z = 4 m, 2 m, 1 m) and simulated with various time steps (∆t = 5 d, 1 d, 0.1 d).
based on the Laplace transform.The hydraulic parameters are specified as kx = ky = kz = 5 m/d and µ = 0.1.To evaluate the influence of mesh size and time step on the numerical accuracy and convergence, this aquifer is meshed by a 3D orthogonal line element with different mesh sizes (Bx = By = Bz = 4 m, 2 m, 1 m) and simulated with various time steps ( t Δ = 5 d, 1 d, 0.1 d).

Figure 3 .
Figure 3.The orthogonal line element mesh.The dimensions of this aquifer are 500 m in length, 12 m in height and 8 m in width.Different from the solid element (tetrahedron and hexahedron), for which that the flow domain is mapped using the nodes and lines.

Figure 3 .
Figure 3.The orthogonal line element mesh.The dimensions of this aquifer are 500 m in length, 12 m in height and 8 m in width.Different from the solid element (tetrahedron and hexahedron), for which that the flow domain is mapped using the nodes and lines.

Figure 4 .
Figure 4. (a) Locations of transient free surface versus time.(b) Iteration steps versus time for different mesh sizes (B x = B y = B z = 4 m, 2 m, 1 m).

Figure 5 .
Figure 5. (a) Locations of transient free surface versus time.(b) Iteration steps versus time for v time steps ( t Δ = 5 d, 1 d, 0.1 d).The mesh size is fixed as Bx = By = Bz = 1 m.

Figure 5 .
Figure 5. (a) Locations of transient free surface versus time.(b) Iteration steps versus time for various time steps (∆t = 5 d, 1 d, 0.1 d).The mesh size is fixed as B x = B y = B z = 1 m.

Figure 6 .
Figure 6.The orthogonal line element mesh.The dimensions of the trapezoidal dam are 34.0 cm in length, 26.2 cm in height and 2.4 cm in width.

Figure 6 .
Figure 6.The orthogonal line element mesh.The dimensions of the trapezoidal dam are 34.0 cm in length, 26.2 cm in height and 2.4 cm in width.

Figure 6 .
Figure 6.The orthogonal line element mesh.The dimensions of the trapezoidal dam are 34.0 cm in length, 26.2 cm in height and 2.4 cm in width.

Figure 7 .Figure 7 .
Figure 7. (a) Locations of transient free surface versus time.(b) Iteration steps versus time for various mesh sizes (B x = B y = B z = 1.2 cm, 0.6 cm, 0.3 cm).The time step is fixed as 1 min.

Figure 8 .
Figure 8.(a) Locations of transient free surface versus time.(b) Iteration steps versus time for various time steps (∆t = 1 min, 0.5 min, 0.25 min).The mesh size is set as 0.6 cm.

Figure 9 .
Figure 9. (a) The orthogonal line element mesh.(b) Locations of transient free surface versus time.The dimensions of this flume are 1 m in length, 1.2 m in height and 0.1 m in width.The mesh size of 0.02 m and the time step of 0.5 min are used in this sand flume model.

Figure 9 .
Figure 9. (a) The orthogonal line element mesh.(b) Locations of transient free surface versus time.The dimensions of this flume are 1 m in length, 1.2 m in height and 0.1 m in width.The mesh size of 0.02 m and the time step of 0.5 min are used in this sand flume model.

Figure 10 .
Figure 10.(a) The orthogonal line element mesh.(b) Locations of steady free surface.The well is 48 inch in height, and the radii of the well and total sand tank are 4.8 inch and 76.8 inch.The well model is divided into three domains meshed with non-uniformed sizes (coarse mesh 4 m, medium mesh 2 m, fine mesh 1 m).

Figure 10 .
Figure 10.(a) The orthogonal line element mesh.(b) Locations of steady free surface.The well is 48 inch in height, and the radii of the well and total sand tank are 4.8 inch and 76.8 inch.The well model is divided into three domains meshed with non-uniformed sizes (coarse mesh 4 m, medium mesh 2 m, fine mesh 1 m).

Figure 11 .
Figure10bshows the steady free surface location after 30 min.For comparison, the experimental data from Hall[32], numerical predictions from Cooley[33] and Clement et al.[34] are also plotted.The predicted free surface is very close to the measured and other numerical results.Besides the above isotropic case, two anisotropic cases are additionally applied to illustrate the influence of hydraulic anisotropy on the free surface variation[35,36].The anisotropic I and II cases are specified as k x = 10k y = k z = 9.57 in/min and k x = k y = 10k z = 9.57 in/min, respectively.The comparisons of free surface variation with time from the

Author
Contributions: Writing-original draft, Y.C., Q.Y., Z.Y. and Z.P.All authors have read and agreed to the published version of the manuscript.