Pumping Well Layout Scheme Design and Sensitivity Analysis of Total Critical Pumping Rates in Coral Island Based on Numerical Model

: Groundwater on small coral islands exists in the form of freshwater lenses that serve as an important water resource for local inhabitants and ecosystems. These lenses are vulnerable to salinization due to groundwater abstraction and precipitation variation. Determination of the sustainable yield from freshwater lenses is challenging because the uncertainties of recharge and hydrogeological characteristics make it difﬁcult to predict the lens response to long-term pumping. In this study, nine pumping well layout schemes along a line are designed using the orthogonal experimental design method, and an optimal well layout scheme is determined by multi-index range analysis and comprehensive balance analysis method. The total critical pumping rates of the freshwater lens corresponding to different schemes are calculated by numerical simulation, and the sensitivity of the total critical pumping rates to hydrogeological parameters is analyzed. The results show that the calculation of the total critical pumping rates needs to be combined with the speciﬁc well layout scheme with consideration to the length of well screens, the number of wells and the distance between wells. The difference in total critical pumping rates between different schemes can be up to three times. The uncertainty of hydrogeological parameters has a great impact on the total critical pumping rates. Within the range of a 30% reduction in parameters, α and K are the key risk factors of pumping; within the range of a 30% increase in parameters, α , n e and K are the key risk factors; α - n e combined changes had the greatest impact. The management of freshwater lenses and the assessment of sustainable yield will continue to be important tasks for coral islands in the future, and this study can help with the sustainable exploitation of island freshwater lenses.


