Numerical Solution of the Two-Dimensional Richards Equation Using Alternate Splitting Methods for Dimensional Decomposition

: Research on seepage ﬂow in the vadose zone has largely been driven by engineering and environmental problems a ﬀ ecting many ﬁelds of geotechnics, hydrology, and agricultural science. Mathematical modeling of the subsurface ﬂow under unsaturated conditions is an essential part of water resource management and planning. In order to determine such subsurface ﬂow, the two-dimensional (2D) Richards equation can be used. However, the computation process is often hampered by a high spatial resolution and long simulation period as well as the non-linearity of the equation. A new highly e ﬃ cient and accurate method for solving the 2D Richards equation has been proposed in the paper. The developed algorithm is based on dimensional splitting, the result of which means that 1D equations can be solved more e ﬃ ciently than as is the case with unsplit 2D algorithms. Moreover, such a splitting approach allows any algorithm to be used for space as well as time approximation, which in turn increases the accuracy of the numerical solution. The robustness and advantages of the proposed algorithms have been proven by two numerical tests representing typical engineering problems and performed for typical properties of soil.


Introduction
The subsurface flow under unsaturated conditions can be described by the Richards equation, whose numerical solution in the case of two-dimensional (2D) or three-dimensional (3D) problems is still not a trivial task [1]. This problem becomes especially significant when a large computational domain with a high spatial resolution and long simulation period is required. The difficulties are also enhanced by the strong non-linearity of the constitutive relations describing the hydraulic properties of porous media [2].
For solving the 2D and 3D Richards equations, a variety of numerical algorithms based on the finite difference (FD) [3][4][5], finite element (FE) [6][7][8], and finite volume (FV) methods [9][10][11][12] have been proposed. Among the aforementioned methods applied for solving multi-dimensional flow problems, it is necessary also to distinguish coupled algorithms, such as the control-volume FD scheme [13], the mixed FE scheme [14][15][16], and the approach based on semi-discretization, in the form of the method of lines (MoL) and transversal method of lines (TMoL) [17]. An interesting group of algorithms used to solve the Richards equation also provide the methods using the Kirchhoff transformation. Such approach for variably saturated 2D flow in layered soil has been applied by Suk and Park [18] as well as by Berninger et al. [19].
The methods based on FD are especially attractive due to their simplicity for implementation in numerical code, however, these methods are mainly limited to a regular domain. Thus, for complex 2D and 3D geometries with curvilinear boundaries, the FE and FV methods are particularly suitable with regard to their flexibility for spatial discretization. However, it is worth adding that algorithms derived for an unstructured and curvilinear grid using finite elements or volumes are often less accurate in comparison to a structured one. To eliminate the above disadvantage, the original equation posed on an irregular grid can be transformed into a problem which is solved using a computational rectangular grid. Such an approach for solving the 2D Richards equation has been proposed by Arrarás et al. [20] as well as by Deng and Wang [21].
The nature of the flow in the unsaturated zone states that while solving a given problem, a significant resolution of the spatial grid is often required. This situation is further complicated when the multi-dimensional Richards equation is considered. Then, the resulting system of algebraic equations contains a coefficient matrix with a very wide bandwidth of non-zero elements. The form of this bandwidth mainly depends on the structure of the assumed grid as well as on the numbering of nodes. It is worth adding that the bandwidth is the main factor that determines the efficiency of the numerical algorithm, in such a way that increasing the bandwidth involves a longer computational time.
In order to improve computational efficiency, the dimensional splitting method can be applied. Due to the splitting algorithm, the solution of a 2D or 3D equation at a given time step is reduced to the solution of a set of 1D equations for each of the directions of the coordinate system. Consequently, the obtained system of algebraic equations has a tri-diagonal coefficient matrix, which can be solved more efficiently than a system with a very wide bandwidth. It is worth noting that a tri-diagonal matrix can also be obtained using the method in which the 2D Richards equation is solved by the MoL in its x-direction and the TMoL in the y-direction [17].
For solving a multi-dimensional Richards equation, the dimensional splitting in the form of an alternating direction implicit (ADI) algorithm was applied by Weeks et al. [22] for unsaturated flow through inclined porous media, and by An et al. [23] for flow into layered soil. As it has been demonstrated by these authors that the solution obtained for the standard ADI method often does not converge even for a small time step. In particular, this is visible for 3D problems and partially saturated conditions [23]. For this reason, modifications of the standard ADI algorithm are required to ensure a stable solution of the Richards equation. Moreover, the family of ADI algorithms is inflexible because such splitting based on a strictly prescribed spatial and temporal approximation of each 1D equation, and consequently their applicability is limited.
To eliminate the above inconveniences during dimensional decomposition, the fractional step method, also called the Godunov method, can be applied. A fundamental description of this approach is provided by Toro [24] and LeVeque [25] for the 2D shallow water equations, whereas an application for solving the 2D diffusive wave equation is given by Gąsiorowski [26]. As a result of the splitting procedure, the obtained 1D equation for each direction can be solved by any algorithm based on the FD, FE, or VF methods. Moreover, the Godunov method is one of the most effective methods and usually leads to a solution with high accuracy when a linear problem is considered. Unfortunately, the Richards equation is highly non-linear. Consequently, dimensional splitting introduces an additional error into the numerical solution. As a result of the calculation over one time step, the obtained solutions differ from each other as well as from the corresponding unsplit or exact solution. This splitting error depends mainly on the time step and on the direction in which the calculations are started. It is worth adding that the splitting error is independent from the accuracy of the scheme applied for solving 1D equations [25]. Due to the splitting error, the Godunov method is only 1st order accurate with regard to time. In order to increase the accuracy of the splitting process, the Strang method can be used [27]. This method is 2nd order accurate, but requires about twice the amount of computational effort in comparison to the Godunov method, which significantly reduces the efficiency of the decomposition process.
In this study, in order to overcome the abovementioned flaws, we applied the modified Strang method as well as deriving an alternate splitting algorithm. The modification is performed by appropriately changing the sequence for a time integration. As a result, the efficiency of the modified method is very similar to that obtained with the Godunov method and at the same time the adequate accuracy of the numerical solution is ensured.

