Comparison of a simple hydrostatic and a data-intensive 3D numerical modeling method of simulating sea-level rise induced groundwater inundation for Honolulu, Hawai’i, USA

Groundwater inundation (GWI) is a particularly challenging consequence of sea-level rise (SLR), as it progressively inundates infrastructure located above and below the ground surface. Paths of flooding by GWI differ from other types of SLR flooding (i.e., wave overwash, storm-drain backflow) such that it is more difficult to mitigate, and thus requires a separate set of highly innovative adaptation strategies to manage. To spur consideration of GWI in planning, data-intensive numerical modeling methods have been developed that produce locally specific visualizations of GWI, though the accessibility of such methods is limited by extensive data requirements. Conversely, the hydrostatic (or ‘bathtub’) modeling approach is widely used in adaptation planning owing to easily accessed visualizations (i.e., NOAA SLR Viewer), yet its capacity to simulate GWI has never been tested. Given the separate actions necessary to mitigate GWI relative to marine overwash, this is a significant gap. Here we compare a simple hydrostatic modeling method with a more deterministic, dynamic and robust 3D numerical modeling approach to explore the effectiveness of the hydrostatic method in simulating equilibrium aquifer effects of multi-decadal sea-level rise, and in turn GWI for Honolulu, Hawai’i. We find hydrostatic modeling in the Honolulu area and likely other settings may yield similar results to numerical modeling when referencing the local mean higher-high water tide datum (generally typical of flood studies). These findings have the potential to spur preliminary understanding of GWI impacts in municipalities that lack the required data to conduct rigorous groundwater-modeling investigations. We note that the methods explored here for Honolulu do not simulate dynamic coastal processes (i.e., coastal erosion, sediment accretion or changes in land cover) and thus are most appropriately applied to regions that host heavily armored shorelines behind which GWI can develop.