Introduction
The freshwater lens also called the Ghyben-Herzberg lens, was discovered by Ghyben [1] and Herzberg [2] in a freshwater supply study of the coastal areas of Europe in 1889 and 1901, respectively. Beneath an island, fresh groundwater may appear in the form of a "freshwater lens", which covers the seawater in a high permeability aquifer. Freshwater lenses are limited renewable groundwater resource on islands. Coral islands and atolls represent a specific subgroup of the broader category of islands. Bryan [3] listed 408 atolls spread over the oceans of the world, 294 in the Pacific Ocean, and 73 in the Indian Ocean, with the remainder distributed across the globe.
Coral islands are usually characterized by their dual-aquifer configuration. Unconformably deposited above the Pleistocene limestone reef deposit are unconsolidated or poorly consolidated Holocene sediments, primarily consisting of coral sand [4,5]. Pauw, P.S., et al. [6] used analytical solutions as a computationally faster alternative to numerical simulations to quantify saltwater coning below horizontal wells in freshwater lenses. The hydraulic conductivity of coral sand is usually between 0.023-110 m/d, and the porosity and specific yield are usually between 0.4-0.55 and 0.012-0.310, respectively, which show that coral sand has the characteristics of higher water capacity but lower specific yield compared with terrigenous sand with the same particle size range [7].
Research on freshwater lenses has been conducted using numerical simulation and physical experiments, and these results have enhanced the understanding of groundwater resources management in coral islands. Bailey et al. [8] implemented Ayers and Vachers' [4] inclusive conceptual model for atoll island aquifers in a comprehensive numerical modeling study, to evaluate the response of the freshwater lens to selected controlling climatic and geologic variables. Chui and Terry [9] simulated the damage and recovery of the freshwater lens of the atoll after storm wave scouring. Davood et al. [10] used 2D and 3D simulation methods to model the freshwater lens of Kish Island in Iran's Persian Gulf arid area and discussed the dynamics of freshwater lenses and the potential threat of climate-induced seawater intrusion. Stofberg et al. [11] employed numerical simulation to give an appraisal of the risk of freshwater lenses disappearing in the Dutch coastal regions. Vincent et al. [12] developed a script-based numerical groundwater model of a freshwater lens to evaluate the potable water supply management scenarios. Briggs et al. [13] used numerical simulation to research the driving factors of freshwater lens dynamics observed in Palmyra Atoll National Wildlife Refuge from 2008 to 2019. Numerical models have also been used to research the formation, evolution, and dynamic process of freshwater lenses in coral islands [7,[14][15][16]. Physical experiments have also been carried out to simulate the formation and exploitation of freshwater lenses [17][18][19]. A study on optimization of pumping rate and recharge to a small coral island aquifer conducted by Pallavi Banerjee and V.S Singh [20] showed that a true reflection of lens response throughout the year could only be accomplished by accounting for seasonal variability of aquifer recharge. This study also demonstrated the feasibility of using numerical models to manage fresh groundwater on islands. There is also research using real experiments on coastal aquifers. Comte J C et al. [21] focused on the hydrogeological outcomes of the research, framed within the principal socio-environmental issues identified. Bourhane A et al. [22] revealed a strong potential for both further developing coastal aquifers and initiating the prospection of inland aquifers.
While the role of real experiments on coastal aquifers has existed, there is still a study deficit on the impact of pumping on freshwater lenses, especially withoptimal well layouts and the uncertainty of hydrogeological parameters. Therefore, research is still needed to guide the management of freshwater lenses to ensure the provision of appropriate quality water.
In this study, a numerical model coupled with water flow and solute transport is utilized to simulate the up-coning caused by pumping in a coral island. An optimal well layout scheme is determined by the orthogonal experimental design method. The total critical pumping rates of the freshwater lenses are calculated by numerical simulation, and the sensitivity of the total critical pumping rates to hydrogeological parameters is analyzed. The objective of this paper is to serve as a resource for the development and utilization of island freshwater lenses, which have important practical applications.

Study Area
The coral island considered in this study is located to the southeast of Hainan Island and is the largest in the Xisha Archipelago, South China Sea, with a total area of 2.1 km 2 . The climate is tropical marine monsoon, hot and humid, with an annual average temperature of 26.5°C. The annual average precipitation of 1505 mm is sufficient; however, the seasonal distribution is irregular. In a year, there are apparent dry and wet seasons. The rainy season lasts from June to November, accounting for 80% of the annual rainfall; whereas the dry season is usually from December to May.
The topography is low and flat, with elevations of less than 6 m. It is difficult to form surface water due to the gentle slope of the terrain and the high permeability of the surface sediments. Since the middle Holocene, the reef flat, beach, sand dike, sand mat, and depression have been formed under different sedimentary dynamic conditions. The depressions in the middle of the island are 2-3 m lower than the surrounding sandbank, which are developed by lagoons. The sand mat is the main body of the island, slightly higher than the depression, with flat terrain. The sand dike is surrounded by the island and is formed by the debris sand of coral shells. The narrow beach slopes to the sea. The reef flat, which is much larger than the island, is hidden under the sea. The depth of groundwater in the coral island is about 0.3-2.9 m. Geological profiles based on 10 boreholes showed that the subsurface can be divided into a Holocene and an underlying Pleistocene unit. The contact between the Holocene and underlying Pleistocene sediments, known as the "Thurber Discontinuity" or the "Holocene-Pleistocene Unconformity" (HPU) [23], occurs at a depth about 20 m below the ground surface. The Pleistocene coral reef limestone is characterized by well-developed pores and dissolution cavities, high permeability, and easy seawater circulation. The Holocene granular sediments are mainly composed of coral clastic medium sand and coral clastic gravel sand. The hydraulic properties of the Holocene and Pleistocene units of the coral island are obviously different. The hydraulic conductivity (K) of the Pleistocene limestone is significantly higher than that of the Holocene sediments. Lens thickness for atoll islands in the Pacific and the Indian Ocean generally ranges from a few meters to 20 m, limited by the Thurber Discontinuity, as any freshwater below the discontinuity becomes salinized [24]. The maximum thickness of the freshwater lens measured on the study island is about 8-13.5 m.
A conceptual sketch of the coral island is shown as a diagram in Figure 1. The diagram includes a simplified geological setting and shows the dual nature of the aquifer system that results from the Holocene sediments overlying the Pleistocene karstified limestone. The recharge source of the groundwater is limited, and precipitation is the only natural source.
Water 2021, 13, x 3 of 16 The topography is low and flat, with elevations of less than 6 m. It is difficult to form surface water due to the gentle slope of the terrain and the high permeability of the surface sediments. Since the middle Holocene, the reef flat, beach, sand dike, sand mat, and depression have been formed under different sedimentary dynamic conditions. The depressions in the middle of the island are 2-3 m lower than the surrounding sandbank, which are developed by lagoons. The sand mat is the main body of the island, slightly higher than the depression, with flat terrain. The sand dike is surrounded by the island and is formed by the debris sand of coral shells. The narrow beach slopes to the sea. The reef flat, which is much larger than the island, is hidden under the sea. The depth of groundwater in the coral island is about 0.3-2.9 m. Geological profiles based on 10 boreholes showed that the subsurface can be divided into a Holocene and an underlying Pleistocene unit. The contact between the Holocene and underlying Pleistocene sediments, known as the "Thurber Discontinuity" or the "Holocene-Pleistocene Unconformity" (HPU) [23], occurs at a depth about 20 m below the ground surface. The Pleistocene coral reef limestone is characterized by well-developed pores and dissolution cavities, high permeability, and easy seawater circulation. The Holocene granular sediments are mainly composed of coral clastic medium sand and coral clastic gravel sand. The hydraulic properties of the Holocene and Pleistocene units of the coral island are obviously different. The hydraulic conductivity (K) of the Pleistocene limestone is significantly higher than that of the Holocene sediments. Lens thickness for atoll islands in the Pacific and the Indian Ocean generally ranges from a few meters to 20 m, limited by the Thurber Discontinuity, as any freshwater below the discontinuity becomes salinized [24]. The maximum thickness of the freshwater lens measured on the study island is about 8-13.5 m.
A conceptual sketch of the coral island is shown as a diagram in Figure 1. The diagram includes a simplified geological setting and shows the dual nature of the aquifer system that results from the Holocene sediments overlying the Pleistocene karstified limestone. The recharge source of the groundwater is limited, and precipitation is the only natural source.

Numerical Model
A two-dimensional numerical density-dependent groundwater flow and solute transport model is developed to represent the profile of a coral island using SEAWAT [25], which can effectively reduce the computational workload and improve the efficiency of calculation. The source code for SEAWAT was developed by combining MODFLOW and MT3DMS into a single program that solves the coupled flow and solute-transport equations. MODFLOW was modified to solve the variable-density flow equation by reformulating the matrix equations in terms of fluid mass rather than fluid volume, and by including the appropriate density terms. This interpretative model is used to simulate the dynamics of freshwater lenses of coral islands with similar hydrogeological conditions. The model is run for 100 years, during which time the freshwater lens develops; the model

Numerical Model
A two-dimensional numerical density-dependent groundwater flow and solute transport model is developed to represent the profile of a coral island using SEAWAT [25], which can effectively reduce the computational workload and improve the efficiency of calculation. The source code for SEAWAT was developed by combining MODFLOW and MT3DMS into a single program that solves the coupled flow and solute-transport equations. MODFLOW was modified to solve the variable-density flow equation by reformulating the matrix equations in terms of fluid mass rather than fluid volume, and by including the appropriate density terms. This interpretative model is used to simulate the dynamics of freshwater lenses of coral islands with similar hydrogeological conditions. The model is run for 100 years, during which time the freshwater lens develops; the model reaches a steady state (i.e., the morphology of lens stays relatively stable) within 50 years. The model includes 22 layers with a depth of 50 m to improve on the vertical resolution in the model simulation and is discretized into 41 columns in the horizontal direction with refinement at the boundaries. Constant head boundaries are defined along both sides of the domain to simulate the sea level, with density specified at 1025 kg/m 3 , representative of typical seawater composition. Constant concentration boundaries are assigned to the same grid cells as the constant head boundaries with chloride concentrations of 19000 mg/L. Assuming that the loose sediments of the coral island are initially saturated by seawater, the initial chloride concentration of the entire model domain is specified at 19,000 mg/L. The ground surface is provided by the topographic map, and the lithology is based on boreholes. The model grid is shown in Figure 2.
Water 2021, 13, x 4 of 16 reaches a steady state (i.e., the morphology of lens stays relatively stable) within 50 years. The model includes 22 layers with a depth of 50 m to improve on the vertical resolution in the model simulation and is discretized into 41 columns in the horizontal direction with refinement at the boundaries. Constant head boundaries are defined along both sides of the domain to simulate the sea level, with density specified at 1025 kg/m 3 , representative of typical seawater composition. Constant concentration boundaries are assigned to the same grid cells as the constant head boundaries with chloride concentrations of 19000 mg/L. Assuming that the loose sediments of the coral island are initially saturated by seawater, the initial chloride concentration of the entire model domain is specified at 19,000 mg/L. The ground surface is provided by the topographic map, and the lithology is based on boreholes. The model grid is shown in Figure 2. Hydraulic conductivity of the principal aquifer is determined based on the lithology and previous studies [7,15]. Holocene sediments were divided into 12 model layers: 1-3 layers represent coral clastic medium sand, with the horizontal hydraulic conductivity (K) of 60 m/d; 4-7 layers represent coral clastic gravelly sand, with the horizontal hydraulic conductivity (K) of 110 m/d; 8-12 layers represent coral breccia, with the horizontal hydraulic conductivity (K) of 150 m/d. Pleistocene limestone, which has good permeability because of many pores and dissolution cavities of different sizes, includes 10 model layers (13)(14)(15)(16)(17)(18)(19)(20)(21)(22), and the horizontal hydraulic conductivity (K) is 1000 m/d [26]. The effective porosity (ne) is 0.25-0.45 and the specific yield (Sy) is 0.1-0.2.
Recharge is applied to the top layer of the model with a chloride concentration of 0 g/L to simulate the monthly average recharge to the aquifer. Recharge is the only input of freshwater to the hydrogeological system and, therefore, is the main mechanism by which the simulated freshwater lens develops in the model [27]. The monthly net recharge is calculated using the precipitation recharge coefficient (α). Before the model reaches a steady state, the pumping wells are not activated. After the model reaches a steady state, the wells are activated for continuous pumping for 20 years to simulate the saltwater upconing. The hydrogeological parameters assigned to the SEAWAT model are summarized in Table 1. Hydraulic conductivity of the principal aquifer is determined based on the lithology and previous studies [7,15]. Holocene sediments were divided into 12 model layers: 1-3 layers represent coral clastic medium sand, with the horizontal hydraulic conductivity (K) of 60 m/d; 4-7 layers represent coral clastic gravelly sand, with the horizontal hydraulic conductivity (K) of 110 m/d; 8-12 layers represent coral breccia, with the horizontal hydraulic conductivity (K) of 150 m/d. Pleistocene limestone, which has good permeability because of many pores and dissolution cavities of different sizes, includes 10 model layers (13)(14)(15)(16)(17)(18)(19)(20)(21)(22), and the horizontal hydraulic conductivity (K) is 1000 m/d [26]. The effective porosity (n e ) is 0.25-0.45 and the specific yield (S y ) is 0.1-0.2.
Recharge is applied to the top layer of the model with a chloride concentration of 0 g/L to simulate the monthly average recharge to the aquifer. Recharge is the only input of freshwater to the hydrogeological system and, therefore, is the main mechanism by which the simulated freshwater lens develops in the model [27]. The monthly net recharge is calculated using the precipitation recharge coefficient (α). Before the model reaches a steady state, the pumping wells are not activated. After the model reaches a steady state, the wells are activated for continuous pumping for 20 years to simulate the saltwater up-coning. The hydrogeological parameters assigned to the SEAWAT model are summarized in Table 1.

Orthogonal Experimental Design
The application of orthogonal experimental design is efficient, fast and economic, and has been widely used in many fields [28][29][30][31]. Orthogonal experimental design is an experimental design method to study multiple factors and multiple levels. According to orthogonality, some representative experimental scenarios are selected from the comprehensive experiments. When three or more components are involved in the experiment, and there may be interactions between them, the workload increases and gets difficult to implement. Orthogonal experimental design is undoubtedly a better choice to solve this problem and can achieve equivalent results with a minimum number of tests.  3 3 )). In the three factors and three levels orthogonal table format, such as in Table 3, the first column represents the number of the experiments, and the specific values in the table represent the combination of different levels of factors. Orthogonality means that each level of one factor can meet once with each level of another factor. The orthogonal table has two characteristics: • Different numbers in each column appear at the same time. For example, in the three factors and three levels orthogonal table, any column has 1, 2, and 3, and the number of occurrences in each column is equal.

