Symmetrical Rank-Three Vectorized Loading Scores Quasi-Newton for Identification of Hydrogeological Parameters and Spatiotemporal Recharges

: In a multi-layered groundwater model, achieving accurate spatiotemporal identification and solving the ill-posed problem is the vital topic for model calibration. This study proposes a symmetry rank three vectorized loading scores (SR3 VLS) quasi-Newton algorithm by modifying the Levenberg–Marquardt algorithm and importing a rank three structure from Broyden–Fletcher–Goldfarb–Shanno algorithm for identification of hydrogeological parameters and spatiotemporal recharge simultaneously. To accelerate directional convergence and approach a global optimum, this study uses a vectorized limited switchable step size in the transmissive groundwater inverse problem. The Hessian approximation rank three uses high and low-rank factor loading scores analyzed from simulated storage fluctuation between adjacent iterations for calculation and matrix correction. Two numerical experiments were designed to validate the proposing algorithm, showing the SR3 VLS quasi-Newton reduced the error percentages of the identified parameters by 1.63% and 9.65% compared to the Jacobian quasi-Newton. The proposing method is applied to the Chou-Shui River alluvial fan groundwater system in Taiwan. Results show that the simulated storage error decreased rapidly in six iterations, and has good head convergence as small as 0.11% with a root-mean-square-error (RMSE) of 0.134 m, indicating that the proposing algorithm reduces the computational cost to converge to the true solution. research topic in model calibration. This study proposes an SR3 VLS quasi-Newton algorithm by modifying the Levenberg–Marquardt algorithm and importing a rank three structure from the Broyden–Fletcher–Goldfarb–Shanno algorithm to solve the constrained identification of hydrogeological parameters and spatiotemporal recharge patterns simultaneously associated with an iterative approach in a large-scale multi-layered groundwater model. The recharge parameters include surface and boundary recharge, and the hydrogeological parameters consist of riverbed leakage, aquifer hydraulic, and aquitard leakance conductivities. To accelerate directional convergence, approach a global optimum, and achieve a well-posed problem, this study uses a vectorized limited switchable step size in the inverse problem of transmissive groundwater. Accordingly, to achieve the Hessian approximation, solve nonlinear ill-posed problems, and detect the global minimum area, this study imports a rank three structure, which uses a high-rank and a low-rank factor loading score analyzed from the simulated storage difference between two adjacent iterations to calculate the correction matrix for the Hessian. Moreover, it modifies the traditional algorithm from the scalar line search to vectorized multi-order derivative exact double false position bracketing using SVD-related rank and depth information for solving the vectorized step size corresponding to different hydrogeological parameters and recharges. The proposing methodology is verified by numerical experiments and is applied in a complex practical multi-layered groundwater system. This study designed two numerical experiments to verify the capability of identifying multiple parameters and recharges through the proposing modified algorithm. Experiment 1 identified six zonal hydraulic conductivities, using a single factor loading score to approximate the Hessian. Except for six hydraulic conductivities, the parameters to be identified in Experiment 2 increased six surfaces recharges and six boundaries recharges that uses overlaid high–low factor loading scores. Results show that the iteratively objective function footprint of the three tested quasi-Newton algorithms (Jacobian, LMA, and SR3 VLS) all successfully descend at a quadratic convergence rate to search for an optimum. The identified average error percentage of the Jacobian quasi-Newton, LMA, and the SR3 VLS quasi-Newton


Introduction
In arid and semi-arid areas, if there are insufficient surface water facilities to store runoff, groundwater provides an important water supply source [1]. Conjunctive use planning for surface water and groundwater plays an important role in total water resource management and can be optimized through a sound groundwater management model based on mathematical system analysis. To develop such a model, it is necessary to understand and to be able to estimate the spatiotemporal patterns of pumpage, net recharge, and hydrogeological parameters of the groundwater system, in which the associated simulation model and optimization model are mostly nonlinear. The governing equation of three-dimensional groundwater flow can be expressed as [2] ( ) ( ) ( )