Richards Equation
In hydrological application, the Richards equation is often considered in a mixed form, which takes into account the water content and the water pressure head [11]. Such a form of equation is mass-conservative and for the 2D flow it can be written as follows: where θ(h) is volumetric water content, (-); h is pressure head (negative for the unsaturated zone), (m); t is time, (s); x, z are spatial variables, (m); and K x (h) and K z (h) are unsaturated hydraulic conductivity in x and z direction, respectively, (m·s −1 ). The first and second terms on the right side of Equation (1) are associated with diffusion transport due to capillarity, while the third term represents advective transport resulting from gravity. The Richards equation must be completed with the relationship between the pressure head, water content, and hydraulic conductivity, describing the hydraulic properties of the porous media. For this purpose, the Mualem-van Genuchten functions [28,29] can be applied: where S e is effective saturation, (-); 0 ≤ S e ≤ 1, K(h) is hydraulic conductivity, (m·s −1 ); C(h) is specific capacity, (m −1 ); K s is saturated hydraulic conductivity, (m·s −1 ); θ s and θ r are saturated and residual water content, respectively, (-); α is the parameter related to the average pore size, (m −1 ); and n and m are dimensionless parameters related to pore size distribution, m = 1 − 1/n, (-). The hydraulic conductivity can be also described using the exponential function in the form of Gardner relation as follows [18,30]: where h g is a scaling parameter. According to the assumed Gardner relation, the effective saturation and specific capacity take the following forms [18]: For such constitutive functions (Equations (5)-(7)), an exact solution of the Richards equation has been derived for the case of 2D as well as 3D [31].

Splitting Methods for Dimensional Decomposition
In order to perform the dimensional splitting, the Richards equation (Equation (1)) can be rewritten to the following operator form: where the differential operators A and B correspond to the terms in z and x directions, respectively. Assuming flow in an isotropic soil, these operators are expressed as: Equation (8) can be considered as a general form of 2D Richards equation, where the constitutive function describing the hydraulic conductivity and other relations take the appropriate form depending on the assumed model, i.e., the Mualem-van Genuchten, the exponential or the Brooks-Corey model [32].