Multi-index Range Analysis
Range analysis is an analysis method used for orthogonal experiment results, which is carried out after the data collection. Range (R) is the difference between the maximum and minimum values, which is used to represent the measures of variation in statistical data. Range analysis determines the primary and secondary order of different factors affecting the experimental results by comparing R. The larger the R, the larger the swing range of the factor in different levels of data, indicating that the factor has a more significant impact on the results. Range analysis results can be recorded in the range analysis table (Table 4). If there are multiple indicators, range analysis is required for each indicator. The operation process is as follows. Firstly, the experiment results of each factor at the same level of a single indicator are summed, respectively, and the sum obtained is denoted as K. Then, the mean values of K1, K2, and K3 are calculated, respectively, which are denoted as _ K1, _ K2, _ K3, and the optimal level of each factor is determined by _ K. R is then calculated and the importance of factors are sorted accordingly. Finally, the above operation is repeated for each indicator. For example (Table 4), for factor A, when the level is 1, _ K is the largest, which means level 1 is the optimal level of factor A for indicator 1. Since RB > RA > RC, the influence of factors on indicator1 is ranked as BAC.

Comprehensive Balance Analysis
The comprehensive balance analysis method is to investigate the influence of each factor on each index, after which analysis and comparison can be performed to determine the best level, to obtain the best test scheme. The operation process is as follows. Firstly, the multi-index range analysis is carried out to obtain a better combination of several factors. In this process, the primary and secondary factors of each indicator are identified, and the level of each factor is determined. Then, the overall optimal factor level combination is determined based on weighing the impact on other indicators when the optimal level of a factor is selected.