Methodology
This study designs the depictions of the Section Methodology composing of five sub-sections. Section 2.1 shows the proposing optimization model for model calibration in a multi-layer groundwater system. Section 2.2 expresses the proposing SR3 VLS quasi-Newton algorithm by modifying the notable Levenberg-Marquardt algorithm and importing a rank three structure from the BFGS algorithm for solving the constrained identification of hydrogeological parameters and spatiotemporal patterns of recharges simultaneously in a large-scale transmissive groundwater model. Section 2.3 expresses the vectorized multi-order derivative exact double false position bracketing using SVD-related rank and depth modified from the traditional scalar line search. Section 2.4 depicts the procedures of the methodology. Section 2.5 expresses the visualized SVD of the change in storage between k th and (k−1) th iteration for Hessian approximation.

Proposing Optimization Model for Model Calibration in a Multi-Layer Groundwater System
Traditionally, for model calibration purposes, a set of the hydrogeological parameters or source/sink patterns ( η) are optimized from a limited number of observations ( o b s , , t x y h ). If the classical least-squares error (LSE) is used to represent the output simulated groundwater head error, the objective function can be expressed as the minimization problem in Equation (2a) [37]. In this study, the order differs substantially between the storage coefficient in confined aquifers and the specific yield in unconfined aquifers, which are the dominant independent groundwater head fluctuation variables, thus, this study devises the LSE of the simulated and observed storage hydrograph (  sim  obs  sim  obs  , , ,  To formulate the optimization model, the zonation approach [38,39] is used to partition the parameter space according to the observation well location. Assume that the spatiotemporal pumpage pattern can be estimated beforehand, and the storage coefficients and specific yield values can be obtained from pumping tests. Then ( 1) The optimization model contains four constraints: (1) the spatiotemporal pattern of multiple recharges should subject to mass balance, (2) the simulation of groundwater head should obey governing equation (1), (3) the upper and lower bounds of in situ hydrogeological parameters, and (4) the allocated upper and lower bounds of surface recharge. In a real case, the identification of hydrogeological parameters and recharges often form a large-scale highly nonlinear problem.

Modification of Traditional Algorithms in Nonlinear Programming for Proposing SR3 VLS Quasi-Newton Algorithm to Solve Hydrogeological Parameters
In nonlinear programming, Newton's method attempts to find the roots of J' by constructing a sequence k η from the initial guess 0 η that converges towards some value * η satisfying The above-mentioned iterative scheme can be generalized to several dimensions by replacing the derivative with the gradient,   k J  η and the reciprocal of the second derivative with the inverse of the Hessian matrix, The iterative scheme obtained on doing so is as follows: If all second-order partial derivatives of J exist and are continuous over the function domain, then the Hessian matrix Η of J is usually defined as follows: At a local minimum, the Hessian is positive semi-definite. The geometric interpretation of the Gauss-Newton method is that at each iteration,  [17], which avoids slow convergence in the direction of a small gradient. Therefore, the identity matrix was replaced with the diagonal matrix comprising diagonal elements of T D D J J , which is used to solve linear ill-posed problems, as expressed in Equation (8-d To accelerate convergence, jump over local minima, and achieve a well-posed problem, this study uses a vectorized limited switchable step size in the inverse problem of transmissive groundwater whose parameters include both hydrogeology and source simultaneously. Therefore, we achieve In the proposing SR3 VLS quasi-Newton algorithm in Equation are symmetric rank-one matrices; however, their sum is a rank three matrix. The row vectors of the low-rank and high-rank loading score matrix, namely   , , According to this design, this study uses the orthogonal component of loading scores   to correct the direction to better approximate the Hessian. That is, the modified algorithm uses the vector-expanded matrix k θ and k ω to scale each high-rank component of the loading scores according to the orthogonal properties eigenvalue k j  to approximate the Hessian and explore different local minima positions, so there is larger movement along the directions for each parameter where the gradient is smaller. Accordingly, this study designs the symmetry rank three loading scores   to solve nonlinear ill-posed problems, thus, properly handling the approximated Hessian elements by scale and correction at all iterations. The calculation of the scale matrix k θ and k ω are formulated as follows: In large-scale groundwater modeling with a large number of parameters and numerical grids, the transient simulation time increases substantially [41]. To imitate the real complex alluvial groundwater hydrogeology in practice, meticulously detailed parameterization of aquifer's/aquitard's shape is unavoidable. Hence, computing and storing the full Hessian matrix directly may be infeasible, unavailable, and too expensive to compute at all iterations in hydrogeological parameter identification. Furthermore, the convergence of the Jacobian quasi-Newton is not guaranteed. The approximation required to ignore the second-order derivative terms may be valid in two conditions [17], for which convergence is to be expected: the function values m r are small in magnitude, at least around the minimum; and the functions are only slightly nonlinear, so that is relatively small in magnitude. Moreover, when the (non-negative) Marquardt damping factor  (scalar) optimized by a line search changes, the shift vector must be recalculated. Furthermore, in cases with multiple minima, the algorithm converges to the global minimum only if the initial guess is close to the final solution [20]. However, for well-behaved functions and reasonable starting parameters, the LMA tends to be slightly slower than the Jacobian quasi-Newton. Besides, letting    and k  in BFGS. Hence, this study uses the designed preserved approximate Hessians with positive definite and numerically stable selections from the various combinations of low-high rank loading scores along with the vectorized limited switchable step sizes to detect the area of the global minimum and approach the global optimum in a multi-minima problem.

Modification from Traditional Scalar Line Search to Vectorized Multi-Order Derivative Exact Double False Position Bracketing Using SVD-Related Rank and Depth
The traditional double false position provides the exact solution for a root-finding problem, which determines a scalar The fact that false position method always converges, and has advantages that minimize slowdowns, recommends it when speed is needed [22], in which the solved For the vectorized step size root finding problem, which is more difficult and complicated, the double false position can be used and written algebraically in the form: determine at the  -th iteration in the exact optimization procedure of the vectorized step size. If In a conventional iterative nonlinear programming for parameter identification, the information   k J η in Equation (6) where l ζ is the l th unit vector; This study uses symmetry rank three loading score quasi-Newton associated with vectorized step size for solving massive parameters in transmissive groundwater. This study specializes the step size for each parameter with a switchable limited scheme; thus, this approach helps reduce the difference between different kinds of parameters (transmissivity and source/sink) and improves the convergence rate and searching direction by modifying the parameter step size. To achieve small error in the large parameters proposing. This study proposes a systematic strengthening approach of vectorized multi-order derivative exact double false position bracketing using SVD-related rank and depth, which guarantee convergence for solving the vectorized step size corresponding to different kinds of hydrogeological parameters and recharges. The detailed procedure is listed in Step 4 of Section 2.4.

Procedures of the Methodology
Step 1: Set the iteration counter 0 k  of identified parameters and make an initial guess 0 η for the minimum.
Step 1-1: In hypothetical experiment 1, the initial guess of k η is set as a relatively large value.
Step 1-2: In hypothetical experiment 2, the initial guesses of surf Step 1-3: In practical application, the spatiotemporal recharge pattern is estimated by multiplying the temporal pattern calculated by using the lumped storage hydrograph fluctuation and the previously estimated pumpage according to Hsu et al. [16], with the initial guess of the spatial pattern estimated according to the factor loadings calculated from the computed distributed inflow hydrographs (  storage + pumpage). The initial estimate of the hydrogeological parameters is assigned according to the in situ pumping tests.
Step 2: Compute the second-order derivative descent direction through the symmetry rank three storage-based loading scores.
Step 2-1: Simulate the spatiotemporal groundwater head pattern by using the MODFLOW from pumpage, recharge and parameters, boundary conditions, and initial conditions. The modular MODFLOW-2005, which is a computer code that solves the groundwater flow equation, was used to simulate the finite-difference flow through aquifers [42].
Step 2-2: Compute lumped and distributed groundwater storage hydrographs, based on the in situ storage coefficient and specific yield (Sy) estimated from pumping tests, drilled hydrogeological structure, and observed/simulated groundwater head hydrographs.
Step 2-3: Use various multi-rank loading score combinations calculated from the simulated storage difference Step 3: Choose Step 3-1: Set the root-finding iteration counter 0   and make an initial interval ,0 ,0 LB UB , , for the minimum.
Step 3-2: Use a series of potentially faster double false position methods including the proposing vectorized multi-order derivative exact double false position bracketing using SVD-related rank and depth, vectorized Anderson−Björk algorithm, vectorized Illinois algorithm, and vectorized conventional false position method to calculate Step 3-4: If the shrinking rate drops below 0.5, the bisection is used to calculate , MP, k l   .
Step Step 4: Update k k and return to Step 2-1.
Step 5: Until the descent gradient of objective function

SVD of the Change in Storage between k th and (k-1) th Iteration for Hessian Approximation
First, we normalize the data (the average of each storage fluctuation variable equals to 0, and the variance equals to 1) and divide all elements of each row of X by the sample standard deviation of the variable where   s s s s 2  2  sim  sim  1  , , , , , .., SVD can produce a series of datasets consisting of non-correlated new variables (partial simulated storage difference by a parameter factor) for the original analyzed spatiotemporal variables (simulated change in storage between k th and (k−1) th iteration). Considering the SVD of the deviation matrix X: The principal component analysis provides a description model for the simulated spatiotemporal change in storage between the k th and (k−1) th iterations by using SVD. The deviation matrix X can be expressed by the factor score  C and the factor loading F as follows: (17) According to the mathematical meaning, the above formula can be interpreted as: there are L independent parameters caused the change in storage following perturbation by

Experiment 1-Identification of a Single Type of Parameter
The descent footprint of the objective function of Experiment 1 is illustrated in Error! Reference source not found. (a). Results show that the objective function of the three tested algorithms all can successfully descend to search for an optimum value, in which the LMA can descend more rapidly than the SR1 VLS quasi-Newton algorithm and the Jacobian quasi-Newton. However, the SR1 VLS quasi-Newton identifies more accurate parameter values in five of the six zonal parameters, as shown in Error! Reference source not found.. The identified error percentage of the Jacobian quasi-Newton at iteration five across multiple zonal hydraulic conductivity ranges from 0.18% to 3.75% with an average of 2.39%, those of LMA range from 0.15% to 3.50% with an average of 1.78%, and those of the SR1 VLS quasi-Newton range from 0.09% to 2.65% with an average of 0.76%. Moreover, the VLS quasi-Newton algorithm only made 25 number of transient simulation times among five iterations for calculation of the Hessian matrix and step sizes, but that of the Jacobian quasi-Newton and LMA requires 39 number of simulation times for calculation. Hence, the VLS quasi-Newton only needs 63.64% of simulation times comparing to the Jacobian quasi-Newton and LMA in a six-zonal parameter system to achieve and converge to an optimum. The detailed optimization process and discussion of Experiment 1 are described in Section 2 of the Supplementary Materials.

Experiment 2-Identification of Multiple Types of Parameters Simultaneously
The descent footprint of the objective function of Experiment 2 is illustrated in Error! Reference source not found.(b). Results show that the objective function of the three tested algorithms all can successfully descend at a quadratic convergence rate to search for an optimum, although the SR3 VLS quasi-Newton algorithm with partial aid of SR1 VLS structure descends more slowly than the Jacobian quasi-Newton before iteration 3 and descends more rapidly than LMA. However, the SR3 VLS quasi-Newton iteratively can identify more precise values for hydraulic conductivity, surface recharge, and boundary recharge, as shown in Error! Reference source not found. (a), (b), and (c), respectively. As the step size of the LMA and Jacobian quasi-Newton is a scalar, the parameter correction vector tends towards highly sensitive surface recharges. The average identified error percentages of the Jacobian quasi-Newton, LMA and the SR3 VLS quasi-Newton at iteration five across multiple zonal hydraulic conductivities are 6.79%, 53.13%, and 0.04%, respectively; those across multiple surfaces recharges are 1.75%, 5.76%, and 0.09%, and those across multiple boundary recharges are 20.62%, 52.31%, and 0.08%. This demonstrates that the set objective functions (LSE of groundwater storage) accompanied by the high-low factor loading scores for direction calculation, which is corrected/enlarged/scaled by the vectorized switchable step sizes, can detect the area of multiple minima, approach the global optimum, and prevent the ill-posed problem. Moreover, the average error of hydraulic conductivity identified by the Jacobian quasi-Newton between Experiment 2 and Experiment 1 rises 4.40%, while Experiment 2 increases twelve parameters over Experiment 1, indicating that the ill-posed occurrence risk rises when the total number of parameters increases. The detailed optimization process, the number of transient simulation times, and discussion of Experiment 2 are described in supplementary materials Section 3.

Study Area
The basin area of the Cho-Shui River alluvial fan which locates at the middle-west of Taiwan is 2079 km 2 . According to terrain, geology, and stratigraphy, the Cho-Shui River alluvial fan can be divided into the apex, middle, and distal fans, as shown in Error! Reference source not found. (a). The aquifer stratification in the alluvial fan is shown in Error! Reference source not found. (b). It shows that rainfall, river, and irrigation water recharge the groundwater system through the unconfined aquifer at the fan apex and that the groundwater then flows west to the coast. The groundwater system in the alluvial fan is approximately 330 m thick, which can be divided into four aquifers and three aquitards. The specific yield (Sy) of the unconfined aquifer at the apex of the fan is 0.137 to 0.237, and the storage coefficient of the confined aquifer varies from 10 −4 to 10 −3 [16]. Accordingly, the constructed conceptual groundwater flow model of the Cho-Shui River alluvial fan using MODFLOW-2005 is shown in Error! Reference source not found. (c). Reference source not found. (a). Results show that the annual pumpage for irrigation during 2012-2014 was 16.54×10 8 , 15.67×10 8 , and 15.80×10 8 m 3 , in consecutive years, while that for non-irrigation was 7.72×10 8 , 7.99×10 8 , and 7.46×10 8 m 3 , respectively. The annual total pumpage was 23.73×10 8 m 3 , of which irrigation pumpage accounted for 67.4% of the total. Figure 7 illustrates the spatial pattern of the month with maximum pumpage (sub-figures (a), (d), and (h)), the month with minimum pumpage (sub-figures (c), (f), and (i)), and the month closest to the median (sub-figures (c), (f), and (i)) during 2012-2014. Results show that the maximum pumpage in 2012 and 2013 (2.09×10 8 and 2.13×10 8 m 3 ) was in March (at the end of the dry season) and February (starting the mixed cultivation period of paddy and dry farming). Because of large water demand and insufficient surface water supply, the water sources for each demand during this month relied mainly on groundwater pumping. Moreover, the minimum pumpage during 2012-2014 (1.33×10 8 , 1.54×10 8 , and 1.18×10 8 m 3 ) were in November and December, which are the end of the wet season and the dry farming cultivation period, respectively. Moreover, in the unusually dry year in 2014, there was only one typhoon Fung-Wong, which occurred at the end of September to provide rainfall and surface water. Therefore, the maximum pumpage was in August (2.30×10 8 m 3 ), with the pumpage being closest to the mean value in April (2.08×10 8 m 3 ). Overall, high-intensity pumping is mainly concentrated in the distal alluvial fan. During the dry season and the month with high water demand, the high-intensity pumping areas extend to the interior of the middle alluvial fan.

Initially Estimated Spatiotemporal Pattern of Groundwater Recharge
After analyzing the principal components of the monthly spatiotemporal inflow and transmission during 2012-2014, the overlain multi-rank factor loading for initial groundwater recharge spatial allocation in F1, F2, F3, and F4 is from PC1 to PC5, which is determined by practical in situ conditions and the minimal value and maximal objective function gradient, as shown in Error! Reference source not found.. This figure shows that the recharge in F1 and F2 of the Cho-Shui River alluvial fan has higher loadings in the proximal fan. This is because the lithology in the apex is mostly composed of high permeability gravel, which presents a gradually decreasing trend for permeability from the apex to the tail [16]. Furthermore, the loadings' spatial pattern in F3 is larger near the northern and eastern regions. This suggests that the recharge source of F3 is mainly from recharge across the boundary, while that from aquitard leakance conductivities is less [36]. The boundary recharge in F4 is mainly distributed in the northeastern region, while the leakage from F3 is mainly concentrated in the middle and distal fans. The above analysis shows that a multi-rank loading combination can initially estimate the general trend of the spatial pattern of surface recharge and boundary recharge.

Calibration Process and Discussion
During the iterative calibration process, the high-rank components for Hessian approximation moved along the corrected direction for surface recharge identification ranging from 2 to 4, with the average of F1, F2, F3, and F4 being 2.37, 2.78, 2.53, and 2.29, respectively. For boundary recharge identification, components ranged from 5 to 9, with the average of F1, F2, F3, and F4 being 6.62, 8.43, 7.29, and 5.37, respectively; and those for hydraulic conductivity identification ranged from 8 to 16, with the average of F1, F2, F3, and F4 being 10.78, 13.61, 11.37, and 9.14, respectively. After 24 iterative calculations, the convergent process of the error descent of each aquifer is shown in Figure   1. Because of the complex groundwater recharge mechanism and hydrogeological structure in F1, the initial simulated error percentage was greater than in the other three deep aquifers. After iteration #6, the error percentage of F1 begins to converge, showing a steady decline, which is lower than the other three deep aquifers, with the average error descent rates from F1 to F4 of −0.72%, −0.46%, −0.91%, and −0.87%, respectively. Finally, after 24 iterations, the stopping criterion was met, and the error percentages from F1 to F4 are 0.11%, 23.58%, 28.99%, and 29.55%, respectively. Results showed that due to an insufficient number of groundwater head monitoring wells and low frequency of the pumping test in deep aquifers (F2, F3, and F4), the groundwater budget estimated from groundwater head fluctuation was highly uncertain. Therefore, the error percentage in F2, F3, and F4 are larger than in F1.   Error! Reference source not found. illustrates the spatial pattern of root mean square error (RMSE) between simulated and observed groundwater heads across multiple aquifers in the Cho-Shui River alluvial fan. The average groundwater head RMSE of F1, F2, F3, and F4 is 0.134 m, 0.677 m, 0.708 m, and 0.916 m, respectively. Results show that the simulated groundwater head accuracy of F1 is better than that of F2, F3, and F4, because of the high uncertainty in the hydrogeological structure of the deep aquifer. However, the correlation coefficient of F2, F3, and F4 can reach 0.659, 0.544, and 0.609, respectively, allowing the simulated hydrograph to capture the groundwater head fluctuation trend. The above-identified and simulated results illustrate that the rank three combined high-low-rank factor loading scores analyzed from the simulated storage fluctuation between adjacent iterations for the Hessian approximation not only can reduce the computation efforts but also can converge to the global minimum solution subject to the true in situ hydrogeology. The areas with larger RMSE in F1 are mainly distributed in the apex and part of the distal fan because the monitoring, estimation, and simulation of boundary recharge are difficult.
The comparisons among the simulated and observed groundwater head hydrographs across multiple aquifers are illustrated in Error! Reference source not found., in which MT station (Error! Reference source not found.a), TK station (Error! Reference source not found.b), and TH station (Error! Reference source not found.c) are located in the distal fan, middle fan, and proximal fan, respectively. Results illustrate that the RMSEs of Aquifer 2 (Error! Reference source not found.-2) in stations MT, TK, and TH are as large as 3.78, 7.40, and 3.67 times that of Aquifer 1 (Error! Reference source not found.-1), respectively; and those of Aquifer 3 (Error! Reference source not found.-3) is as large as 4.67, 3.60, and 4.87 times that of Aquifer 1 (Error! Reference source not found.-1), showing that the confined aquifer with insufficient monitoring network (e.g., middle fan of Aquifer 2 and proximal fan of Aquifer 3) would be hard to detect the actual flow condition and reduce the identified accuracy in groundwater head and parameters. From Error! Reference source not found.-4 which illustrates the quartile analysis of the simulated error across multiple aquifers, this study discovers that the groundwater head at Aquifer 1 across all sub-fans can achieve closely full fitting in groundwater head. However, the transient simulated errors at Aquifers 2 and 3 increase in which the standard deviation is 1.26m and 0.88m in MT station, respectively; that is 1.52m and 0.92m in TK station, and that is 0.61m and 1.83m in TH station. Hence, this study discovers that the exchange and flow mechanism in distal and middle fan at Aquifer 2 is more complex than that at Aquifer 3; Contrarily, the mechanisms in the proximal fan at Aquifer 3 is more complex than Aquifer 2. The current monitoring network is incapable to detect the full inflow, outflow, and groundwater flow accurately at Aquifers 2 and 3.

Conclusion and Suggestion
In a multi-layered complex groundwater model, achieving accurate spatiotemporal identification for different kinds of parameters and solving the ill-posed problem is the vital research topic in model calibration. This study proposes an SR3 VLS quasi-Newton algorithm by modifying the Levenberg-Marquardt algorithm and importing a rank three structure from the Broyden-Fletcher-Goldfarb-Shanno algorithm to solve the constrained identification of hydrogeological parameters and spatiotemporal recharge patterns simultaneously associated with an iterative approach in a large-scale multi-layered groundwater model. The recharge parameters include surface and boundary recharge, and the hydrogeological parameters consist of riverbed leakage, aquifer hydraulic, and aquitard leakance conductivities. To accelerate directional convergence, approach a global optimum, and achieve a well-posed problem, this study uses a vectorized limited switchable step size in the inverse problem of transmissive groundwater. Accordingly, to achieve the Hessian approximation, solve nonlinear ill-posed problems, and detect the global minimum area, this study imports a rank three structure, which uses a high-rank and a low-rank factor loading score analyzed from the simulated storage difference between two adjacent iterations to calculate the correction matrix for the Hessian. Moreover, it modifies the traditional algorithm from the scalar line search to vectorized multi-order derivative exact double false position bracketing using SVD-related rank and depth information for solving the vectorized step size corresponding to different hydrogeological parameters and recharges. The proposing methodology is verified by numerical experiments and is applied in a complex practical multi-layered groundwater system.
This study designed two numerical experiments to verify the capability of identifying multiple parameters and recharges through the proposing modified algorithm. Experiment 1 identified six zonal hydraulic conductivities, using a single factor loading score to approximate the Hessian. Except for six hydraulic conductivities, the parameters to be identified in Experiment 2 increased six surfaces recharges and six boundaries recharges that uses overlaid high-low factor loading scores. Results show that the iteratively objective function footprint of the three tested quasi-Newton algorithms (Jacobian, LMA, and SR3 VLS) all successfully descend at a quadratic convergence rate to search for an optimum. The identified average error percentage of the Jacobian quasi-Newton, LMA, and the SR3 VLS quasi-Newton algorithm across multiple zonal hydraulic conductivities in Experiment 1 are 2.39%, 1.78%, and 0.76%, respectively; and those across multiple surfaces recharges in Experiment 2 are 1.75%, 5.76%, and 0.09%. This demonstrates that the set objective function (LSE of groundwater storage) and the high-low combination of factor loading scores employed can detect the exact global minimum area, in which the vectorized limited switchable step sizes were striving to accelerate convergence along the corrected direction to leap from local optimum and approach a global optimum iteratively. Hence, the identified parameters by using the SR3 VLS quasi-Newton can converge to the true solution by solving nonlinear ill-posed problems. Moreover, the adopted PC Rank Depth pattern for identification illustrates that surface recharge is the largest and most direct factor affecting the fluctuation of groundwater storage with shorter travel times, and boundary recharge greatly influences the storage fluctuations with longer response lag time.
The established method was applied to the groundwater system of the Chou-Shui River alluvial fan in Taiwan. The simulated period is from January 2012 to December 2014. Results show that the RMSE during the calibration process can quickly converge in six iterations, with a decreased rate of −59.73%, −17.10%, −17.62%, and −23.96%, by considering the transmissive flow differences between single perturbation and multiple corrections of parameters simultaneously. Moreover, the average simulated storage error percentage of aquifer 1 is as small as 0.11%, with an RMSE of 0.134 m. This shows that the developed methodology not only validly detects the flow tendency and error source to quantify the spatiotemporal correction values but also precisely simulates the peak and fluctuation in groundwater head.
While calculating the multi-rank factor loading score analyzed from the simulated storage difference using SVD, if the monitoring well pattern is poor or overdense, the decomposed high-low loading scores for the rank three Hessian approximate may be far away from the real matrix. To achieve accurate and beneficial identification of parameters, optimal monitoring network design is the source of performance limitation of the developed SR3 VLS algorithm.

Supplementary Materials:
The following are available online at www.mdpi.com/2073-4441/12/4/995/s1, Figure S1: PC Rank Depth adoption diagram for identification of hydraulic conductivity during iterations across multiple parameter zones, Figure S2. Footprint of iterations vs. hydraulic conductivity in Experiment 1 and Experiment 2 optimized by (a) Jacobian quasi-Newton method; (b) LMA; (c) SR3 VLS quasi-Newton algorithm, Figure S3. PC Rank Depth adoption diagram for identification of recharge during iterations across multiple zones in Experiment 2, Figure S4.

Conflicts of Interest:
The authors declare no conflict of interest.