Godunov Method
At first, let us consider the simplest splitting algorithm, in the form of the Godunov method. To this order, the spatial domain is overlaid with a rectangular mesh with N rows and M columns giving the total of M × N nodes ( Figure 1). According to the Godunov approach, over each time step the Richards equation (Equation (8)) can be rewritten as an equivalent sequence of the two following equations: where (1) and (2) are superscripts indicating the stage in directions z and x, respectively; ∆t is the time step; i and j are node indices with regard to directions z and x, respectively; and M and N are the number of nodes with regard to directions z and x, respectively. As a result, the splitting process is carried out in two stages using 1D equations. In the first stage, Equation (11) is solved in the z direction over the time step ∆t with the initial condition h (1) (x,z,t) = h 0 . The procedure is performed for the consecutive rows j = 1, 2, . . . , N (see Figure 1a). Then, using the obtained results as the initial condition h (2) (x,z,t) = h (1) (x,z,t + ∆t) in the second stage, Equation (12) is solved for x direction also over the time step ∆t. In this case, the calculations are performed with the consecutive columns i = 1, 2, . . . , M (Figure 1b). After the last stage, the values of the pressure head at the time level t + ∆t are obtained, i.e., h(x,z,t + ∆t) =h (2) (x,z,t + ∆t). Since the Richards equation is non-linear, the dimensional decomposition applied to its solution introduces a splitting error. As a result of the accuracy analysis, the Godunov method is only 1st order accurate with respect to the time. Consequently, in this case, the splitting error takes the following form [25]: where O(∆t 3 ) denotes the higher terms omitted in the Taylor series expansion. As can be seen, the splitting error depends mainly on the value of the time step ∆t as well as on the relation between the differential operators A and B. When the operator multiplication, i.e., product AB is equal to BA, then the applied method leads to a 2nd order accurate solution. This case takes place for the linear problem. However, with regard to the high non-linearity of the Richards equation, the operator multiplications differ from each other (AB BA). As a consequence, this introduces the splitting error, which decreases an accuracy of the final solution.

Standard and Modified Strang Method
In order to eliminate or to minimize the splitting error, the Strang method can be applied for the dimensional splitting [27]. Using this method, the decomposition process is carried out in three stages, such that the 2D Richards equation (Equation (8)) is decomposed into the following set of 1D equations: Similarly, as it was presented for the Godunov method, the corresponding calculations are performed using the grid presented in Figure 1. At the first stage, the subproblem in z direction (Equation (14)) is solved over a half time step of the length ∆t/2 with the initial condition h(x,z,t) = h 0 . The solution of this stage h (1) (x,z,t + ∆t/2) is used as the initial condition at the second stage for the solution of Equation (15) in x direction over a full time step ∆t. Finally, at the third stage, Equation (16) is integrated for another half time interval ∆t/2 again in z direction with the initial condition h (3) (x,z,t) = h (2) (x,z,t + ∆t). After the last stage at the time level t + ∆t, the solution h(x,z,t + ∆t) = h (3) (x,z,t + ∆t) is obtained. A schematic sketch of the Strang algorithm is shown in Figure 2a. The Strang method is 2nd order accurate with respect to the exact solution. Consequently, this method ensures a solution with a smaller splitting error than the Godunov method, even if the operator multiplications (AB or BA) differ from each other. However, while increasing accuracy up to the 2nd order by using the Strang method, at the same time, the computational efficiency of this method is significantly reduced in comparison with the sequential Godunov method (Equations (11)-(12)). As it is easy to check, the Strang method requires approximately 50% more calculations (with stages S = 3N T , where N T is the number of time steps) than the 1st order Godunov method (with stages S = 2N T ). Considering a reduction of computational work, Strang [27] proposed a modification of the algorithm in such a way that it becomes more effective while simultaneously maintaining 2nd order accuracy. In order to examine the introduced modification, let us consider the solution of the original equation (Equation (8)) at time t + ∆t using an unsplit method. This solution can be rewritten in the form [25,33]: Similarly, one can rewrite the solution for the Strang method: h(x, z, t + ∆t) = e 0.5∆tA e ∆tB e 0.5∆tA h 0 (18) and for the Godunov method: Strang observed [27], that while solving the differential equation on a time period (0, T), the calculations are usually not performed for a single large step ∆t, but for many steps that are systematically added up and lead to the total period of time T. As a consequence, applying the Strang method for N T steps (so that N T ·∆t =T) we obtain the following solution [25]: Then using the substitutions e 0.5∆tA e 0.5∆tA = e ∆tA the final solution can be rewritten as follows: h T = e 0.5∆tA e ∆tB e ∆tA e ∆tB e ∆tA . . . e ∆tB e 0.5∆tA h 0 As it can be seen, due to this modification, the sequence of integration is performed with a full time step ∆t for A and B operators. Only at the beginning as well as at the end of the period T, the algorithm is executed with a half time step (Figure 2b). These observations allow the construction of the modified Strang method, which, applied to the directional splitting of the 2D Richards equation, leads to the following sequence of calculations: where the total number of stages S corresponding to the number of time steps N T can be calculated by the following relation: Figure 3 shows the computational efficiency for the Strang splitting methods with respect to the number of time steps N T . To this purpose, the relative increase of the computational effort with reference to the Godunov method was estimated using the following formula: where S ref is the number of stages performed using the Godunov method (S = 2N T ), while S is the number of stages for the standard Strang method (S = 3N T ) or for the modified Strang method (S according to Equation (27)). As illustrated in Figure 3, the standard Strang method requires 50% more calculations than the Godunov method, while the modified Strang method yields a slight increase in the calculation effort (25%) only at the initial calculation period, and it quickly decreases to 5% after 10 time steps. A further increase in efficiency is approximately asymptotic and after 50 steps the growth of computational effort reaches a value below 1.0%.

