Accurate Empirical Calculation System for Predicting the Seepage Discharge and Free Surface Location of Earth Dam over Horizontal Impervious Foundation

In the numerical analysis of the Laplace equation, which is the governing equation of the seepage phenomena of homogeneous, isotropic earth dams, it has been confirmed that numerical analysis with high accuracy is possible by using the interpolation finite difference method (IFDM). In a previous paper, based on this numerical analysis method, the equivalent Kozeny (KZ) flow method was proposed as a new empirical method to calculate the seepage discharges and free surface locations of earth dams. Although this method is generally a highly accurate method compared with the empirical method of A. Casagrande, owing to calculating the seepage problems within a few percentages of discharge relative errors, several additional studies are necessary. By integrating the finding of this study to the previous literature, an empirical seepage calculation system with high accuracy, the equivalent KZ flow method, is created. Owing to the finally proposed empirical method, called “interpolation-equivalent KZ flow method”, the discharge and free surface location can be predicted with high accuracy in a wide range.


Introduction
Conventionally, high-accuracy numerical analyses of seepage problems of earth dams have been implemented using the finite element method. Numerical analyses using the finite difference method (FDM) have been limited to cases where the calculation domains are comparatively simple. However, by applying the interpolation finite difference method (IFDM), two-and three-dimensional elliptic partial differential equations (PDEs) over complex domains with high speed and high accuracy can be freely solved [1][2][3] The IFDM is composed of two kinds of methods [3] (1) Algebraic Polynomial Interpolation Method (APIM), where finite difference schemes are formulated instantaneously over equally or unequally spaced grid points, and (2) Boundary Polynomial Interpolation Method (BPIM), where numerical calculations are executed only over equally spaced grid points, and the boundary interpolation is performed to match the boundary conditions. In the paper, the numerical calculations are executed by the IFDM-BPIM.
This method is adopted as a numerical analysis method for solving confined/unconfined seepage problems. In such problems, because of the adoption of IFDM-BPIM using the quadratic interpolation at the boundary, it has been shown that highly accurate numerical calculations are possible even in the FDM [4,5].

Governing Equation of Numerical Calculation
In our previous papers [4,5], numerical calculation methods based on the IFDM-BPIM were explained in detail. Here, the governing partial differential equation (PDE) and its FDE are verified. In the IFDM, the PDE for the potential calculation is the following parabolic PDE: where h is the total head (position head, y, +pressure head, p) and t denotes time; then, ∂h/∂t is the pseudo-acceleration term [5,15]. The x-direction is horizontal and the y-direction is vertical. Under a Eng 2020, 1 62 fixed boundary condition, as t → ∞ , ∂h/∂t → 0, the equation reduces to the Laplace equation, ∇ 2 h = 0. The following expression describes the stream function, s: The numerical calculation methods for Equations (1) and (2) are the same; consider Equation (1). The Laplace equation, ∇ 2 h = 0, is usually calculated using the successive over relaxation (SOR) scheme. On the other hand, Equation (1) is usually solved by the forward time, centered space (FTCS) scheme. The FDE of the FTCS scheme is expressed as follows [4]: (∆y) 2 (3b) In the FTCS scheme, calculated variables are replaced simultaneously. However, if we are only interested in the convergence state, successive displacement of variables accelerates the convergence. This scheme is the time marching successive displacement (TMSD) scheme [2,4,15]. Changes from the former to the latter are carried out easily and instantaneously. Moreover, it has been confirmed that the TMSD scheme is exactly equivalent to the SOR scheme [16] limited as using the second-order centered space finite difference scheme [2,15]. In the paper, all calculations are carried out using the TMSD scheme.
In Equation (3), ϕ i,j is the time-interval adjustment factor, and l wi,j is the wall-distance factor.
In the definition of ϕ i,j ≡ l wi,j p w in Equation (3), p w = 1 is defined under the condition that the calculation element of concern is the near-wall element with 0 < l wi,j < 1, and the boundary condition is the Dirichlet condition. In all other cases, ϕ i,j = 1 (namely, p w = 0) and Equation (3) reduce to the usual FDE [4]. The time interval is defined as follows: where α b is the acceleration factor of the TMSD scheme, and α bmax is the theoretical maximum acceleration factor. The stability condition 0 < α b < α bmax (= 2) ensures that the calculation has a convergent solution. It is advantageous to use the optimum acceleration factor, α b = α bopt , to achieve the fastest convergence [2]. However, in this study, α b = 1.80 is adopted because this value guarantees stable convergence without exception.