Pumping Wells Layout Scheme Design of the Freshwater Lens
"Sustainable exploitation" is an extension of "allowable exploitation", based on the concept of sustainable development (UNCED, 1992) [32]. In 1999, Alley et al. [33] defined sustainable groundwater development as "long-term and permanent development and utilization of groundwater without serious social, economic and environmental consequences". The concept of sustainable exploitation highlights environmental factors and emphasizes the renewable and sustainable utilization of groundwater resources. This study proposes that the essence of sustainable exploitation for the groundwater system of coral islands is that after a change in source and sink terms due to human intervention, the freshwater lens system changes from the original dynamic equilibrium to a new dynamic equilibrium, and the latter is acceptable. Because the variation in the groundwater table is at centimeter-level and the variation in the concentration isoline is at meter level in many scenarios of this simulation, this study regards the saltwater up-coning height as the main limiting factor of sustainable exploitation (Figure 3). A pumping well extracting from the freshwater lens can cause the saltwater to move upwards towards the well, a phenomenon known as saltwater up-coning [6,34]. A schematic diagram of the saltwater up-coning process is shown in Figure 4. In this study, when the groundwater numerical model reaches the steady-state, the corresponding pumping rate when the rising height just reaches the bottom of the well screen is defined as the critical pumping rate (QC). The sum of QC of all pumping wells can then be expressed by the total critical pumping rate (QT), and the sustainable yield cannot exceed QT. Tm is the maximum thickness of the freshwater lens under pumping conditions, and TC is the central thickness (Figure 4). The coordinate corresponding to the Tm changes with pumping conditions, but the position corresponding to the TC will not change. Due to pumping, the TC doesn't coincide with the Tm. In this study, the length of the well screens, the number of wells, and the distance between wells are selected as the factors of the design scheme of pumping wells, and each factor has three levels ( Table 5). One of the typical well layout schemes is shown in Figure 4.  A pumping well extracting from the freshwater lens can cause the saltwater to move upwards towards the well, a phenomenon known as saltwater up-coning [6,34]. A schematic diagram of the saltwater up-coning process is shown in Figure 4. In this study, when the groundwater numerical model reaches the steady-state, the corresponding pumping rate when the rising height just reaches the bottom of the well screen is defined as the critical pumping rate (Q C ). The sum of Q C of all pumping wells can then be expressed by the total critical pumping rate (Q T ), and the sustainable yield cannot exceed Q T . T m is the maximum thickness of the freshwater lens under pumping conditions, and T C is the central thickness ( Figure 4). The coordinate corresponding to the T m changes with pumping conditions, but the position corresponding to the T C will not change. Due to pumping, the T C doesn't coincide with the T m . In this study, the length of the well screens, the number of wells, and the distance between wells are selected as the factors of the design scheme of pumping wells, and each factor has three levels ( Table 5). One of the typical well layout schemes is shown in Figure 4.
In this study, we use a chloride concentration of 250 mg/L as the upper threshold for fresh groundwater, which is the World Health Organization (WHO) [35] recommendation for drinking water in 2011. However, other researchers have also used 500 mg/L and 600 mg/L as the upper threshold of chloride concentration [14,36,37]. When the model runs stably without pumping, the freshwater flow is stressed by saltwater when the boundary is moved on both sides, whereas saltwater moves upward driven by buoyant effects in the middle of the model [38]. Three variables are selected as evaluating indicators in the numerical test: total critical pumping rate (Q T ), the central thickness of the freshwater lens (T C ), and maximum thickness of the freshwater lens (T m ) (Figure 4). The orthogonal experimental design scheme and its corresponding numerical test results are shown in Figure 5 and Table 6.
TC doesn't coincide with the Tm. In this study, the length of the well screens, the number of wells, and the distance between wells are selected as the factors of the design scheme of pumping wells, and each factor has three levels ( Table 5). One of the typical well layout schemes is shown in Figure 4.    In this study, we use a chloride concentration of 250 mg/L as the upper threshold for fresh groundwater, which is the World Health Organization (WHO) [35] recommendation for drinking water in 2011. However, other researchers have also used 500 mg/L and 600 mg/L as the upper threshold of chloride concentration [14,36,37]. When the model runs stably without pumping, the freshwater flow is stressed by saltwater when the boundary is moved on both sides, whereas saltwater moves upward driven by buoyant effects in the middle of the model [38]. Three variables are selected as evaluating indicators in the numerical test: total critical pumping rate (QT), the central thickness of the freshwater lens (TC), and maximum thickness of the freshwater lens (Tm) (Figure 4). The orthogonal experimental design scheme and its corresponding numerical test results are shown in Figure  5 and Table 6.  The numerical experiment results show that different pumping well layout schemes have a great impact on QT. The QT of test 3 is 4.2 times that of test 8 (Table 6). Of course, QT is not the only factor to judge whether the scheme is the best or not. It is also necessary to comprehensively analyze the thickness of the freshwater lens under the pumping scheme. If the thickness of the freshwater lens becomes too shallow, the current pumping  The numerical experiment results show that different pumping well layout schemes have a great impact on Q T . The Q T of test 3 is 4.2 times that of test 8 (Table 6). Of course, Q T is not the only factor to judge whether the scheme is the best or not. It is also necessary to comprehensively analyze the thickness of the freshwater lens under the pumping scheme. If the thickness of the freshwater lens becomes too shallow, the current pumping scheme is not reasonable, although a large pumping rate may be obtained. In principle, the sustainable yield shall not exceed Q T . The determination of Q T can provide a reference for the decision-making of sustainable yield.
According to the calculation results and difference analysis of orthogonal tests (Table 7), RB > RA > RC for Q T and T m , RA > RB > RC for T C . The results indicate that the number of wells, and the length of well screens, are the most critical factors in the design of the well layout scheme, followed by the distance between wells. Although the more wells, the greater the Q T , if the number of wells is increased without considering the distance between wells, the T C will be greatly reduced and the risk of damage to the freshwater lens will increase. Therefore, the number of wells needs to be limited. The longer the well screen and the deeper the buried depth at the bottom of the well screen, the smaller the Q T , so the design of the length of the well screen needs to be cautious. If long-term pumping is implemented under the conditions of short well screens with Q T , it is likely to reduce the T C to less than 1/3 of the original thickness, which increases the risk of damage to the freshwater lens. Therefore, it is necessary to control the pumping intensity of shallow wells.
In order to reduce T C and T m as little as possible and obtain a larger Q T , the comprehensive balance analysis method is used to select the optimization scheme. Compared with Q T and T m , factor A has a more significant influence on T C . It can be seen from Table 4. that A2 is the best choice. The influence of A on Q T and T m occupies second place among the three factors, with A1 and A3 being the best choices. However, when A1 is selected, T C and T m are very thin, so A1 is not considered; when A2 is selected, T m is 2.6% less than that of A3; when A3 is selected, not only is T C reduced by 4%, but Q T is also reduced by 6% compared with A2, so A2 is considered comprehensively. The influence of factor B on Q T ranks first among A, B, and C. However, when B3 is taken, T C is less than 1/3 (2.5 m) of the original thickness. When B2 is taken in turn, there is the same problem, so B1 is taken finally, which is also the best level for T C . According to Table 4, the influence of factor C on Q T , T C , and T m takes third place, and C1 and C3 can be taken into consideration. When C1 is selected, T C and T m increase by 25.5% and 9.4%, respectively; When C3 is selected, Q T will increase by 39.5% compared with C1, but T C and T m are only reduced by 20.3% and 8.6%, respectively, so C3 is a better choice. As a result, the best scheme is A2B1C3, that is, the linear well layout scheme with a 3 m well screen length, two wells, and 200 m of the distance between wells (Test 9). This distance accounts for about 10% of the island width. Under this scheme, Q T is 0.33 m 3 /d. T C is 4.6 m, which accounts for 60% of the original thickness, and T m is 9.8 m, which accounts for 86% of the original thickness.