Alternate Splitting Method
As it can be seen, the modified Strang method differs from the Godunov splitting method only due to the fact that at the beginning and at the end of the computational period T, the algorithm is executed with a half time step (Figure 2b). Moreover, an appropriate result can be expected with the application of the modified Strang method if the time step is constant. However, in practical applications, the Richards equation is solved using algorithms with a variable time step, usually taking very small values (even about 10 −5 s) at the beginning of the simulation. Consequently, in such conditions, the solution obtained by means of the modified Strang method will converge to a solution with the accuracy obtained using the Godunov method. For this reason, an alternate modified algorithm for dimensional splitting is proposed. Based on the modification of splitting for the Strang method, it seems that an increase in accuracy can be achieved by successively reversing the order of integration at each time stage (Figure 4), according to the following algorithm: ∂θ (4) ∂t = Ah (4) over ∆t for j = 1, . . . , N (32) It follows from the above formulations that at a given time level (t = n·∆t), the calculations start in direction z (operator A), while at the next time level (t = (n + 1)·∆t) the reverse sequence is performed starting in direction x (operator B). All integrations are carried out with the full time step ∆t, and consequently, the computational efficiency corresponds to that obtained with the Godunov method. It is worth noting that such a sequence of calculations is valid for the appropriate number of time steps N T corresponding to an even number of stages S performed over the total simulation period T.