Concept of Equivalent KZ Flow
In the previous paper [5], the concept of equivalent KZ flow was presented in detail. Because the theory employed in this study is intertwined with the concept, this section begins by confirming the basic contents presented in the previous paper. Kozeny (1931) [17] derived an analytical solution of the seepage phenomenon with a discharge angle α = 180 • (see Figure 1) by conformal mapping and showed that the free surface becomes a parabola. A. Casagrande (1937) [18] introduced this theory to his method-the basic parabola method. Even with the proposed equivalent KZ flow method, the Kozeny flow is considered as the basis of the theory; however, there are several variations in defining this flow, which is presumed to be beyond the range intended by Kozeny (1931) [17]. Furthermore, considering the simplicity of frequently appearing the expression, the Kozeny flow, it was referred to as the KZ flow in the previous paper [5].   [5]. In this figure, is defined as follows: The free-surface parabola (FSP) of this flow is expressed in the coordinate system (i.e., KZ coordinate system) specified in Figure 1, as follows: On the other hand, the boundary-potential parabola (BPP) is expressed as follows: When = 0, = + 2 ⁄ ; therefore, | | = 2 = 2.5 ⁄ in Figure 1. Both are parabolas with focus at the origin, O. The potential and stream functions in the domain are expressed by the following analytical solution [5].
In the unconfined seepage problem, it is extremely infrequent to obtain a simple theoretical solution in a finite domain using the conformal mapping method. Therefore, in order to confirm the accuracy of numerical analysis, the KZ flow is an effective numerical analysis target.
Interestingly, in the KZ flow, let us refer to the ratio, = ⁄ , as the basic Ratio. Then, it is clear that two KZ flows with the same basic aspect ratio are similar. The validity of this proposition is directly confirmed by deriving the dimensionless potential, h/H from Equation (8) and the dimensionless stream function s/H from Equation (9). They become as follows: ℎ * = + 1 − * + * + * , ℎ * = ℎ , * = , ⁄ * = ⁄ ⁄ (10)  Figure 1 shows an example of the flow net of a KZ flow [5]. In this figure, S is defined as follows: The free-surface parabola (FSP) of this flow is expressed in the coordinate system (i.e., KZ coordinate system) specified in Figure 1, as follows: On the other hand, the boundary-potential parabola (BPP) is expressed as follows: When y = 0, x = d + S/2; therefore, GG d = S/2 = 2.5 in Figure 1. Both are parabolas with focus at the origin, O. The potential and stream functions in the domain are expressed by the following analytical solution [5].
In the unconfined seepage problem, it is extremely infrequent to obtain a simple theoretical solution in a finite domain using the conformal mapping method. Therefore, in order to confirm the accuracy of numerical analysis, the KZ flow is an effective numerical analysis target.
Interestingly, in the KZ flow, let us refer to the ratio, R b = d/H, as the basic Ratio. Then, it is clear that two KZ flows with the same basic aspect ratio are similar. The validity of this proposition is directly confirmed by deriving the dimensionless potential, h/H from Equation (8) and the dimensionless stream function s/H from Equation (9). They become as follows: Eng 2020, 1 So, if we reduce Figure 1 to 1/10, this is the KZ flow at H = 1, and the discharge also becomes 1/10. This concept of similarity is important in this paper. The potential calculation is obtained irrespective of the permeability coefficient under the assumption of homogeneous material. In the stream function, Equation (9) must be multiplied by the permeability coefficient in order to correspond to the real discharge. However, in this paper, the result of setting the permeability coefficient k c = 1 is formally assumed; then, no problem occurs and the discharge of the KZ flow is expressed as q = S = y 0 , the ordinate value of point D in Figure 1.
Equation (6) is expressed as x = ay 2 + c. That is, if two independent conditions are given, this expression is determined. The KZ flow is determined by the point A(d, H) and the focal point O(0, 0). This is described as KZ0(A, O)-the basic expression of the KZ flow. There are three other definitions of the KZ flow under the condition that point A is given: (1) the area of the KZ flow, A K , is given, and this is described as KZ1(A, A K ); (2) the other arbitrary point B(xB, yB) is given, and this is described as KZ2(A, B); (3) the discharge is given, and this is described as KZ3(A, q) where q = S; finally (iv) both the discharge and discharge point are given, and this is described as KZ4(q, B). For KZ1(A, A K ), the area of KZ flow, A K , must be determined. In Figure 1, the KZ flow area is divided into two: A(I) and A(II). Using Equation (6), area A(I) is obtained as follows: Using Equation (7), we have the following: Therefore, the total area is as follows: This yields the following: Then, the focal point, O 1 , is determined, that is, the KZ flow is determined. For KZ2(A, B) in the KZ coordinate system, the parabola is expressed as x = ay 2 + c; thereafter, the condition that points A(d, H) and B(xB, yB) are on the parabola determines the values of a and c. From parabola theory, the focal point and discharge are O 2 (c + (1/4a), 0) and S = 1/2a, respectively.
For KZ3(A, q), q(= S) is given. In the parabola x = ay 2 + c, d = aH 2 + c because point A is on the parabola. The discharge q = 1/2a is given. Then, the parabola and focal point O 3 (c + (1/4a), 0) are determined.
Finally, for KZ4 (q, B), in the parabola x = ay 2 + c, xB = ayB 2 + c because point B is on the parabola. The other description is the same as that of KZ3(A, q), and the focal point O 4 is determined.
All the definitions, (1), (2), (3), and (4), are simply different expressions of KZ0(A, O); however, each of them has a significant function in the equivalent KZ flow method.
In Figure 1, when the drain is inclined, such as OG , the analytical solution of KZ flow becomes invalid. The free surface location moves downward as it approaches the discharge point, B , which is not on the FSP. As an approximate rule of thumb, we found that the concept of area-equivalent KZ flow is effective [5]. The deficit area, A de (OFB ), cut by the inclined line OG is determined depending on the discharge angle, α = ∠GOG . Let the equivalent area be defined as A e = A k − A de . Then, the equivalent KZ flow having the entire area, A e , is determined with point A as a fixed point. By substituting A e for A k in Equation (15), the equivalent distance, d e , is determined; specifically, the equivalent origin is determined, and the area-equivalent KZ flow, KZ1(A, A e ), is defined. In the previous study, KZ1(A, A e ) proved to be a substantially more advantageous empirical method than A. Casagrande's basic parabola method. However, with the goal of approximately 3% permissible error in this study, this rule of thumb must be reconsidered, especially in the range of α < 90 • . Figure 2 shows the numerical calculation result of the case shown in Figure 1 (Fukuchi, 2018). AlongP * 1 P * 11 , the theoretical solution of BPP is discretely given. The position coordinates of the points of P * k , k = 1, 2, · · · , 13 and the length of the horizontal drain are input. In the definition of P k , k = 1, 2, · · · , P 35 , P k = P * k , k = 1, 2, · · · , 11 and the points P k , k = 12, 13, · · · , 35 are determined in the calculation process. In particular, the notation P * k is used for defining the dam configuration; on the other hand, the notation P k is directly used for numerical calculation. In Figure 2, the equipotential lines are 0.1, 0.2, · · · , 0.9H (division number, N p = 10), and the streamlines are 0.2, 0.4, · · · , 0.8q (division number, N s = 5), indicated by solid lines. The theoretical equipotential lines are indicated by broken lines with N p = 20 and N s = 10 for streamlines. In the IFDM calculation, the active calculation elements consist of three kinds of element: normal element (colorless), near-wall element (yellow), and dummy element (gray). Numerical calculation is performed at grid points, which are centroids of each calculation element. The numerical solution agrees well with the theoretical solution. From the result of the potential calculation, the flow velocity in the horizontal direction (u) is defined at each grid point, such as u i, j = −k c h i+1, j − h i−1, j /(2∆x). Assume that k c = 1. In the range of |G d O|, 15 individual discharge values are calculated. Each discharge is calculated by the parabolic numerical integration method [1,5], then the value at either end is obtained by extrapolating the corresponding quadratic curve; this ensures sufficient accuracy. The discharge becomes q n = 4.995(−0.09%) ± 0.058(= 3σ); σ is the standard deviation. Although slight dispersions of individual discharges are inevitable, the relative error of the calculated mean discharge is −0.09% with respect to the theoretical value q t = 5.00; this is highly accurate, being less than the 0.1% (absolute) relative error. Such features of the discharge calculation are common in the IFDM calculation.
In the numerical calculation process, it is necessary to estimate the boundary points on the right side of point O. It is assumed that the curveP 23 P 21 P 13 is expressed by x = ay 2 + c. This equation is derived from P 23 (x 23 , y 23 ) and P 21 (x 21 , y 21 ). Each point of P 20 , P 19 , · · · , P 14 is on the curve. The y-coordinate of each boundary point is defined as y 21 × (i/N), i = N − 1, N − 2, · · · , 0. Here, N = 8; the division number, N, is defined such that the calculation points are set to have appropriate intervals [5]. However, this is an example of α = 180 • ; if α ≤ 90 • , only one point is estimated. In this case, two equations for estimating point B are conceived: (1) linear extrapolation, y = ax + b and (2) quadratic extrapolation, y = ax 2 + bx + c. Generally, it makes no large difference as to which one is used; however, in the study, extrapolation is performed by (3) the KZ flow extrapolation.   Dirichlet condition is also given; however, when the discharge angle is 90°< α <180°, it is a position-dependent variable and is designated (DV). Along 4 ̅̅̅̅̅ , although the condition is the Neumann condition, it is denoted as NS to indicate that it is a free surface part. Regarding the stream function, the Neumann condition (N) is specified along 1 2 ̅̅̅̅̅̅ , and the Dirichlet condition (0.00) is designated along 2 3 ̅̅̅̅̅̅ . Along 3 4 ̅̅̅̅̅̅ , a Neumann condition orthogonal to the equipotential line is set; it is generally expressed as NV because the equipotential line at the boundary point of concern must be individually calculated according to each position. The discharge should be given as Dirichlet condition along 4 ̅̅̅̅̅ , but there is no problem specifying 1.00 as the Dirichlet condition. By multiplying the numerically calculated discharge, (as the result of the potential calculation), the correct conclusion is always obtained.  The global coordinate values of input points and boundary conditions of the potential and stream function calculations are shown in the input data. In the potential calculation, along P 1 P 2 ≡ AP * 2 , the Dirichlet condition (8.00) is given; along P 2 P 3 ≡ P * 2 O , the Neumann condition (N), where the equipotential line is orthogonal to the boundary line, is given. Along P 3 P 4 , the Dirichlet condition is also given; however, when the discharge angle is 90 • < α < 180 • , it is a position-dependent variable and is designated (DV). Along P 4 A, although the condition is the Neumann condition, it is denoted as NS to indicate that it is a free surface part. Regarding the stream function, the Neumann condition (N) is specified along P 1 P 2 , and the Dirichlet condition (0.00) is designated along P 2 P 3 . Along P 3 P 4 , a Neumann condition orthogonal to the equipotential line is set; it is generally expressed as NV because the equipotential line at the boundary point of concern must be individually calculated according to each position. The discharge should be given as Dirichlet condition along P 4 A, but there is no problem specifying 1.00 as the Dirichlet condition. By multiplying the numerically calculated discharge, q n (as the result of the potential calculation), the correct conclusion is always obtained.     Figure 3b, the point A c is given so that the areas of A(II) ( Figure 1) and triangle AA c G d are equal; this practically coincides with P * 2 . Then, the entrance angle, θ = θ c , is defined as the criterion entrance angle. The starting point C of FSP(c) matches point A. Moreover, the free surface of numerical calculation, FS(n), and FSP(c) coincide and the discharges would be the same. It is confirmed that the results are q n = 2.51 ± 0.06, q c (= y c ) = 2.52 (0.5%). On the other hand, in A. Casagrande's basic parabola method, it becomes q Ca = 2.44(−3.1%) [5]. It is understood that his method determining point C a , C a A = 0.3 P * 2 G d is incorrect. Figure 4 is a general case where point P * 2 does not match point A c ; fundamental results are shown in the caption. In A. Casagrande's method, the position of point C a is determined based on P * 2 G d ; however, to be correct, it should be determined on the basis of L c = P * 2 A c , as shown in Figure 4b. As a result of the detailed examination, in the definition CA = m v L c , m v is determined as follows [5]: Therefore, point C is determined, and the basic KZ flow, KZb(C, O), is defined. Because q b = 1.85 (−0.1%) for the numerically calculated discharge, q n = 1.86 ± 0.05, it can be seen that the calculation is highly accurate. (In this paper, the exact solution is defined only for the KZ flow and "rectangular dam flow" (described later). Therefore, the relative errors in the empirical methods below are evaluated by the value of the numerical calculation result.) If the entrance angle is θ = 90 • , R c = L c /H < 0, then Equation (16a) is used. Through the calculations up to this point, the discharge q and discharge point B were calculated. Let us define this as the empirical QB estimation. The estimation of A. Casagrande is of course an empirical QB estimation.  A remaining problem is to express accurately the free surface from FSP(b). In Figure 4b, let the deviation width of FSP(b) at point A be ; moreover, = | |. In this case, matches . Then, the deviation must be corrected between 0 ≤ < ; that is, the FSP(b) must be modified to conform to the FS(n), called the modified FS(b). The modification equation is defined in the previous paper [8] as follows: where ℎ * * is the modification value, * = / is the dimensionless distance, and * ≡ ⁄ is the dimensionless half-decrease distance ( is the half-decrease distance). In this expression, at * = 0, it becomes ℎ * = ; at * = 1 , it becomes ℎ * = 0 . The dimensionless half-decrease distance, * , generally depends on the ratio, ≡ ⁄ ; * is defined as follows: A remaining problem is to express accurately the free surface from FSP(b). In Figure 4b, let the deviation width of FSP(b) at point A be δH; moreover, Then, the deviation must be corrected between 0 ≤ l < L; that is, the FSP(b) must be modified to conform to the FS(n), called the modified FS(b). The modification equation is defined in the previous paper [8] as follows: where δh * l * is the modification value, l * = l/L is the dimensionless distance, and L * h ≡ L h /L is the dimensionless half-decrease distance (L h is the half-decrease distance). In this expression, at l * = 0, it becomes δh * 0 = δH; at l * = 1, it becomes δh * 1 = 0. The dimensionless half-decrease distance, L * h , generally depends on the ratio, R h ≡ L/H; L * h is defined as follows: The detail is described in the previous paper [5]. Equation (18) is formulated based on the numerical calculation results of the rectangular dam with R b ≥ 0.2. It should be noted that in the case of other cross sections, some deviation may occur. In Figure 4b, it is confirmed that the green line AF b : modi f ied FS(b) practically coincides with the calculated free surface FS(n). The discharge was calculated and at the same time the free surface was defined by modifying FSP(b). Let us define this as the empirical QF estimation. The empirical QB estimation becomes the corresponding empirical QF estimation by using the modification Equation (17). The parabola FSP(b) is called the "basic parabola" in the equivalent KZ flow method as well.
When α = 50 • in Figure 4, the calculation result is shown in Figure 5. The numerical calculation yields q n = 2.00 ± 0.01 and yB n = 2.84. In A. Casagrande's method, the calculation gives q Ca = 1.73 (−13.3%) and yB Ca = 2.49(−4.3%). As α decreases, the results given by A. Casagrande's method rapidly deteriorate. The position of point C in Figure 5b is the same as that in Figure 4b. Using the area-equivalent KZ flow, KZ1(C, A e ), the results are q 1 = 1.96(−2.0%) and yB 1 = 3.29(5.7%). In terms of discharge, the result is considerably better than A. Casagrande's method; however, as a whole, their (absolute) relative errors are too large. In the range 90 • ≤ α ≤ 180 • , KZ1(C, A e ) yields good results; however, as α decreases from 90 • , it gradually yields poor conclusions [8]. Here, it should be noted that q 2 = 2.01(0.3%) in the discharge point-equivalent KZ flow, KZ2(C, B n ). The position of point C varies by nature depending on the change in α; however, it can be observed that the change is practically negligible in this example. Figure 5c shows the numeric-equivalent KZ flow, KZn(q n , B n ). The calculated free surface FS(n) is almost precisely expressed by the modified FS(n). In the figure, the equivalent vertical face, C E C E , is shown, which is set in order that the area of A(II) of KZn(q n , B n ) is equal to the area of the rectangle C E C E C n C n . The discharge (q E ) and discharge point (yB E ) of the equivalent cross section, C E C E OG, which are shown in detail later, practically coincide with those of KZn(q n , B n ). The numerical calculation yields the results q E = 1.99(−0.7%) and yB E = 2.79(−0.6%); their relative errors are both within 1%. The equivalent aspect ratio of this equivalent trapezoidal dam is R E = C E O /H = 2.130. The concept of equivalent aspect ratio has an important role in the following discussion.