Uncertainty Analysis of Q T and Main Influence Factors
Due to the uncertainty of hydrogeological parameters, Q T in the above-mentioned pumping wells layout scheme is also uncertain. Hydraulic conductivity (K), vertical dispersion coefficient (DI), precipitation infiltration coefficient (α), effective porosity (n e ), and specific yield (S y ) are the influencing factors: Sensitivity analysis is used to analyze the uncertainty of Q T and identify the main risk factors. Sensitivity can measure the impact of the change in influencing factors on the indicators. The global sensitivity analysis can fully account for the effect of parameter interaction on the simulation results. In this study, the method proposed by Morris [39] in 1991 was adopted. By adjusting only one parameter at a time, the influence of each parameter on the overall results is calculated in turn, and the importance of each parameter is effectively prioritized. The influence of parameter changes on the overall results is measured by sensitivity (S i ), expressed as: where S i is the sensitivity of function y(x) to variable X i ; y 0 is the value of Q T before the parameter changes; X 0 is the initial parameter. Among them, the larger the value of |S i | is, the more sensitive the variable X i is, and the greater the influence on the function is.
When K, DI, α, n e , and S y change from 0 to 30% respectively, the |S i | of α is the largest, ranging from 0 to 2.020, therefore, changes in α will have a greater impact on Q T compared with other four parameters with the same amplitude. |S i | of n e is in the second place, ranging from 0 to 0.566, and the influence of n e is much smaller than that of α. The other three parameters have less influence, so the focus should be on the changes in α and n e . When K, DI, α, n e, and S y change from -30% to 0, the |S i | of α and K are significantly greater than the other three parameters, ranging from 0 to 1.202 and 0 to 0.949, respectively. Therefore, the focus should be on the changes of α and K on the determination of Q T . In conclusion, from −30% to 30%, α, K, and n e are the main influencing factors in the determination of Q T .
The results of the global sensitivity analysis of α-K, α-n e , K-n e and α-K-n e show that the comprehensive influence of α-n e is the largest, and the corresponding |S i | ranges from 0 to 2.182. When the parameter amplitude is −30% to 0, the comprehensive influence of α and n e on Q T is less than the sum of their respective influence, but when the parameter amplitude is 0 to 30%, the comprehensive influence of α-n e on Q T is greater than the sum of their respective influence. The analysis results of local sensitivities of all parameters and joint sensitivities of main influencing parameters are shown in Table 8 and Figure 6. ranging from 0 to 0.566, and the influence of ne is much smaller than that of α. The other three parameters have less influence, so the focus should be on the changes in α and ne. When K, DI, α, ne, and Sy change from -30% to 0, the |Si| of α and K are significantly greater than the other three parameters, ranging from 0 to 1.202 and 0 to 0.949, respectively. Therefore, the focus should be on the changes of α and K on the determination of QT. In conclusion, from -30% to 30%, α, K, and ne are the main influencing factors in the determination of QT.
The results of the global sensitivity analysis of α-K, α-ne, K-ne and α-K-ne show that the comprehensive influence of α-ne is the largest, and the corresponding |Si| ranges from 0 to 2.182. When the parameter amplitude is -30% to 0, the comprehensive influence of α and ne on QT is less than the sum of their respective influence, but when the parameter amplitude is 0 to 30%, the comprehensive influence of α-ne on QT is greater than the sum of their respective influence. The analysis results of local sensitivities of all parameters and joint sensitivities of main influencing parameters are shown in Table 8 and Figure 6.  The increase in α and n has positive impacts on QT, but the increase in K has negative impacts ( Figure 7). The reason for this is that α affects the amount of the net recharge. A larger α means that the freshwater lens can obtain a larger hydraulic gradient to maintain The increase in α and n has positive impacts on Q T , but the increase in K has negative impacts ( Figure 7). The reason for this is that α affects the amount of the net recharge. A larger α means that the freshwater lens can obtain a larger hydraulic gradient to maintain or even increase the thickness of freshwater, and has a larger maximum permissible rising height of saltwater, resulting in a larger Q T . With the increase in K, the thickness of the freshwater lens will become thinner, and the maximum allowable rising height of saltwater will decrease, which requires pumping rates to be relatively reduced to protect the wells from salinization.
Water 2021, 13, x 13 of 16 or even increase the thickness of freshwater, and has a larger maximum permissible rising height of saltwater, resulting in a larger QT. With the increase in K, the thickness of the freshwater lens will become thinner, and the maximum allowable rising height of saltwater will decrease, which requires pumping rates to be relatively reduced to protect the wells from salinization. Under the most unfavorable parameter combination (α and n decrease by 30%, K increases by 30%), QT is 0.174 m 3 /d, which is only 52.7% of the original (0.33 m 3 /d), so the original calculated QT is not absolutelyreliable; under the most favorable parameter combination (α and n increases by 30%, K decreases by 30%), QT is 0.573 m 3 /d, which is 73.63% higher than the original QT. Therefore, the range of QT is from 0.174 to 0.573 m 3 /d. The identification of the main influencing factors can help identify the key risk sources that affect the determination of QT, from the standpoint of hydrogeological parameters. The results show that the key risk sources of QT are the uncertainties of α, n, and K. The influence of parameters on the determination of QT is related to the variation range of the parameters, within the range of a 30% reduction in parameters, α is the most important risk factor of pumping, and K is the second; within the range of a 30% increase in parameters, α is also the most important risk factor, followed by n and K. The uncertainty of the hydrogeological parameters has a great impact on QT, but their spatial variabilities haven't been taken into consideration, which may also influence the results, and is meant for relevant analysis.
The 2D model assumes that recharge is stored in the profile and there is no flow movement perpendicular to the profile. This generalization may lead to a slight overestimation of reserves. However, the overestimation is generally within a reasonable range. Numerical models can be used as a powerful tool to predict and evaluate the water resources of coral islands, but it is difficult to avoid errors in comparison to reality, especially in the conceptualization of hydrogeological conditions of coastal aquifers. The determination of sustainable yield should be combined with the specific well layout schemes, and the uncertainty of hydrogeological parameters should also be taken into account, because of the complex hydrogeological conditions of coral islands.  Under the most unfavorable parameter combination (α and n decrease by 30%, K increases by 30%), Q T is 0.174 m 3 /d, which is only 52.7% of the original (0.33 m 3 /d), so the original calculated Q T is not absolutelyreliable; under the most favorable parameter combination (α and n increases by 30%, K decreases by 30%), Q T is 0.573 m 3 /d, which is 73.63% higher than the original Q T . Therefore, the range of Q T is from 0.174 to 0.573 m 3 /d. The identification of the main influencing factors can help identify the key risk sources that affect the determination of Q T , from the standpoint of hydrogeological parameters. The results show that the key risk sources of Q T are the uncertainties of α, n, and K. The influence of parameters on the determination of Q T is related to the variation range of the parameters, within the range of a 30% reduction in parameters, α is the most important risk factor of pumping, and K is the second; within the range of a 30% increase in parameters, α is also the most important risk factor, followed by n and K. The uncertainty of the hydrogeological parameters has a great impact on Q T , but their spatial variabilities haven't been taken into consideration, which may also influence the results, and is meant for relevant analysis.
The 2D model assumes that recharge is stored in the profile and there is no flow movement perpendicular to the profile. This generalization may lead to a slight overestimation of reserves. However, the overestimation is generally within a reasonable range. Numerical models can be used as a powerful tool to predict and evaluate the water resources of coral islands, but it is difficult to avoid errors in comparison to reality, especially in the conceptualization of hydrogeological conditions of coastal aquifers. The determination of sustainable yield should be combined with the specific well layout schemes, and the uncertainty of hydrogeological parameters should also be taken into account, because of the complex hydrogeological conditions of coral islands.