Solution of 1D Equations
The specificity of the proposed method used for dimensional decomposition allows the obtained 1D equation to be solved by means of any algorithm based on the FD, FE, or FV methods. In this study, for spatial approximation, the FD method is applied using a combination of backward and forward differences for the diffusion term while using a forward difference for the advection term when an equation in z direction is approximated. Due to the spatial discretization, the partial differential equation is reduced to a system of time-dependent ordinary differential equations (ODE). For this reason, the obtained ODE is solved using the two-level implicit scheme [34,35]. Finally, as a result of the spatial and temporal approximation, the following FD scheme corresponding to the 1D Richards equation in z direction (Equation (11)) is obtained: where i is the index of the internal node; i = 2,3,...,M−1,M is the total number of nodes; n is the index of the time level; ∆z is the space step; ∆t is the time step; θ n i is the nodal value of the water content; h n i is the nodal value of the pressure head; K n i is the nodal value of hydraulic conductivity; K n i+1/2 is the nodal value of hydraulic conductivity averaged using the arithmetic mean [36][37][38], K n i+1/2 = 0.5(K i + K i+1 ); b is a parameter; b = 1 if the advection term is incorporated into the equation for z direction, otherwise b = 0; and η is the weighting parameter ranging from 0 to 1.
In Equation (36), the weighting parameter η determines the stability and accuracy of the numerical scheme [35]. Its particular value gives a well-known algorithm used for the solution of ODE. For η = 1, Equation (36) corresponds to the implicit Euler (backward) scheme, which is usually used in the case of the Richards equation. Such an approximation of the time derivative is of the 1st order of accuracy, and consequently, the numerical diffusion can cause excessive smoothing of the propagating front in the final solution [36,38]. For η = 1/2, Equation (36) becomes the implicit trapezoidal scheme (also known as the Crank-Nicholson scheme), which ensures an approximation with the 2nd order of accuracy when a linear equation is considered. However, the trapezoidal scheme can generate numerical oscillations, which lead to instabilities.
The system of algebraic equations (Equation (36) for i = 2, 3, . . . , M−1) is non-linear. For its solution, Celia et al. [39] proposed the modified Picard algorithm, which minimizes the mass balance error. According to this method, a nodal value of the water content at the next time n + 1 and iteration k + 1 is approximated using the following formula: where n is the time level index; k is the iteration index; C i n+1,k is the nodal value of specific capacity.
Using the above approximation, the modified Picard scheme applied to Equation (36) gives the following algebraic non-linear equation: where ∆h n+1,k+1 is the correction of the pressure head,∆h n+1,k+1 = h n+1,k+1 − h n+1,k ; and f k i is the right side term at the kth iteration described as follows: Equation (38), written for each node i = 2, 3, . . . , M−1 and completed with boundary conditions, can be expressed using the matrix notation: where D is the tri-diagonal matrix of the size M × M; ∆h k+1 is the correction vector, . . , f M ] T is the vector with components described using Equation (39). While solving the set of non-linear equations (Equation (40)), the iteration process is continued until the following criterion of convergence is satisfied for all nodes: where ε is the convergence tolerance, whose value is specified. The presented procedure is repeated for each column j (where j = 1, 2, . . . , N), so that to obtain a solution for whole domain in z direction. The same algorithm is also applied for the 1D equations in x direction (Equation (12)) for each row i (where i = 1, 2, . . . , M). For this reason, in Equation (38), index i and space interval ∆z are replaced with j and ∆x, respectively, while the matrix D is of dimension N × N. Moreover, since an advective term vanishes for the equation in x direction, parameter b is equal to 0.
In order to improve the computational efficiency, the time step is frequently adjusted according to the current number of iterations performed with the iterative procedure. For this reason, the following algorithm can be used [6,11]: Thus, if the number of iterations N iter for ensuring convergence at each time level is higher than the predetermined maximal value N max , then the time step is reduced by the factor f r < 1. On the other hand, when the number of iterations is less than N min , then the time step is multiplied by the factor f m > 1. Usually, satisfactory results for different flow conditions are obtained with N min = 3, N max = 7, f r = 0.6-0.85, and f m = 1.1-1.25 [11,23].

Numerical Tests
For evaluating the splitting algorithms, two numerical tests have been performed where infiltration into dry isotropic soil is considered. In order to check the accuracy of different splitting strategies, the numerical solutions are compared to an exact solution or the reference solution computed for a refined computational grid with significantly smaller spatial and time steps. A comparison is performed in terms of the root mean square (RMS) error computed at the end of the simulation. This error is defined as follows: where RMS is the root mean square error calculated at the final simulation time T; h i is the calculated nodal value of the pressure head at time T; h re f i is the reference or exact nodal value of the pressure head at time T; and M is the number of nodes.
The obtained results are also checked with regard to the mass conservation principle. In this study, the relative mass balance error is used in the following form: where ∆M is the relative mass balance error; and V T and V 0 are the total final (at time t = T) and total initial (at time t = 0) amount of water in the domain given by the formulas: where F b is the total net flux through the boundaries of the domain: Moreover, in the study, the computational efficiency is evaluated based on the total number of iterations N iter performed during the simulation time T. It is worth adding that the efficiency depends on many factors. For example, the efficiency can be increased by increasing the time step or/and spatial step. Consequently, changing the time step and the resolution of the grid, at the same time also affects the accuracy of the final solution. For this reason, it seems that the computational efficiency should be evaluated using both the total number of iterations (or the computational time) and the accuracy of the solution [40,41]. In this study, for the estimation of the computational efficiency, the following formula is proposed: where Ef is the computational efficiency; RMS is the solution accuracy calculated using the RMS error (Equation (43)); and N iter is the total number of iterations N T .