Calculation Accuracy of Conventional Empirical Methods
As for the calculation accuracy of A. Casagrande's method, because of confirming the calculation accuracy for each calculation example in the previous paper [5], it was generally found that nonnegligible deviations from the numerical calculation results occur. In this paper, the methods of Schaffernak-Van Iterson and L. Casagrande regarding the seepage phenomenon of the trapezoidal dam with a vertical entrance face are discussed. These methods are included in the calculation system of A. Casagrande's method [18]. The governing formulae and numerical calculation techniques of both these methods are shown in detail in Appendix A and Appendix B. Regarding discharge angle,

Calculation Accuracy of Conventional Empirical Methods
As for the calculation accuracy of A. Casagrande's method, because of confirming the calculation accuracy for each calculation example in the previous paper [5], it was generally found that non-negligible deviations from the numerical calculation results occur. In this paper, the methods of Schaffernak-Van Iterson and L. Casagrande regarding the seepage phenomenon of the trapezoidal dam with a vertical entrance face are discussed. These methods are included in the calculation system of A. Casagrande's method [18]. The governing formulae and numerical calculation techniques of both these methods are shown in detail in Appendices A and B. Regarding discharge angle, the scope of the Schaffernak-Van Iterson method was set to α < 30 • , and the scope of L. Casagrande's method was set to α ≤ 90 • . Here, apart from such an empirical rule, the theoretical characteristics when both equations are applied to the seepage of a trapezoidal cross section α ≤ 90 • are considered.
In general, the discharge, q, and y-coordinate, yB, of discharge point B can be expressed by the following functional definitions: Eng 2020, 1