Summary and Concluding Remarks
Groundwater resources of coral islands are critical for water supply and have great economic and social value. This study provides an optimal well layout scheme by using the orthogonal experimental design and numerical simulation and analyzes the sensitivity of the pumping rate to the uncertainty of hydrogeological parameters. On this basis, the main influencing factors affecting the determination of pumping rates are selected, and the possible risks associated with the uncertainty of hydrogeological parameters are suggested. The results show that: (1) The determination of Q T should be based on the specific well layout, with consideration of the length of well screens, the number of wells, and the distance between wells, because the calculation results of Q T corresponding to different well layouts can differ by three times. (2) The orthogonal experimental design is adopted to analyze the important elements of well layout, and thereby the results based on multi-index extremum range analysis and comprehensive balance analysis show that the number of wells and the length of well screens are the most critical factors in the design of well layout scheme, followed by the distance between wells. (3) In this study, the optimal well layout is selected from nine schemes. The length of the well screen of the scheme is 3 m, the number of wells is two and the distance between wells is 200 m, which accounts for about 10% of the island width. For the optimal scheme, Q T is 0.33 m 3 /d, T C is 4.6 m, which accounts for 60% of the original thickness, and T m is 9.8 m, which accounts for 86% of the original thickness. (4) The sensitivities of the hydrogeological parameters with respect to Q T are significant.
Within the range of a 30% reduction in parameters, α and K are the key risk factors of pumping; within the range of a 30% increase in parameters, α, n and K are the key risk factors; α-n combined changes had the greatest impact. Under conditions of the most unfavorable parameter combination, Q T is only 52.7% of the original value.
The management of freshwater lenses and the assessment of sustainable yields, together with the uncertainty of hydrogeological factors, will continue to be the primary activity of coral islands. This paper only focuses on the variation and combination of parameters, without considering the spatial variability of hydrogeological parameters in the basic model. Spatial variability of the hydrogeological parameters may also have an important influence on the optimization of well layout and the sensitivity of pumping rate and is recommended for future work.