Test 1-Simulation of Infiltration Using the Gardner Model
In the first test, the hydraulic properties in the Richards equation are described using the Gardner functions given by Equations (5)-(7) with the following parameters: residual water content θ r = 0.15, saturated water content θ s = 0.45, saturated hydraulic conductivity K s = 10 −5 m/s, and scaling parameter h g = 2 m. The calculations are performed for the rectangular domain with the dimensions L x = 1 m and L z = 2.5 m. At the top of the domain for z = 0 as a boundary condition, the following distribution of the pressure head is imposed: Thus, the pressure head on the top is equal to zero in the center (L x /2, L z /2) and decreases to the value h b on the left (x = 0, z = 0) and right (x = L x , z = 0) edges (Figure 5a). On all boundaries of the domain, a constant pressure head is assumed: The initial condition at the time t = 0 corresponds to a uniform distribution of the pressure head for the whole domain and is specified as h(x,z,t = 0) = h 0 = −10 m. For such initial−boundary conditions, the analytical solution of the 2D Richards equation has the following form [31]: Numerical tests were performed for a grid with a constant spatial length ∆x = ∆y = 0.025 m. The analytical solution was compared with the results obtained using dimensional splitting algorithms in the form of the Godunov method (G), Equations (11)-(12); the Strang method (S), Equations (14)- (16); and the modified Strang method (mS), Equations (22)- (26). The calculations were carried out for two variants, depending on the sequence in which the 1D equations are solved at each time level. In the first variant, at the initial stage, the calculation always starts with z direction (ZX), while for the second variant, a reverse sequence is performed with the initial stage in x direction (XZ). Simulations were driven for two values of the weighting parameter: η = 0.5 (the implicit trapezoidal scheme) and η = 1 (the implicit Euler scheme) in the FD scheme (Equation (36)) and using a constant time step ∆t equal to 5, 10, 20, 40, and 50 s. While solving the system of non-linear algebraic equations, the tolerance of the iterative process ε was set to 0.001 m (Equation (41)). Such assumed value of tolerance ensures high accuracy while maintaining the number of iterations at an appropriate level usually not exceeding twenty iterations at a given time step.
A graphical illustration of the results calculated at the end of the simulation T = 1000 s is presented in Figures 5 and 6. Figure 5a displays the pressure head distribution in the form of a contour map obtained for the exact solution, whereas in Figure 5b-d, the corresponding isolines of water pressure, calculated using various splitting methods with the time step ∆t = 50 s, are compared with the exact solution. Moreover, in Figure 6, the water saturation profiles along the vertical z direction at the cross-section x = 0.5 are presented. It is apparent that the Godunov splitting method (Figures 5b and 6b) reproduces solutions where the position of the wetting-drying front is underestimated (G-ZX) or overestimated (G-XZ), depending on the sequence of calculations. The difference between the exact and numerical solutions is the largest for the calculations performed with the sequence XZ. The best results are obtained with the sequence ZX for both the standard (Figures 5c and 6b) and the modified Strang method (Figure 5d). The calculations performed with the sequence XZ result in a greater deviation from the exact solution, but in this case, the difference is smaller than for the Godunov splitting method. It seems that the discrepancy between the results obtained with the ZX and XZ sequences is caused by the fact that the 1D Richards equations after dimensional decomposition differ by the advection term resulting from gravity. Consequently, while splitting, the final result depends on the sequence in which the 1D equations are numerically integrated. The above results are also confirmed by the accuracy analysis in terms of the RMS errors presented in Figure 7. RMS errors are calculated using the solutions along the axis z at the cross-section x = 0.5 m obtained for various constant time steps ∆t equal to 5, 10, 20, 40, and 50 s. Moreover, the RMS error analysis is carried out for the solutions obtained using the implicit trapezoidal scheme with the weighting parameter η = 0.5 (Figure 7a) and the implicit Euler scheme with η = 1 (Figure 7b). Moreover, the standard and modified Strang splitting methods give very similar values of RMS errors if the trapezoidal scheme is applied (Figure 7a). The obtained results indicate that the Strang method and its modification are approximations of the 2nd order of accuracy and they minimize the splitting error when dimensional decomposition is applied for the Richards equation. However, when the implicit Euler scheme with η = 1 is used (Figure 7b), the solutions are burdened with a higher RMS error, and the consistency between the standard and modified Strang splitting methods is also not observed, as was the case for the solution with η = 0.5. These results confirm that in addition to the appropriate splitting technique, the accuracy of the resulting solution is also affected by the choice of the time integration scheme for solving ODE. The applied trapezoidal scheme is an approximation of the 2nd order of accuracy (with respect to the time) and, combined with the Strang methods, leads to the best results. While the implicit Euler algorithm, in conjunction with the Godunov splitting method, ensures only the 1st order of accuracy leading to a worse solution. Figure 8 presents the time evolution of the mass balance error. It is apparent that for the solution obtained using the standard Strang method with the weighting parameter η = 0.5 (Figure 8a), the mass principle is satisfied regardless of the simulation time t (the errors ∆M are less than 0.5%). Unfortunately, other splitting algorithms do not satisfy the mass conservation principle, because the mass balance errors systematically increase during the simulation time, reaching even several percent. In these cases, the mass balance error takes a negative or positive value depending on the sequence of computations. In Figure 8b, the results for the solution with the implicit Euler scheme (η = 1) are presented. As it can be seen, the mass balance values significantly differ from those obtained for the trapezoidal scheme. Thus, this test shows that besides the method of dimensional splitting, the time integration scheme also determines the mass conservation property. Since the implicit Euler scheme is the 1st order of accuracy, it generates a numerical diffusion which causes an additional deviation of the mass balance in a non-conservative algorithm. It is worth adding that this effect is also observed in other solutions of non-linear equations such as the kinematic wave equations or the shallow water equations, which are written in a non-conservative form [26,42]. The next calculations were carried out to check the efficiency of the splitting algorithms. In Figure 9a, the dependence of the number of iterations N iter on the time step is presented. This experiment confirms the earlier conclusions that the standard Strang method requires a greater number of time steps (and therefore number of iterations) than the Godunov method as well as the modified Strang method. It is also worth noting that for the modified Strang method, the number iterations necessary to achieve a convergence also depends on the sequence of splitting. However, when the accuracy of the solution is also taken into account for the estimation of computational efficiency (according to Equation (48)), the Godunov method is less effective than the Strang methods (Figure 9b). This is especially apparent for greater values of the time step ∆t ≥ 20 s. The best results with regard to such computational efficiency were achieved by means of the modified and standard Strang methods with the sequence ZX of calculations.