Introduction
Sea-level rise (SLR) presents inevitable challenges for low-lying coastal municipalities (Hallegatte et al 2013, Hinkel et al 2014. Even as SLR projections evolve, many researchers have concluded that flooding will grow progressively damaging within decades Strauss 2017, Sweet et al 2018). This is especially true for regions where rates of SLR exceed the global mean (i.e., East and Gulf Coasts of the US) (Sweet et al 2017). Hightide flooding is already problematic at these locations resulting in drainage failure, road closure, and the deterioration of municipal infrastructure (Sweet and Park 2014).
uppermost 100 to 200 m of the aquifer is unconfined (Macdonald et al 1983, and influenced directly by rainfall and near-shore sea-level fluctuations produced by tides, wave set-up, and longer period sealevel variations (Ponte 1994, Wu et al 1996, Yin et al 2001, Gonneea et al 2013, Habel et al 2017.
Groundwater in the unconfined aquifer is not potable owing to elevated salinity and urban contamination, and therefore is only extracted for use in small-scale irrigation and cooling towers (Whittier et al 2010). Flow in the unconfined aquifer is driven by the pressure gradient from the underlying confined aquifer into the unconfined aquifer, and through surficial recharge by rainfall, leakage of water-conveyance infrastructure, and small-scale irrigation (Engott et al 2017). The subsurface geology comprises post-erosional volcanics, alluvial debris, artificial fill, and reefal carbonates related to Quaternary sea-level high stands (Stearns and Vaksvik 1935, Ferrall 1976, Munro 1981, Finstick 1996, Izuka et al 2018.

Sea-level rise projections
The IPCC Fifth Assessment Report (Church et al 2013) estimates that global mean sea level could reach magnitudes ranging from 0.52 to 0.98 m by 2100 (relative to 1986-2005) under Representative Concentration Pathway 8.5, the 'business as usual scenario.' However, studies incorporating ice-shelf hydrofracturing and icecliff collapse mechanisms, triggered under high-emissions scenarios, indicate the potential for higher SLR late this century (i.e., Kopp et al 2017). Additionally, owing to global gravitational effects (i.e., ice fingerprinting), SLR particular to the Hawaiian Islands will likely exceed the global mean (Spada et al 2013, Kopp et al 2014. For the purpose of this comparison study, when simulating GWI, we consider a SLR magnitude of 1 m, consistent with the intermediate scenario of Sweet et al (2017).

Methods
Model accuracy is assessed by comparing simulated present-day groundwater levels produced using the hydrostatic and 3D numerical methods to groundwater-level observations compiled within the study area from monitoring wells. We consider Honolulu's mean sea-level (MSL) tide stage and MHHW tide stage (0.33 m+MSL) in model construction (National Oceanic and Atmospheric Administration NOAAa 2017).
The approach used to characterize areas vulnerable to GWI was adopted from similar studies (Rotzoll and Fletcher 2013, Cooper et al 2015, Habel et al 2017 such that GWI is characterized in locations where groundwater elevations exceed the ground-surface elevations. Ground-surface elevations were simulated using a digital elevation model (DEM) that was constructed by merging rasterized 2013 NOAA DEM tiles. The tiles were produced by NOAA using LiDAR ground-return elevation data referenced to local mean sea level (LMSL) (National Oceanic and Atmospheric Administration NOAAb 2017). The elevation data has 0.15 m and 1 m vertical and horizontal resolution, respectively (Office for Coastal Management 2019).

Groundwater elevation data and subsets
Observations of groundwater elevations were compiled for use in 3D model calibration and to test the accuracy of model simulations. Groundwater-level data available for this study included 247 discrete water-level observations obtained from Hawai'i Department of Health Leaky Underground Storage Tank records (State of Hawaii Department of Health DOH 2018), and 73 sets of continuous water-level observations compiled from local hydrologic studies. For model calibration, 193 discrete water-level measurements were spatially subsampled and 49 continuous records were temporally subsampled from the original dataset, and the remaining measurements were used for cross-validation analysis.
Cross-validation compares simulated groundwater levels with two observation subsets representative of the MSL and MHHW scenarios. Compared to oscillations in ocean water levels, groundwater-level oscillations are attenuated, decreasing in amplitude and increasing in temporal lag as oscillations propagate inland. The influence of ocean oscillations was quantified for each set of continuous observations using the methods of Habel et al (2017) in which temporal lag was evaluated by cross correlating tidal signals at the Honolulu tide station with tidal signals observed in the groundwater data, and tidal efficiency was calculated using linear leastsquares regression of lag-corrected groundwater time series to tidal-signal data. To accommodate the tidalresponse phenomenon, data in each subset were chosen with consideration of the tidal elevation such that the tide was within 4 cm of the respective MSL or MHHW tidal scenario at the lag-corrected time of data collection. Subsamples consist of 43 discrete measurements representing the MSL scenario, 11 discrete measurements representing the MHHW scenario, and 24 continuous observations used to represent both scenarios. All discrete observations were corrected for anomalous sea-surface height using the methods of Habel et al (2017). The 24 continuous measurements were processed to represent MSL by calculating average recorded water levels during times with negligible rainfall and non-anomalous sea-surface levels. Measurements were processed to represent MHHW by calculating the average reconstructed tidal influence produced during the MHHW tide stage. Reconstructions of tidal influence were produced using UTIDE (Codiga 2011).

Hydrostatic model construction
In hydrostatic modeling, the groundwater hydraulic gradient and attenuated tidal response with distance inland are ignored (i.e., Marcy et al 2011, Strauss et al 2012. Thus, we set the water table to an elevation equal to simulated sea level, referenced to the MSL datum. For example, for cases in which sea level is simulated at an elevation of 0 m, the water table inland is assumed to be located at an elevation of 0 m across the entire study area. For a simulated SLR of 1 m, the water table everywhere is assumed to be at an elevation of 1 m.

3D Numerical model construction
Groundwater simulation using the 3D numerical method with MODFLOW 2005 (Harbaugh 2005) follows Habel et al (2017), but has been extended to represent a larger study area encompassing the Primary Urban Center of Honolulu. The numerical model simulates steady-state conditions of the water table at current mean sea level and a 1-m increase in sea level. Simulated subsurface hydrogeologic conditions are based on conditions determined in regional studies (Ferrall 1976, Munro 1981, Finstick 1996, Oki 2005, Rotzoll and El-Kadi 2007. The model was calibrated using discrete and continuous water-level measurements discussed previously (figure 1).
The extended model consists of 48,483 active, 100-m uniform grid cells and three layers, representing unconsolidated caprock (model layer 1), consolidated caprock (model layer 2), and basalt (model layer 3) hydrogeological units, respectively. The inland boundary is defined by the 0 m elevation contour (figure 2) that represents the uppermost extent of the basalt aquifer (Rotzoll and El-Kadi 2007). The seaward boundary is defined by the 200-m depth contour of 2013 US Army Corps of Engineers LiDAR bathymetry data (National Oceanic and Atmospheric Administration NOAAb 2017). The top of the unconsolidated caprock unit (model layer 1) is defined by mosaicked 2013 NOAA LiDAR topography data and 2013 US Army Corps of Engineers LiDAR bathymetry data (NOAAb 2017), with a specified thickness of 10 m based on the approximate depth in which consolidated caprock material has been encountered (Ferrall 1976, Munro 1981, Finstick 1996. The consolidated caprock unit (model layer 2) extends from the base of the unconsolidated caprock unit to the uppermost extent of the basalt unit (model layer 3) as defined by elevation data that represents the uppermost extent of the basalt aquifer (Rotzoll and El-Kadi 2007). Details of the flow in the basalt aquifer unit are beyond the scope of this study. Therefore, the basalt aquifer is represented by a thin unit that extends an arbitrary 1 m below the base of the consolidated caprock unit and was included exclusively to simulate flow from the basalt into the caprock aquifer in which model layers 1 and 2 represent the caprock aquifer.
The model domain is bounded on the bottom and the sides by no-flow boundaries (with the exception of the inland boundary for the bottom unit); the upper-boundary is a specified recharge boundary; the inland lateral boundary of the bottom unit (model layer 3) is a specified head boundary to simulate flow from the basalt unit to the upper caprock units (figure 2). Specified-head values were based on simulations of confined groundwater flow in southern O'ahu (Rotzoll and El-Kadi 2007). Seaward of the 0-m land-surface elevation contour, current sea-level conditions were simulated using a specified general-head boundary at the ocean bottom with a conductance of 10 m 2 d −1 . A 1-m increase in sea level was simulated by re-evaluating the general-head boundary seaward of the 1-m elevation contour, and by increasing the elevation of hydraulic head from 0 m to 1 m; conductance was not changed.
Well locations (figure 2) and withdrawal rates available from the State Commission on Water Resource Management were adopted from existing groundwater-flow models representative of the Honolulu aquifer (Rotzoll and El-Kadi 2007) and only those wells pumping from the caprock were considered following the application of Habel et al (2017). Well withdrawal rates were defined as the arithmetic mean of respective pumping rates from 1996 to 2005. Recharge data were acquired from the mean annual water-budget model for the Island of O'ahu, Hawai'i (Engott et al 2017), which simulated hydrological processes including rainfall, fog interception, irrigation, direct runoff, return flows from septic systems, and evapotranspiration.
Hydraulic-conductivity values for model layer 1 representing the unconsolidated caprock unit were estimated using the nonlinear inverse modeling utility, PEST in which Tikhonov preferred homogeneous regularization was used (Doherty and Hunt 2010). As part of this approach, pilot points were established on 500-m grid across the study area, totaling 361 points. All post-calibration values applied to the unconsolidated caprock unit were within the range of values previously observed for the study area, ranging from 0.001 to 854 m d −1 (Finstick 1996)   Following model calibration, the simulated mean residual water level and root-mean-squared error were 0.04 m and 0.12 m, respectively.
The specific limitations of the modeling methodology are summarized in Habel et al (2017) in which the main limitations include: • The model is steady-state and thus does not assess time-dependent hydrological processes such as variations in recharge, pumping rates, boundary flows and groundwater storage, and aperiodic short-term changes in sea level by phenomena such as storm-surges, tsunamis, etc; • MODFLOW-2005 assumes a uniform density of water, and thus does not incorporate the influence of density-driven fluid flow such as mixed seawater and freshwater flows; • The model does not consider flow that occurs in the unsaturated zone or surface-water flow, evaporation from surface-water sources, and ponding or routing of waters that occurs once groundwater has breached the ground surface; • The model does not consider dynamic changes in landscape (i.e., erosion) produced by SLR.

Tidal application
For the hydrostatic method, tidal influence was assessed by elevating the groundwater level to the tide stage elevation being simulated, thus natural attenuation of the tidal signal in an inland direction was not considered and water levels were not adjusted when simulating the MSL tide stage. As such, this method would tend to overestimate the inundation derived from groundwater when considering tides as it neglects attenuated head responses to tidal forcing. For the 3D numerical method, tidal influence was evaluated by performing regression analyses to compute analytical solutions that assign tidal efficiency as an exponential function of distance from the coastline following the assessment in Habel et al (2017). Tidal efficiencies representing six subzones were calculated by comparing tidal amplitudes recorded at the NOAA Honolulu Tide station (National Oceanic and Atmospheric Administration NOAAa 2017) to those observed at continuous monitoring wells. The six subzones represent the eastern, middle, and western extent of the study area, in which the western extent was further divided into four subzones. Subdivisions were made in the western extent to represent unique patterns of observed tidal efficiencies that correlated with distinct geologic strata that constitute the shallow geology (i.e., fill, beach deposits, lagoon and reef deposits, and Honolulu Volcanics) (Sherrod et al 2007). Average post-calibration hydraulic conductivities within each subzone defined for the unconsolidated caprock unit are 198.6 m/d (eastern), 88.8 m d −1 (middle), 64.2 m d −1 (western-fill), 87.1 m d −1 (western-beach deposits), 184.7 m d −1 (western-lagoon and reef deposits), and 157.0 m d −1 (western-Honolulu Volcanics). Tidal efficiencies are expressed as increases in piezometric head considering the tidal half amplitude (h 0 ) as follows: h(x)/h 0 =e −0.002x (eastern), h(x)/h 0 =e −0.005x (middle), h(x)/h 0 =0.55e −0.0007x (western-fill), h(x)/h 0 =0.43e −0.001x (western-beach deposits), h(x)/h 0 =0.44e −0.0004x (western-lagoon and reef deposits), and h(x)/h 0 =0.43e −0.0002x (western-Honolulu Volcanics) in which h(x) is the increase in piezometric head (m), and x is the distance from the shoreline (m). Based on the estimated diffusivities for these units, the computed tidal efficiencies are reasonable and consistent with the hydraulic properties. Similar to Habel et al (2017), the analytical solutions were applied within the respective boundaries of each of the six subzones to a raster grid as a function of the distance of each grid cell to the modeled coastline. Raster values representing tidal efficiency were calculated by setting h 0 to 0.33 (the MHHW tide elevation above the MSL datum in meters), which were subsequently summed with water-table raster data from the model output to generate the tidally influenced watertable height considering the MHHW tide stage.

Test of model
As part of cross-validation analysis, error statistics were calculated for the 3D numerical and hydrostatic methods based on the residual values between simulated and observed groundwater levels representing MSL and MHHW scenarios under current conditions (table 1).
In simulating MSL, the RMSE and systematic error (bias) of hydrostatic residuals reveal profoundly low estimates of groundwater level. This finding is reinforced by the negative maximum residual, indicating that at every comparison point, the method underestimates groundwater level. Performance improves markedly when simulating MHHW. This is evident from the significant reduction in RMSE and bias; however, the bias remains negative indicating overall underestimation of groundwater level. As expected, the RMSE representing the 3D numerical model is similar to that calculated in Habel et al (2017) considering both tidal scenarios.
Simulated groundwater levels are illustrated in cross-section along shore-normal transects (figure 3). Crosssection locations represent regions that feature different hydraulic gradients with distinctive subsurface geologies. Transect A represents a calcium carbonate platform comprising mainly Pleistocene skeletal limestone; Transect B represents alluvium and fill; Transect C represents limestone and fill and rises relatively abruptly in elevation. Results indicate that, although the hydrostatic simulations do not thoroughly align with the more robust model, the RMSE overlaps with results of the 3D numerical simulations (other than on transect C landward of approximately 1500 m). This is true for both evaluated tide scenarios, however, for the hydrostatic simulations the RMSE of the MHHW scenario is half that of the MSL scenario.

SLR flood simulations
Illustrations of GWI considering a 1-m SLR scenario are presented in figure 4. The hydrostatic method reproduces 65 percent (MSL) and 88 percent (MHHW) of the inundated area depicted by the 3D numerical method. In the MHHW case, 14 percent of the area inundated by the hydrostatic method lies outside of the 3D numerical-method inundation area. Thus, the total inundated area simulated using the hydrostatic method comes within 2 percent of the total inundated area simulated using the 3D numerical method. However, the flooded footprint of the two methods for the MHHW case differs by 26 percent as indicated by the uniquely flooded areas for the two methods (no overlap between two methods).  Figure 3. Composite water table simulations considering the influence of MSL (left column) and MHHW (right column) generated using hydrostatic (red) and 3D numerical (blue) simulations, along shore-normal transects (see figure 1 for location). Dotted lines denote RMSE departure from observed water levels. Lines A, B and C represent differing hydrogeologic conditions within the study area. Note that water levels along transect A dip at the inland extent due to the backshore presence of Pearl Harbor.

Discussion
As anticipated, our results reinforce 3D numerical modeling as the more robust of the two methods. Data assimilation by the 3D method provides for better representation of observed water levels. However, when referenced to MHHW, we find the hydrostatic approach produces simulations that are usefully accurate, with specific caveats. Improved hydrostatic simulation of MHHW relative to MSL results from two offsetting and unrealistic assumptions inherent in the method. These are (1) that an aquifer has a hydraulic gradient equal to zero, and (2) that tidal signals do not attenuate as they move through an aquifer. These assumptions produce errors that are oppositely sensed. This explains why the hydrostatic method underestimates groundwater elevations (RMSE of 40 cm, Bias of -33 cm) in the MSL scenario, as only the negative bias of the first assumption is introduced. This results from the fact that in most areas in the MSL case the actual groundwater hydraulic gradient is oriented toward the coast.
Commonly used hydrostatic simulations conveniently and fortuitously are referenced to local MHHW (i.e., NOAA SLR Viewer). When referenced to MHHW, the hydrostatic method can produce reasonable estimations of groundwater elevation in low-lying coastal regions, especially given the limited effort required for model construction. Municipalities can employ the method as a first-cut approach towards revealing vulnerabilities to GWI.
However, the two methods produce localized differences in flood simulation (figure 4) that result from the ability of the 3D method to capture unique hydrological conditions (i.e., recharge, and conductivity) that influence tidal efficiency and head. These differences illustrate the inability of the hydrostatic method to produce high-quality simulations that can be used as the basis for fine-scale decision making. Thus, we do not recommend use of the hydrostatic method alone for such endeavors. Rather, we advise it be used as an indicator of exposure to GWI at the municipal scale that could be used to inform decisions of whether methods of greater accuracy and precision are necessary. We recognize these findings apply specifically to the Honolulu area, since tidal ranges, topography, and hydraulic gradients vary regionally. However, because the minimum elevation of coastal groundwater generally exceeds that of local mean sea level in coastal regions (Turner et al 1997), it is reasonable to assume that where a coastal plain aquifer exists, the hydrostatic method (specifically considering a MSL tide stage) provides, at worst, a minimum estimate of groundwater elevation, and in turn GWI.
Use of the 3D numerical approach is more appropriate when (a) modeling coastal regions that feature particularly complicated conditions such as those that host extensive extraction/injection wells, (b) conducting modeling efforts that consider specific tidal scenarios (i.e. lower stages of the tide, extreme high tide); or (c) developing engineering techniques to mitigate flooding from GWI (i.e., implementation of extraction wells). We also note that, although the 3D numerical approach is more rigorous and widely applicable than the hydrostatic approach, the numerical approach may be unreliable if sufficient data are not available to constrain and evaluate model performance.
We also recognize that neither modeling approach presented here simulates dynamical coastal processes ( i.e., coastal erosion, sediment accretion or changes in land cover (Lentz et al 2015, Anderson et al 2018) that drive evolution of the landscape as sea level rises (FitzGerald et al 2008). Hence, our conclusions are most appropriately applied to regions, or environments that are less impacted by dynamical coastal processes (i.e. heavily developed shorelines that have been structurally hardened).

Summary and conclusion
Numerous coastal municipalities around the world face impacts from SLR flooding. The impacts are wide ranging and include disruptions in daily commerce, progressive failure of critical infrastructure, and intensified socio-economic burdens. In an effort to manage SLR impacts, informed, adaptive management is crucial and necessitates specific consideration of the various components of flooding including GWI. The GWI component is often overlooked in vulnerability studies, yet it is arguably the more challenging to manage as it includes complete saturation of the ground that is difficult to mitigate. This type of flooding can evade coastal armoring ( i.e., seawalls, revetments, levees) and overwhelm traditional drainage conveyances, rendering them ineffective.
To spur consideration of GWI in policy and planning, a data-intensive 3D numerical method was developed by Habel et al (2017) to specifically simulate SLR induced GWI; however, its accessibility is limited by data requirements to produce robust simulations. Here we investigate applicability of the more simple and accessible hydrostatic method in simulating GWI. The hydrostatic method is commonly used to produce flood simulations considering a direct marine source; however its applicability towards simulating GWI had not previously been explored.
Comparison of the hydrostatic method to 3D numerical modeling reveals each method's ability to replicate present day groundwater levels at MSL and MHHW stages of the tide, and similarities of GWI simulations considering 1 m SLR. For Honolulu the hydrostatic method produces groundwater level and GWI simulations that are comparable to the more physically based method, specifically when referencing the local MHHW tide stage (generally typical of flood studies). Hydrostatic simulations produce a RMSE of 20 cm during the MHHW tide stage, compared to 10 cm produced by the 3D numerical method. Further, hydrostatic simulation of GWI in a scenario of 1 m SLR at MHHW reproduces 88% of the inundated area simulated using the 3D numerical method. However, because neither method has been designed to simulate dynamic landscape changes their use should be limited to settings or environments that are less impacted by dynamic coastal processes that accompany change as a result of SLR (i.e. regions that host widespread coastal armoring).
Though use of data-assimilating numerical modeling methods are more appropriate in cases where high accuracy simulations are necessary, we find that use of the hydrostatic method (specifically when referencing the local MHHW tide stage) is suitably accurate as a first-cut approach in identifying municipal vulnerabilities to GWI. which is sponsored by the University of Hawaii Sea Grant College Program, SOEST, under Institutional Grant