71
Here, θ is the entrance angle, R b (= d/H) is the basic aspect ratio, and α is the discharge angle ( Figure A1). These equations show the (geometrical) similarity law of seepage problems. From now on, let us set H = 1.0. Under the conditions θ = 90 • and α ≤ 90 • , let us investigate the characteristics of the calculation error when using the two types of empirical methods: the Schaffernak-Van Iterson method and L. Casagrande's method. As already mentioned, the similarity law holds in KZ flow. This must be true for all the empirical methods as well.
Occasionally, the overall physical phenomena image becomes distinct by considering the extreme states. Under the condition R b = 2.0, the numerical analyses are conducted, shown in Figures 6-9, and the two extreme states of Figures 6 and 9 are obtained. In Figures 7-9, the calculation results using A. Casagrande's method are also presented. Figure 6 shows the one extreme state, where the discharge angle (α) is the minimum: α = α min = arctan(1/R b ) 26.6 • ; accordingly, the entrance point (A) and discharge point (B) agree. This problem is a special confined seepage problem, and the numerically calculated discharge becomes q n = 0.500. In the figure, equipotential lines and streamlines are shown. The calculation area is X[0.0, 2.0], Y[0.0, 1.0], and the finite difference width is ∆x = ∆y = 0.05. In the calculation examples of H = 1.0 presented in this paper, the calculation element disposition and finite difference width are always used in this specification. Except for Figure 6, illustrations of calculation element dispositions are not presented. Despite the theoretical importance of this seepage problem, no description pertaining to its analytical solution in the leading references (e.g., Polubarinova-Kochina, 1962 [19], Harr, 1991 [20]) is found. However, this problem is regarded as part of the seepage phenomenon, particularly in the case where the head of P 1 P 2 is 1.00 and that of P 3 P 4 is 0.00 in the confined seepage problem of domain P 1 P 2 P 3 P 4 ; accordingly, no contradiction occurs. Therefore, its theoretical solution is q t = 1/R b = 0.500. Equation (A2), which is the discharge expression in the Schaffernak-Van Iterson method, would also be effective in this case. Then, this equation becomes q SV = a sin α tan α = a/a tan α = 1/R b = 0.500 (0.0%), which agrees with the theoretical solution. Even in L. Casagrande's method, Equation (A12) would be an effective definition of the discharge; accordingly, q LC = a sin 2 α = sin α = 0.447(−10.6%). Each method represents a one-dimensional approximation of the seepage phenomenon. The physically reasonable equation in the extreme state of Figure 6 is evidently the former expression. When R b → ∞ , q LC → q SV (sin α → tan α) , and, generally, q LC < q SV holds true. Figures 7 and 8 show the results when α = 30 • and α = 60 • , respectively. In the figures, the values of R u ≡ P * 1 P * 4 , referred to as up-side aspect ratio, are shown. In the process of R u → 0 , it has been confirmed that the discharges and discharge points from numerical calculation and empirical methods (Schaffernak-Van Iterson and L. Casagrande methods) change rapidly. Figure 9 shows the other extreme pattern. The discharge point calculated by the Schaffernak-Van Iterson method coincides with P * 3 . The free surface distinctly deviates from the numerically calculated free surface FS(n). Nevertheless, its discharge, q SV = 0.250(0.0%), agrees with the theoretical solution, q t = 0.250 (Appendix A). In this method, q SV is equal to the theoretical solution at both extreme states of Figures 6 and 9. Even when the free surface is not correct, if a strict discharge is defined in this method, then a proper result can be obtained by the discharge-equivalent KZ flow, KZ3(C, q SV ). However, in the range α min α 90 • , such a result is not guaranteed. In L. Casagrande's method in Figure 9, q LC = 0.222(−11.4%) and the discharge point separate upward. In general, when paying attention to the free surface, the accuracy is superior to that of the Schaffernak-Van Iterson method, but the discharge error always occurs depending on R b and α (almost the same relative error when α = α min ). A. Casagrande (1937) [18] recommends L. Casagrande's method as an alternative method when his basic parabola method cannot be used. However, it cannot be used from a strict viewpoint.          Tables 1 and 2 summarize the results of the above empirical methods. From the viewpoint of error-within-a-few-percent estimation, none of the empirical methods shown here can be adopted. Therefore, we must establish some new empirical rules.    (16a)), which is almost equal to |AC |. In Figure 10b, when α = 120°, |AC | = 0.154; moreover, it is confirme d that the location of C changes depending on α. On the other hand, by definition, |AC | is a quantity unique to the value of R under a fixed value of θ, so it is the same as in Figure 10a. In Figure 10a, the modification method of FSP(n) is the same as that in Figure 4b. The only difference is that the modified value becomes negative. Tables 1 and 2 summarize the results of the above empirical methods. From the viewpoint of error-within-a-few-percent estimation, none of the empirical methods shown here can be adopted. Therefore, we must establish some new empirical rules.    Figure 10 shows the calculation examples of R b = 1.00. In Figure 10a, when α = 60 • , AC n = 0.136. The starting point of the parabola defined by the basic parabola of the basic KZ flow is expressed as AC = 0.135 (using Equation (16a)), which is almost equal to AC n . In Figure 10b, when α = 120 • , AC n = 0.154; moreover, it is confirmed that the location of C n changes depending on α. On the other hand, by definition, AC is a quantity unique to the value of R b under a fixed value of θ, so it is the same as in Figure 10a. In Figure 10a, the modification method of FSP(n) is the same as that in Figure 4b. The only difference is that the modified value becomes negative. Using the system adopted in Figure 10, and are systematically changed to develop the tables for , , and . The results are listed in Tables 3-5. Figure 11a-c show the range of ≤ 3.50 listed in the Tables. In Figure 11a,b, the theoretical values in the case of = are also shown. In Figure 11c, each value of − | | is shown by a red solid line. When = 1.00 and ≥ 60°, it is found that − | | ≅ − | |. When ≥1.25, this practically holds true regardless of α. In Figure 11c, it is noted that if | | is adopted instead of −| |, the above simple relationship is not found. Parts of the parenthesized value (#. ###) listed in Tables 3-5 are those that cannot be numerically calculated under the fixed value of ∆ = ∆ = 0.05 (relatively it is too large, namely, the calculation element (near-wall element) cannot be defined on the right side of point O, and calculation becomes impossible). Usually, these values may be unnecessary; however, they can be estimated from the values at α = 90°. In Tables 3 and 5, each value would be regarded as the value at α = 90°.

Creating Basic Tables
On the other hand, in Table 4, each value (#. ###) is linearly decreasing from the value at α = 90° to 0 at α = 180° (Figure 11b). Using the system adopted in Figure 10, R b and α are systematically changed to develop the tables for q n , yB n , and C n . The results are listed in Tables 3-5. Figure 11a-c show the range of R b ≤ 3.50 listed in the Tables. In Figure 11a,b, the theoretical values in the case of α = α min are also shown. In Figure 11c, each value of R b − AC is shown by a red solid line. When R b = 1.00 and α ≥ 60 • , it is found that R b − AC R b − AC n . When R b ≥1.25, this practically holds true regardless of α. In Figure 11c, it is noted that if AC n is adopted instead of − AC n , the above simple relationship is not found. Parts of the parenthesized value (#. ###) listed in Tables 3-5 are those that cannot be numerically calculated under the fixed value of ∆x = ∆y = 0.05 (relatively it is too large, namely, the calculation element (near-wall element) cannot be defined on the right side of point O, and calculation becomes impossible). Usually, these values may be unnecessary; however, they can be estimated from the values at α = 90 • . In Tables 3 and 5, each value would be regarded as the value at α = 90 • . On the other hand, in Table 4, each value (#. ###) is linearly decreasing from the value at α = 90 • to 0 at α = 180 • (Figure 11b).   There is a slight discrepancy between FS(n) and modified FS(n). As previously described, Equation (18) is identified by the rectangular dam and is calculated as L * h = 0.42 · · · . If L * h = 0.20 is specified in Equation (17) as an individual value in this case, both of them visually match. However, in this paper, we consider this degree of deviation as an acceptable deviation.  Tables 3-5, the discharge, discharge point, and starting point of the FSP corresponding to any two conditions, and α , can be estimated using interpolation methods. There are several two-dimensional mathematical interpolation methods-the bilinear and bicubic interpolation methods (Web-Interpolation, 2017) may be fundamental. The details are described in Appendix C. In the interpolation in Tables 3-5, the bicubic interpolation method is used as a whole; however, in the range of 90°< α < 180°, the bilinear interpolation method would give sufficient accuracy.

From the information listed in
By interpolation, three estimation values are obtained: discharge, , y-coordinate value of discharge point, , and starting point, . The discharge point, , is defined from and its discharge angle . The following five methods are conceivable as the interpolation-equivalent KZ flow methods. Under using point C, (i) KZI( , ), (ii) KZI( , ) are thought. Estimation values are , , and , then the three kinds of KZ flow can be obtained: (iii) KZI( , ), (iv) KZI( , ), and (v) KZI( , ). In the interpolation values, the value is a key factor; therefore, in the above combinations, it is conceivable to use (i), (iii), and (iv). Here, we will confirm the results of (i) KZI( , ) and (iii) KZI( , ) directly for a brief description. Figure 12a shows the results of KZn( , ) when = 1.125 and = 65°. As for KZI( , ) and KZI( , ), only the calculated results are shown in the inset table of the figure to avoid the complexity of the illustration; they are within an error of 1%. Figure 12b shows the results of KZn( , ) when = 1.125 and = 135°. The numerical results of KZI( , ) are within an error of 1%, but in the calculation of KZI( , ), the calculated discharge is * = 0.422(−1.5%). In the range ≥ 1.00, all the calculations of KZE( , ) almost become the error-within-a-few-percent estimation. This is because in the range, point C practically almost agrees with point (see Figure  11c). On the other hand, in the range ≤ 1.00 KZI estimation accuracy becomes poor. The examples of Figure 13 with = 0.625 show this fact. In Figure 13a, = 75° and so (≡ | ̅̅̅̅ |) = 0.357. When is small and as → 0, both discharge and discharge point change rapidly. Therefore, the accuracy of interpolation estimation deteriorates. Figure 13b is the example of = 135°. The result of KZI( , ) shows good accuracy. However, the calculated value of the discharge of KZI( , ) is * = 0.673(−6.1%). In conclusion, in order to secure the accuracy in the case of ≤ 1.00, it is necessary to set the calculation points densely; such examples are described later.  Tables 3-5, the discharge, discharge point, and starting point of the FSP corresponding to any two conditions, R b and α, can be estimated using interpolation methods. There are several two-dimensional mathematical interpolation methods-the bilinear and bicubic interpolation methods (Web-Interpolation, 2017) may be fundamental. The details are described in Appendix C. In the interpolation in Tables 3-5, the bicubic interpolation method is used as a whole; however, in the range of 90 • < α < 180 • , the bilinear interpolation method would give sufficient accuracy.

From the information listed in
By interpolation, three estimation values are obtained: discharge, q e , y-coordinate value of discharge point, yB e , and starting point, C e . The discharge point, B e , is defined from yB e and its discharge angle α. The following five methods are conceivable as the interpolation-equivalent KZ flow methods. Under using point C, (i) KZI(C, B e ), (ii) KZI(C, q e ) are thought. Estimation values are q e , yB e , and C e , then the three kinds of KZ flow can be obtained: (iii) KZI(q e , B e ), (iv) KZI(C e , B e ), and (v) KZI(C e , q e ). In the interpolation values, the value yB e is a key factor; therefore, in the above combinations, it is conceivable to use (i), (iii), and (iv). Here, we will confirm the results of (i) KZI(C, B e ) and (iii) KZI(q e , B e ) directly for a brief description. Figure 12a shows the results of KZn(q n , B n ) when R b = 1.125 and α = 65 • . As for KZI(q e , B e ) and KZI(C, B e ), only the calculated results are shown in the inset table of the figure to avoid the complexity of the illustration; they are within an error of 1%. Figure 12b shows the results of KZn(q n , B n ) when R b = 1.125 and α = 135 • . The numerical results of KZI(q e , B e ) are within an error of 1%, but in the calculation of KZI(C, B e ), the calculated discharge is q * e = 0.422(−1.5%). In the range R b ≥ 1.00, all the calculations of KZE(C, B e ) almost become the error-within-a-few-percent estimation. This is because in the range, point C practically almost agrees with point C n (see Figure 11c). On the other hand, in the range R b ≤ 1.00 KZI estimation accuracy becomes poor. The examples of Figure 13 with R b = 0.625 show this fact. In Figure 13a, α = 75 • and so R u ≡ AG = 0.357. When R b is small and as R u → 0 , both discharge and discharge point change rapidly. Therefore, the accuracy of interpolation estimation deteriorates. Figure 13b is the example of α = 135 • . The result of KZI(q e , B e ) shows good accuracy. However, the calculated value of the discharge of KZI(C, B e ) is q * e = 0.673(−6.1%). In conclusion, in order to secure the accuracy in the case of R b ≤ 1.00, it is necessary to set the calculation points densely; such examples are described later.

Application Examples of Basic Tables to General Dams
Here, a usual trapezoidal dam without a drain is considered. Figure 14a shows the numericequivalent KZ flow, KZn( , ) . The method of setting the equivalent vertical face, ′ ̅̅̅̅̅̅ , is explained in Figure 5c. It is assumed that the trapezoidal cross section, ′ 3 * , yields the same discharge and discharge point. However, the vertical face has to be estimated without using the numeric-equivalent KZ flow. This is expressed by using the interpolation-equivalent KZ flow, KZI( , ). To obtain by the functional definition, Equation (20), it is necessary to create a threedimensional table with the parameters , , and ; however, this is not a promising idea because it requires many calculation points. Let be estimated by a two-dimensional interpolation using Table 4. For this purpose, a trapezoidal cross section with a vertical entrance face corresponding to this seepage phenomenon has to be found.
In Figure 14b, first, point C of the basic parabola is located. Thereafter, ( )(≈ ( )) and the position of (≈ ) are determined, and the corresponding equivalent vertical face, ′ ̅̅̅̅̅̅̅ , is

Application Examples of Basic Tables to General Dams
Here, a usual trapezoidal dam without a drain is considered. Figure 14a shows the numeric-equivalent KZ flow, KZn(q n , B n ). The method of setting the equivalent vertical face, C e C e , is explained in Figure 5c. It is assumed that the trapezoidal cross section, C E C E P * 3 G, yields the same discharge and discharge point. However, the vertical face has to be estimated without using the numeric-equivalent KZ flow. This is expressed by using the interpolation-equivalent KZ flow, KZI(C, B e ). To obtain yB e by the functional definition, Equation (20), it is necessary to create a three-dimensional table with the parameters θ, R b , and α; however, this is not a promising idea because it requires many calculation points. Let yB e be estimated by a two-dimensional interpolation using Table 4. For this purpose, a trapezoidal cross section with a vertical entrance face corresponding to this seepage phenomenon has to be found.
In Figure 14b, first, point C of the basic parabola is located. Thereafter, BPP(b)(≈ BPP(e)) and the position of G b (≈ G e ) are determined, and the corresponding equivalent vertical face, C E C E , is defined (see Figure 5c). According to this definition, the equivalent aspect ratio becomes R E = C E P * 3 /H = 2.738. With R b = R E , yB e is determined by finding yB e /H by two-dimensional interpolation from the list in Table 4; FSP(e) is determined from KZI(C, B e ). The numerical calculation result is q n = 1.71 ± 0.01 and yB n = 3.50, and KZI(C, B e ) yields q * e = 1.71(−0.01%) and yB e = 3.58(0.9%). The asterisk of q * e means that the value is calculated from KZI(C, B e ), not directly the interpolated value. The discharge of the cross section, C E C E P * 3 G, is interpolated as q E = 1.72(0.5%). In this example, under the definition of P * 3 (45.00, 0.00)(α = 20.3 • ), it yields q n = 1.18 ± 0.01, yB n = 3.82, and q * e = 1.18(−0.1%), yB e = 3.79(−0.3%); these results have sufficient accuracy.
Eng 2020, 1, FOR PEER REVIEW 23 of 39 defined (see Figure 5c). According to this definition, the equivalent aspect ratio becomes = | ′ 3 * ̅̅̅̅̅̅ |⁄ = 2.738 . With = , is determined by finding ⁄ by two-dimensional interpolation from the list in Table 4; ( ) is determined from KZI( , ) . The numerical calculation result is = 1.71 ± 0.01 and = 3.50 , and KZI( , ) yields * = 1.71(−0.01%) and = 3.58(0.9%). The asterisk of * means that the value is calculated from KZI( , ), not directly the interpolated value. The discharge of the cross section,    Figure 15 shows an example of another interpolation-equivalent KZ flow, KZI(C, B e ), which corresponds to Figure 5. A precise conclusion, q * e = 2.01(0.4%) and yB e = 2.82(−0.1%), is obtained. Many other cases were examined for KZI(C, B e ). Concerning the discharge angle, the applicable range is 20 • ≤ α < 180 • , and regarding the entrance angle, the applicable range is 20 • ≤ θ ≤ 90 • . However, the equivalent aspect ratio, R E , must be smaller than 4.5 (Tables 3-5).  Many other cases were examined for KZI( , ). Concerning the discharge angle, the applicable range is 20 ≤ < 180 , and regarding the entrance angle, the applicable range is 20 ≤ ≤ 90 . However, the equivalent aspect ratio, , must be smaller than 4.5 (Tables 3-5).   By the above formulation, it seems that the target of error-within-a-few-percent estimation has practically been reached. However, the combination of θ and α where R b and R u (R b ≥ R u ) are exceedingly small must be eliminated.

Special Cases When the Basic Table Cannot Be Used
When the up-side aspect ratio, R u , is small in the trapezoidal dam with a vertical entrance face, the calculation accuracy is poor, as shown in Figure 13a. To improve this, a table for q n and yB n is developed by systematically changing R u and R b ; the results are listed in Tables 6 and 7. Because the defined interval is small, bilinear interpolation can be used (Appendix C). The tables are utilized for analyzing the case of R b ≥ R u ; however, for the continuous analysis with R u ≥ 0.20, the values in parentheses (R b < R u ) are also prepared by numerical calculation. Figure 16 shows a calculation example. In the defined region of the tables, KZI(q e , B e ) ensures that the error-within-one-percent estimation is achieved. Table 6. Trapezoidal dam with θ = 90 • , and small value of R u , discharge q n depending on R u and R b . By the above formulation, it seems that the target of error-within-a-few-percent estimation has practically been reached. However, the combination of and where and ( ≥ ) are exceedingly small must be eliminated.

Special Cases When the Basic Table Cannot Be Used
When the up-side aspect ratio, , is small in the trapezoidal dam with a vertical entrance face, the calculation accuracy is poor, as shown in Figure 13a. To improve this, a table for and is developed by systematically changing and ; the results are listed in Tables 6 and 7. Because the defined interval is small, bilinear interpolation can be used (Appendix C). The tables are utilized for analyzing the case of ≥ ; however, for the continuous analysis with ≥ 0.20, the values in parentheses ( < ) are also prepared by numerical calculation. Figure 16 shows a calculation example. In the defined region of the tables, KZI( , ) ensures that the error-within-one-percent estimation is achieved.  Next, when focusing on and analyzing the center core part of an earth dam, dedicated tables are also prepared. The cross section is bilaterally symmetric. The results are summarized in Tables 8 and 9. In this case, the down-side aspect ratio is defined as ≡ * * . The table is Next, when focusing on and analyzing the center core part of an earth dam, dedicated tables are also prepared. The cross section is bilaterally symmetric. The results are summarized in Tables 8 and 9. In this case, the down-side aspect ratio is defined as R d ≡ P * 2 P * 3 . The table is intended for analyzing the result in the case of R d ≥ R u ; however, for the continuous analysis with R u ≥ 0.20, the values in parentheses (R d < R u ) are also prepared. Figure 17 shows a calculation example. The bilinear interpolation is used as well. Using the defined region of the tables, KZI(q e , B e ) almost ensures that the error-within-one-percent estimation is achieved. intended for analyzing the result in the case of ≥ ; however, for the continuous analysis with ≥ 0.20, the values in parentheses ( < ) are also prepared. Figure 17 shows a calculation example. The bilinear interpolation is used as well. Using the defined region of the tables, KZI( , ) almost ensures that the error-within-one-percent estimation is achieved.  Table 9. dam center core, discharge point depending on and . There may be other particular cases to be investigated. For example, the problems pertaining to inclined dam cores and dams with an entrance angle of > 90 have not yet been considered. However, in these cases, it has been confirmed that the numeric-equivalent KZ flow, KZn ( , ), is effective as well. In such cases, if necessary, two-dimensional tables for obtaining the interpolationequivalent KZ flow could also be prepared.

Theoretical Consideration on the Seepage of Rectangular Dams
Finally, we will consider a rectangular dam where = . Figure 18a shows how of the rectangular dam changes depending on . The regression Equation (21) The theoretical discharge is = (2 ) ⁄ = 1 (2 ) ⁄ [19]. Let the discharge point defined by Equation (21) be . Then, the equivalent KZ flow is defined by ( , ). Figure 19a shows that the result of the numerical calculation and the result of ( , ) agree with high There may be other particular cases to be investigated. For example, the problems pertaining to inclined dam cores and dams with an entrance angle of θ > 90 • have not yet been considered. However, in these cases, it has been confirmed that the numeric-equivalent KZ flow, KZn (q n , B n ), is effective as well. In such cases, if necessary, two-dimensional tables for obtaining the interpolation-equivalent KZ flow could also be prepared.

Theoretical Consideration on the Seepage of Rectangular Dams
Finally, we will consider a rectangular dam where R b = R u . Figure 18a shows how yB of the rectangular dam changes depending on R b . The regression Equation (21) is as follows [5]: Eng 2020, 1

83
The theoretical discharge is q t = H 2 /(2d) = 1/(2R b ) [19]. Let the discharge point defined by Equation (21) be B f . Then, the equivalent KZ flow is defined by KZ f q t , B f . Figure 19a shows that the result of the numerical calculation and the result of KZ f q t , B f agree with high accuracy.
Modified FS(f ) calculated from FSP(f ) by Equations (17) and (18) also almost agrees with the free surface FS (n) obtained by numerical calculation.  (17) and (18) also almost agrees with the free surface FS (n) obtained by numerical calculation. It is clear that Equation (21) is a quadratic equation and has an application limit. In the case of > 1.4, the calculation is conducted by using the interpolation-equivalent KZ flow, KZI( , ), based on Table 4; however, there is a convenient method specifically used for rectangular dams. Figure 20a shows the relationship between and in the same data as that shown in Figure 18a. It is certain that as → ∞ , ( , ) → (0,0) . As  Here, the condition that the seepage phenomenon in the vicinity of the discharge point becomes independent of the shape of the entrance face is considered. As shown in Figure 4b, the theoretical solution of the basic KZ flow, KZb( , ), and that of the numerical calculation are different as a whole; however, near the discharge point, the two are practically the same. This indicates that when is increased to some extent, the flow near the discharge point is independent of the shape of the entrance face. This is an example of the case α = 180; nevertheless, the concept seems to hold for any value of α. Let this concept be called "discharge-independent principle". Figure 20a implies that the above concept holds even at α = 90. Theoretically, it is presumed that Equation (23) is the asymptotic line of the function = ( ) as → ∞( → 0). It may be considered that a discharge point specific to the discharge is defined in the range ≥ 1.0. It is clear that Equation (21) is a quadratic equation and has an application limit. In the case of R b > 1.4, the calculation is conducted by using the interpolation-equivalent KZ flow, KZI(C, B e ), based on Table 4; however, there is a convenient method specifically used for rectangular dams. Figure 20a shows the relationship between yB and q in the same data as that shown in Figure 18a. It is certain that as R b → ∞ , (q, yB) → (0, 0) . As R b = 1.4, (q, yB) = (0.357, 0.249). The linear equation passing through these two points is as follows.
This seems to hold in the range 0.8 ≤ R b . Now, q = q t = H 2 /(2d) = 1/(2R b ), therefore This equation is shown in Figure 18a, (II). Although the ranges of Equations (21) and (23) overlap, let us set the application range of Equation (21) as R b < 1.0 and that of Equation (23) as R b ≥ 1.0 for convenience.
Here, the condition that the seepage phenomenon in the vicinity of the discharge point becomes independent of the shape of the entrance face is considered. As shown in Figure 4b, the theoretical solution of the basic KZ flow, KZb(C, O b ), and that of the numerical calculation are different as a whole; however, near the discharge point, the two are practically the same. This indicates that when R b is increased to some extent, the flow near the discharge point is independent of the shape of the entrance face. This is an example of the case α = 180 • ; nevertheless, the concept seems to hold for any value of α. Let this concept be called "discharge-independent principle". Figure 20a implies that the above concept holds even at α = 90 • . Theoretically, it is presumed that Equation (23) is the asymptotic line of the function yB = f (q) as R b → ∞(q → 0) . It may be considered that a discharge point specific to the discharge is defined in the range R b ≥ 1.0. The accuracy of numerical calculation depends on the difference width. In this paper, as H = 1.00 is designated, ∆x = ∆y = 0.05 is used for all numerical calculations. As using ∆ = ∆ = 0.025, the result is = 1.429(0.0%) + 0.004, and = 0.742. There is almost no difference in the discharge compared to the former, but the position of the discharge point is calculated slightly higher. The theoretical solution for the free surface of a rectangular dam was derived by Polubarinova-Kochina [19]. In quantifying the theoretical solution, a systematic iterative calculation is required, which was detailed in the previous paper [4]. Figure 19b shows the comparison between the numerical solution and the theoretical solution. It is observed that the agreement is almost good. It is thought that all of the results in Figure 18a will be slightly different if ∆x = ∆y = 0.025 is used. Note that the theoretical solution only gives the free surface, but it requires much more calculation time than the numerical calculation of IFDM. Finally, let us consider the case of 90 < α < 180. Figure 21 shows the relationship between and depending on α : = ( ) obtained from Tables 3 and 4. It is shown in the range 0.5 ≤ ≤ 3.5. We have already described the case of α = 90, here the line of Equation (22), = ( = 0.7), is also shown. Corresponding to 90 < α < 180, it is thought to be 0.7 > β > 0. At α = The accuracy of numerical calculation depends on the difference width. In this paper, as H = 1.00 is designated, ∆x = ∆y = 0.05 is used for all numerical calculations. As using ∆x = ∆y = 0.025, the result is q n = 1.429(0.0%) + 0.004, and yB n = 0.742. There is almost no difference in the discharge compared to the former, but the position of the discharge point is calculated slightly higher. The theoretical solution for the free surface of a rectangular dam was derived by Polubarinova-Kochina [19]. In quantifying the theoretical solution, a systematic iterative calculation is required, which was detailed in the previous paper [4]. Figure 19b shows the comparison between the numerical solution and the theoretical solution. It is observed that the agreement is almost good. It is thought that all of the results in Figure 18a will be slightly different if ∆x = ∆y = 0.025 is used. Note that the theoretical solution only gives the free surface, but it requires much more calculation time than the numerical calculation of IFDM.
Finally, let us consider the case of 90 • < α < 180 • . Figure 21 shows the relationship between yB and q depending on α: yB = f α (q) obtained from Tables 3 and 4. It is shown in the range 0.5 ≤ R b ≤ 3.5. We have already described the case of α = 90 • , here the line of Equation (22), yB = βq (β = 0.7), is also shown. Corresponding to 90 • < α < 180 • , it is thought to be 0.7 > β > 0. At α = 180 • , yB = 0, Eng 2020, 1 86 but at discharge point (xB, yB), xB = 0.5q, so the discharge point is uniquely determined, see Figure 1. Now, suppose the following relationship holds. 180 , = 0 , but at discharge point ( , ) , = 0.5 , so the discharge point is uniquely determined, see Figure 1. Now, suppose the following relationship holds. In order to confirm the numerical calculation accuracy, we will examine how the above discussion regarding the rectangular dam changes when theoretical solutions are used for both q and yB. Figure 18a is changed to Figure 18b. Figure 20a is changed to Figure 20b. There is no change in the theoretical formation, but there is a 1-2% difference in the calculation of yB. In the numerical calculation, the discharge calculation error remains within 1%, may be around 0.1%. However, the error of yB changes depending on the calculation conditions, especially the difference width. Higheraccuracy numerical calculation is possible if necessary, but it is considered that the numerical calculations performed in this paper have sufficient accuracy from an engineering point of view. In Figure 21, the line = is shown regarding α = 120 and α = 150. Although there is a slight discrepancy in numerical calculations, Equation (24) seems to be valid. This indicates that under the conditions of 90 < α < 180 and ≥ 1.0, the discharge point is uniquely determined corresponding to and α. In the range α < 90, it is easily confirmed that Equation (24) does not hold. In such a case, the discharge independent principle still may hold. This is a problem to be solved in the future.

Summary of This Section
This section details the new empirical method. In order to provide the benefits of applying this method to the field, the whole picture is summarized in a compact manner. As expressed by Equations (19) and (20), the shape determinant of the isotropic seepage phenomenon of earth dam is ( , , , ). However, by applying the similarity law of seepage phenomenon, it is reduced to (1, , , ). Furthermore, by replacing the dam of arbitrary shape with an equivalent verticalentrance-face trapezoidal dam, it is reduced to (1,90, , ). The basic process is as follows: (i) To replace a general dam with an equivalent vertical-entrance-face trapezoidal dam, determine the starting point C of the basic parabola by Equation (16), (ii) Furthermore, the position of the vertical entrance face is determined in consideration of Equation (13) (see Figures 14 and 15). Now we are In order to confirm the numerical calculation accuracy, we will examine how the above discussion regarding the rectangular dam changes when theoretical solutions are used for both q and yB. Figure 18a is changed to Figure 18b. Figure 20a is changed to Figure 20b. There is no change in the theoretical formation, but there is a 1-2% difference in the calculation of yB. In the numerical calculation, the discharge calculation error remains within 1%, may be around 0.1%. However, the error of yB changes depending on the calculation conditions, especially the difference width. Higher-accuracy numerical calculation is possible if necessary, but it is considered that the numerical calculations performed in this paper have sufficient accuracy from an engineering point of view.
In Figure 21, the line yB = βq is shown regarding α = 120 • and α = 150 • . Although there is a slight discrepancy in numerical calculations, Equation (24) seems to be valid. This indicates that under the conditions of 90 • < α < 180 • and R b ≥ 1.0, the discharge point is uniquely determined corresponding to q and α. In the range α < 90 • , it is easily confirmed that Equation (24) does not hold. In such a case, the discharge independent principle still may hold. This is a problem to be solved in the future.

Summary of This Section
This section details the new empirical method. In order to provide the benefits of applying this method to the field, the whole picture is summarized in a compact manner. As expressed by Equations (19) and (20), the shape determinant of the isotropic seepage phenomenon of earth dam is (H, θ, R b , α). However, by applying the similarity law of seepage phenomenon, it is reduced to (1, θ, R b , α). Furthermore, by replacing the dam of arbitrary shape with an equivalent vertical-entrance-face trapezoidal dam, it is reduced to (1, 90 • , R b , α). The basic process is as follows: (i) To replace a general dam with an equivalent vertical-entrance-face trapezoidal dam, determine the starting point C of the basic parabola by Equation (16), (ii) Furthermore, the position of the vertical entrance face is determined in consideration of Equation (13) (see Figures 14 and 15). Now we are ready to use Basic Table 4, (iii) The discharge point B e is determined by Table 4 owing to applying the similarity law, then the empirical QF estimation is completed by KZI(C, B e ).
In order to execute the above calculation precisely, it is necessary to create a dedicated system. However, a rougher approach would be acceptable to approximately know the discharge and the discharge point. In this paper, the bicubic interpolation method is used in the interpolation of the basic table, but the results are not so different even in the bilinear interpolation.
This empirical method using the basic Table 4 has a wide applicable range: R b > 1.0, 20 • ≤ α < 180 • , and 20 • ≤ θ ≤ 90 • . Regarding the three exceptional cases where sufficient accuracy is not ensured by the above method, (i) vertical-entrance-face trapezoidal dam with small upper side, (ii) symmetrical dam center core, and (iii) narrow rectangular dam, dedicated tables were created individually. By creating Tables 6-9 and the insert table of Figure 18, the approximate values of discharge and discharge point can be instantaneously grasped.
Interestingly, the actual soil material is generally anisotropic. In this case, the governing equation of the seepage phenomenon is often expressed by the following equation.
k x is the permeability coefficient in the x direction and k y is the permeability coefficient in the y direction.
In this case, the flow net is created geometrically or by numerical calculation from the flow net created assuming isotropy (Harr, 1991, p. 29 [20]). For this reason, the empirical method proposed in this paper is also applied to such anisotropic seepage dams.

Conclusions and Discussion
In the previous paper [8], the concept of the equivalent KZ flow method was formulated. In this study, this method is developed as a more accurate empirical method. The contents of this paper are summarized as follows.
Conventionally, there were three empirical methods stipulated in the design standards regarding the seepage of earth dam: (i) A. Casagrande's method, (ii) the Schaffernak-Van Iterson method, and (iii) L. Casagrande's method. Theoretical consideration on the method (i) and the evaluation of the error were shown in the previous paper. In this paper, the methods (ii) and (iii) are newly considered and the characteristics of the calculation error are confirmed. The investigations revealed that it is impossible to stably use these methods even if the permissible error is set to 10%. It is clear that the conventional methods essentially need to be reviewed in order to reach the accuracy within a few-percent error in a wide range.
For the natural development of the logic from the previous paper [5] to the present paper, the basic conclusions of the concept of the equivalent KZ flow method presented in the previous paper were briefly summarized. The area-equivalent KZ flow, KZ1(C, A e ), is an important method reached in the previous paper. This method powerfully rebuts the basic parabola method of A. Casagrande. This method obtains good results in the range of 90 • ≤ α ≤ 180 • , but in the range of α < 90 • , there is a drawback that the accuracy becomes worse as α becomes smaller. In order to overcome this difficulty, the concept of the numeric-equivalent KZ flow, KZn(q n , B n ), is important among the several equivalent KZ flows proposed in the previous paper. By using this KZ flow, the free surface shape that fits the numerical calculation can be calculated almost accurately in any shape of trapezoidal dam. A. Casagrande's method estimates the discharge and discharge point, but no method for accurately estimating the free surface has been proposed. This is defined as empirical QB estimation. On the other hand, in the numeric-equivalent KZ flow, the free surface itself can be estimated by modifying the corresponding free surface parabola FSP(n) using a modification equation. This is defined as empirical QF estimation. All of the equivalent KZ flow methods in this paper give the empirical QF estimation. A set of the determinant factors of the seepage of a trapezoidal dam is (H, θ, R b , α), but it became clear that only the case of H = 1.00 should be considered due to the similarity law of seepage. Therefore, in order to determine the values of q and yB by interpolation method, it is sufficient to create a three-dimensional " Table" of (θ, R b , α). It is possible in principle but inconvenient. In this paper, a two-dimensional table is created by fixing θ = 90 • . As long as θ = 90 • , the interpolation is completed, and if the estimated values are q e and yB e , the corresponding interpolation-equivalent KZ flow, KZI(q e , B e ), is obtained. However, in the case of θ < 90 • , an empirical method must be devised separately. A. Casagrande proposed a method to determine the starting point of the basic parabola C a . The starting point changes depending on all the factors of θ, R b , and α, but he defines the basic parabola drawn when α = 180 • and even if α changed, he considered the starting point immutable. The author basically adopted his hypothesis (in this sense, the author's idea is basically dependent on A. Casagrande's method) but with a more rational method for determining the starting point C of the basic parabola.
For the trapezoidal dam with vertical entrance face, in the range R b ≥ 1.00, the starting point C can be applied as a fixed value that does not depend on α (detailed Table 5, Figure 11c). This is extended to a usual homogeneous dam, and highly accurate discharge and free surface are calculated (Figures 14 and 15). That is, an equivalent trapezoidal dam with vertical entrance face, which is almost equivalent to the discharge and the discharge point of the dam itself, is defined, and the discharge point yB e is estimated from Table 4 based on the similarity law. Finally, interpolation-equivalent KZ flow, KZI(C, B e ), calculates the discharge and free surface. It has been confirmed that KZI(C, B e ) almost has an error-within-one-percent estimation in a wide range of 20 • ≤ θ ≤ 90 • and 20 • ≤ α ≤ 180 • .
Even in the case of the trapezoidal dam with vertical entrance face, if R u ≤ R b < 1.00 and as R u → 0 , sufficient accuracy cannot be expected in both KZI(q e , B e ) and KZI(C, B e ). In this case, a special method needs to be adopted. In general, the three parameters q n , B n , and C n rapidly change in the process such that the up-side aspect ratio of the trapezoidal dam becomes R u → 0; therefore, it is necessary to consider focusing on R u . The following two patterns are investigated: (i) For a trapezoidal dam with a vertical entrance face, the equivalent KZ flow is KZI(q e , B e ); the applicable ranges are R u ≤ R b , 0.20 ≤ R u ≤ 0.50, and 0.20 ≤ R b ≤ 1.00. (ii) For the symmetrical dam center core, the equivalent KZ flow is KZI(q e , B e ); the applicable ranges are R u ≤ R d , 0.20 ≤ R u ≤ 0.50, and 0.20 ≤ R d ≤ 1.40.
With the above formulations, the calculation which is practically necessary can be almost covered. In addition to this, it is necessary to consider separately the cases of inclined core or θ > 90 • , but KZn(q n , B n ) is also effective in such cases. If the corresponding two-dimensional table is created, it will be possible to perform highly accurate calculation by the same method as that in (1) and (2) above.
Finally, we examined the case of R u = R b , that is, a rectangular dam. This is already included as a special case in (i) and (ii), but there is a problem to be considered independently. For a rectangular dam, the equivalent KZ flow is KZ f q t , B f ; the applicable range is 0.2 ≤ R b ≤ 1.4; however, it is finally extended to 0.2 ≤ R b → ∞ according to the concept of "discharge-independent principle". That is, when the basic aspect ratio, R b , increases to some extent, the discharge phenomenon in the vicinity of the discharge point becomes a phenomenon peculiar to the discharge amount itself, and has no relationship with the shape of the entrance face. This concept is extended to the case of arbitrary α. By measuring yB, the discharge may be determined by q = (k c )βyB, in which β depends on α.
The accuracy of the numerical calculation used in this paper was confirmed in the previous paper. The theoretical solution for generally obtaining the discharge and discharge point of a trapezoidal dam with vertical entrance face is in the literature [19]), but it is exceedingly difficult. Lo (1971) [9] numerically evaluated them (discharge and discharge point) and showed this graphically. Based on his graphs, the author confirmed several cases of seepage and it was judged that they were in agreement with the author's numerical calculation results within some error, but there was a reading error and it could not withstand rigorous evaluation. In relation to this, the graphic display that enables various estimations was valuable in an era when the calculation method was insufficient, and it was also useful for intuitively grasping the whole image. However, now, it is considered that the numeric information that makes a pair with such a graph is essential, because the computer can instantly perform calculation processing using an interpolation method.
In this paper, the validity of the numerical calculation in the correspondence with the theoretical solution is directly confirmed only for the KZ flow and the rectangular dam flow. Regarding the latter, the numerical calculation result and the theoretical solution agree with the discharge, but there may be an error of about 1-2% in the discharge point. From the viewpoint of the goal of this paper, within a few percent of error, it was judged that there is no problem. Now, let us point out the following. In the calculation of the free surface, there is considerable freedom in setting the initial surface, but many iterative calculations are required until it converges. If the high-accuracy initial surface is set, the calculation is completed by repeating the calculation several times. Theories evolve in stages, from simple to complex, from special to general. As pointed out in Section 5.5, soil materials generally exhibit anisotropy. The general governing equation for anisotropic seepage problems including the case of treating full permeability tensor has already been established (ex. Lei et al. 2015 [21]), and numerical calculation by IFDM-BPIM is possible. Furthermore, a discontinuous line (or plane) of the soil material is formed in the calculation domain. A general numerical calculation system in this case will be built based on the findings confirmed in the previous two papers [4,5] and this paper.
Funding: This research received no external funding.

Conflicts of Interest:
The author declares no conflict of interest.   basic aspect ratio=H/d R E equivalent aspect ratio α discharge angle α b acceleration factor of the TMSD scheme α bmax theoretical maximum acceleration factor of the TMSD scheme (=2.00) α bopt optimum acceleration factor ∆t time difference width ∆t c criteria time difference ∆x

Abbreviations
x-direction difference width ∆y y-direction difference width θ entrance angle θ c criteria entrance angle ϕ i, j time-interval adjustment factor