Test 2-Simulation of Infiltration Using the Mualem-van Genuchten Model
In this case, a solution of the 2D Richards equation with the Mualem-van Genuchten functions is considered. The calculations are performed for typical properties of sand and loam presented by Carsel and Parrish [43] ( Table 1). The computational domain is rectangular with the dimensions L x and L z . For a part of the top boundary (0.46 <x b < 0.54 m), a constant pressure head is imposed h(x b ,z = 0) = 0 m (Figure 10), while on other boundaries, the zero flux condition in the normal direction is applied. The initial condition is specified as h(x,z,t = 0) = h 0 = −10 m. Since this test is more challenging with regard to the reproduction of a solution containing a sharp wetting front, to ensure the appropriate accuracy as well as the efficiency of the solution, the time step is adjusted according to the procedure described by Equation (42). To this order, the predetermined values of the minimal and maximal number of iterations are generally set at N min = 3 and N max = 7, respectively. The minimal and maximal allowable values of the time step (∆t min and ∆t max ) are specified depending on the considered soil. All comparisons were performed over a fixed simulation time T. The results obtained using different splitting strategies were compared to the reference solution computed using a refined computational grid with a smaller spatial step ∆x = 0.01 m and the maximal time step ∆t max = 1 s. Moreover, the convergence error ε in the iteration process was decreased from 0.001 m to a value of 0.0001 m for the reference solution. Simulations were driven using the trapezoidal scheme with η = 0.5. The values of the hydraulic parameters of the soil as well as other numerical parameters used in Test 2 are summarized in Table 1. For each soil, the computations were performed using the Godunov method (G), Equations (11)-(12); the Strang method (S), Equations (14)- (16); and the alternate modified method (mod), Equations (29)- (35). As in the previous test, the simulations were carried out for two variants of a sequence of splitting (ZX or XZ) in which 1D equations are solved.
The results of the simulation obtained for sand and loam are presented in Figures 10-12, whereas Table 2 summarizes the statistics in regards to the accuracy, mass balance, and computational efficiency. Figure 10 shows the position of wetting fronts after the simulation time T obtained using the splitting algorithms with the ZX sequence only, because the solutions for the XZ sequence are very similar with regard to graphical presentation. In order to enhance the difference between various splitting methods, the predetermined values of the minimal and maximal number of iterations is set on N min = 4 and N max = 8, respectively. Due to this treatment, the time step takes a higher value during simulation than for the values previously assumed (N min = 3 and N max = 7), and consequently, in the numerical solution, greater splitting errors are generated. As it can be seen in the case of sand (Figure 10a) as well as loam (Figure 10b), all splitting algorithms reproduce a wetting front which is slightly overestimated. However, the best result with regard to the RMS error is obtained for sand using the alternate modified method (mod-ZX and mod-XZ), while for loam, all algorithms give very similar values of the RMS errors ( Table 2).  In Figure 11, the mass balance error is presented. As with the previous test, the Godunov method (G-ZX, G-XZ) reveals a low mass conservation for both the ZX and XZ sequence. In this case, the error ∆M is stabilized to the value of the order at 5% and 8% for sand and loam, respectively. While the Strang method (S-ZX, S-XZ) and alternate modified method (mod-ZX, mod-XZ) minimize the mass balance error. This is especially apparent in the case of the latter method, which substantially reduces the mass balance error ∆M to a value below 1% and 2% for sand and loam, respectively. Figure 12 displays the number of iterations for different simulation times. The Godunov method as well as the proposed method with alternate modifications are the most effective and require about 50% less iterations than the Strang splitting method. However, considering the computational efficiency estimated according to Equation (48), the best results are obtained using the alternate modified method ( Table 2). Note that with such defined efficiency, the Strang method is the least effective.

Conclusions
In the paper, the dimensional splitting algorithms for solving the 2D Richards equation have been presented. Due to the splitting procedure, the solution of the 2D equation is reduced to the solution of a set of 1D equations for each of the directions. Consequently, the obtained system of algebraic equations has a tri-diagonal coefficient matrix, which can be solved more efficiently than as is the case with unsplit 2D algorithms. All splitting methods presented in the study allow the use of any algorithm for space as well as time approximation, which undoubtedly is a significant advantage over the methods based on ADI algorithms.
The first of the analyzed methods, the Godunov method, generates a significant splitting error, which consequently reduces the accuracy of the obtained solution. As it was shown, this splitting method also leads to a solution with a significant mass balance error.
In order to increase the accuracy of the solution, the standard and modified Strang methods were also applied for dimensional decomposition. Due to these methods, it is possible to significantly minimize the splitting error, therefore ensuring the accuracy of the solution to the 2nd order. However, the standard Strang method is not computationally effective, while the modified Strang splitting method requires the same computational effort as the Godunov method. On the other hand, the modified Strang method does not satisfy the mass conservation principle and it can be used only when the time step is constant. Numerical experiments also indicate that solutions obtained using splitting algorithms can be overestimated or underestimated depending on the sequence in which 1D equations are solved at each time level. This effect is especially apparent for the Godunov method.
In view of the above inconveniences, an alternate splitting method has been proposed which is based on successively reversing the order of integration at each time stage. The conducted experiments showed that the proposed method ensures adequate accuracy of the solution and minimizes the mass balance error. At the same time, it is the most effective of all the analyzed splitting methods, in terms of both the number of iterations and the efficiency determined, taking into account the accuracy of the solution (Equation (48)).
In the analysis, a 2D unsaturated flow for a simple domain geometry using a rectangular grid was considered. However, the proposed splitting method can also be successfully extended to 3D flow problems. Moreover, it seems that this approach can also be used for more complex geometries on an irregular grid if a coordinate transformation method is applied.