A density-dependent multi-species model to assess groundwater flow and nutrient transport in the coastal Keauhou aquifer, Hawai‘i, USA

Fresh groundwater is a critical resource supporting coastal ecosystems that rely on low-salinity, nutrient-rich groundwater discharge. This resource, however, is subject to contamination from point- and nonpoint-sources such as on-site sewage disposal systems (OSDS) and urban developments. Thus, the significance of flow and transport processes near the coastline due to density effects and water circulation in a complex hydrogeologic system was investigated. A three-dimensional, density-dependent groundwater model was developed for the Keauhou basal aquifer (Hawai‘i Island, USA), where hydraulic head, salinity, nutrient concentrations, and submarine spring flux rates were used as calibration variables to best constrain parameters and produce a comprehensive aquifer management tool. In contrast, a freshwater-only model failed to properly simulate nutrient transport, despite the reasonable success in calibrating hydraulic head measurements. An unrealistic value for hydraulic conductivity was necessary for freshwater-only calibration, proving that hydraulic conductivity is a process-based variable (i.e., depends on model conceptualization and the simulated processes). The density-dependent model was applied to assess relative contaminant source contributions, and to evaluate aquifer response concerning water levels and quality due to changing environmental conditions. Nutrients detected in the aquifer are primarily sourced from OSDS, which was supported by a nitrogen isotope mixing model. Additionally, effects of sea-level rise emphasized the complexity of the study site and the importance of model boundaries. While the model is developed and applied for West Hawai‘i, the adapted approaches and procedures and research findings are applicable to other coastal aquifers.


Introduction
Groundwater is a critical natural resource required for daily human use such as potable water supply, industrial uses, and agricultural irrigation (Kemper 2004;Alley 2006). The Keauhou region, located on the west side of Hawai'i Island (USA), has approximately 97% of its freshwater demand supplied by groundwater resources pumped from the Keauhou aquifer system (Fukunaga 2017). As Hawai'i's population increases, future land development and higher water demand are inevitable, along with an expected overall decline in recharge rates due to projected climate change conditions (Elison-Timm et al. 2015). These combined factors will further exacerbate saltwater intrusion (Ferguson and Gleeson 2012), thus emphasizing the need to understand how future circumstances may affect groundwater resources.
Additionally, submarine groundwater discharge (SGD) transports human-derived nutrients to the coast, discharging as point sources (submarine springs) or non-point seepage (diffuse) sources (e.g., Burnett et al. 2003;Moore 2010;Parra et al. 2014). Natural and anthropogenically derived nutrients such as nitrate (NO 3 − ) and phosphate (PO 4 3− ), have the potential to negatively impact drinking water and the biota of coastal ecosystems, which serve as important environmental, economic, food, and cultural resources (Duarte et al. 2010). Excess nutrients often originate from different types of land use and land cover (LULC) and wastewater infrastructuree.g., on-site sewage disposal systems (OSDS); e.g., Knee et al. (2010); Shuler et al. (2017); Delevaux et al. (2018). Hawai'i encompasses several types of LULC, including natural, agricultural, and varying degrees of urbanization . Continuous urban, suburban, rural housing and commercial development alters groundwater quality through wastewater and landscape runoff and is of growing concern, especially where groundwater directly discharges to coastal waters (Knee et al. 2010). In many coastal areas, OSDS associated with development are considered one of the major sources of nutrient contamination (e.g., Lapointe et al. 1990;Harris 1995;Reay 2004). The majority of these OSDS are cesspools, where wastewater effluent is not treated before entering the groundwater system . Most of Hawai'i's OSDS are located near the shoreline , and hence overlooking freshwater-saltwater interactions can lead to inaccurate interpretations of resulting contaminant loads and consequently poorly inform water management decisions.
Multiple studies across Hawai'i demonstrated that elevated nitrogen is the most common coastal water impairment problem (Howarth and Marino 2006). Excess amounts of nutrients such as NO 3 − , are of primary concern because they can have negative impacts on human health and can significantly affect coastal ecosystem productivity and health (Kendall and Aravena 2000;Duarte et al. 2010). Besides natural sources and agricultural and lawn runoff, wastewater is the most frequently implicated source of excessive nitrogen Street et al. 2008). Stable isotopes of nitrate (δ 15 N and δ 18 O of NO 3 − ) are commonly used to distinguish nitrogen derived from various sources (e.g., Kendall and Aravena 2000). It can sometimes be difficult, however, to identify the exact nitrogen source due to the wide range and overlapping nature of these isotopic endmembers. An integrative, modelbased approach can quantify relative nutrient contributions from various sources and identify essential areas of remediation for resource managers (e.g., Shuler et al. 2017).
Water resource managers often turn to numerical models to better understand the groundwater system and predict how aquifers will respond to future climate and LULC scenarios. Outside of Hawai'i, various numerical models have been developed to investigate groundwater contamination from LULC and OSDS in coastal aquifer systems using a freshwater-only approach (e.g., Pang et al. 2006;Shuler et al. 2017;Walter 2008) and a density-dependent approach Murgulet and Tick 2016). Generally, groundwater flow and contamination in Hawai'i are studied using the US Geological Survey (USGS) suite of models (US Geological Survey 2021;El-Kadi and Moncur 2006). Unfortunately, there are only a few Hawai'i-based models developed to investigate LULC and OSDS contamination. Delevaux et al. (2018) and Amato et al. (2020) developed freshwater models to examine LULC and OSDS-derived contamination transport, while Habel et al. (2017) examined groundwater inundation as a result of sea level rise. In comparison, studies such as Souza (1987, 1998) and Gingerich and Voss (2005), have used density-dependent models to examine the freshwater-saltwater transition zone but did not assess OSDS contamination. To avoid complications in solving nonlinear equations (e.g., Holzbecher 1998) and increased data needs, current management of Hawai'i's water quality relies on a freshwater-only modeling approach (e.g., Whittier and El-Kadi 2014). This approach has several limiting factors, including the failure to account for (1) dynamics of the freshwater-saltwater interface and forming a brackish transition zone; (2) inland saltwater influxes required for accurate assessment of groundwater sustainability; (3) accurate parameters based on salinity measurements; (4) potential sea level rise and saltwater intrusion; and (5) an accurate flow field that controls chemical transport (Mezzacapo et al. 2020).
Based on these limiting factors, it is hypothesized that a freshwater-only model fails to properly simulate nearshore solute transport. The objective of this study was thus to investigate the significance of flow and transport processes near the coastline due to density effects and water circulation in a complex hydrogeologic system. This study aimed to ameliorate earlier model approaches and demonstrate the applicability of the model as an aquifer assessment tool. To that end, a three-dimensional (3D) density-dependent multi-species model was developed for the Keauhou basal aquifer and calibrated with all available data, including hydraulic head, salinity, submarine spring flux rates, and nutrient concentrations. To avoid overfitting, the number of calibration parameters was minimized by utilizing independently estimated values and critical input variables were investigated through sensitivity analysis. Further, δ 15 N was used for independent model validation. The relative nutrient contributions predicted by the model were integrated with δ 15 N measurements and an isotope mixing equation to confirm point-and non-point source contamination. A similar approach was utilized for Tutuila, another tropical volcanic island, by Shuler et al. (2017), but has not yet been done in this study area. The study also assessed suitability of a freshwater-only model in simulating the transport of nutrients in the aquifer. Finally, the density-dependent model was applied to test cases, specifically urban development, climate change, and sea level rise.

Study area
Hawai'i Island is the largest island (10,464 km 2 ) in the Hawaiian archipelago and was built by five shield volcanoes: Kohala, Mauna Kea, Hualālai, Mauna Loa, and Kīlauea (Izuka 2011). Basaltic lava flows from these five volcanoes covered the island's surface less than 1 million years ago, making Hawai'i Island the youngest island in the chain (Clague and Dalrymple 1987). The eruption of each volcano was partially simultaneous with its neighboring volcano, so it is likely that lava flows from neighboring volcanoes interfinger with each other, creating complex geologic contacts (Wolfe et al. 1997;Izuka et al. 2018). Basaltic lava typically flows as either pāhoehoe or 'a'ā, which differ in viscosity and texture (e.g., DePaolo and Stolper 1996;Lau and Mink 2006;Izuka et al. 2018). For sequences of pāhoehoe, groundwater primarily travels through spaces, tubes, and vesicles that are formed between the multiple layers of lava flows, resulting in higher porosity and permeability. 'A'ā are comprised of two general components-clinker beds and massive cores. The 'a'ā clinker beds typically demonstrate higher porosity and permeability while the dense 'a'ā cores demonstrate lower porosity and permeability and thus impede groundwater flow across layers (Izuka et al. 2018). Because it is common to see alternating pāhoehoe and 'a'ā flows in cross section, the hydrogeologic parameters typically represent both flow types (Lau and Mink 2006). Located in the northern Kona district, the southern flank of Hualālai is defined as the Keauhou aquifer (426 km 2 ; Fig. 1). The Keauhou aquifer is geologically divided into a basal aquifer and a high-level aquifer with a discontinuous head profile characterized by a near abrupt change. Water levels near the coast are approximately 1 m above mean sea level (msl) but dramatically rise to more than 100 m above msl approximately 4-8 km inland (Fig. 1). It is believed that a subsurface geologic structure reduces groundwater flow from the high-level aquifer to the coast (Oki 1999). Efforts to characterize this structure include modeling of gravity survey data that indicate the presence of a dense substructure under Hualālai, located approximately 4 km southwest of the eruptive vents that map Hualālai's southeast rift zone (Kauahikaua et al. 2000;Flinders et al. 2013). Although the geologic structure has been conceptually hypothesized (Oki 1999, their Fig. 6), the exact nature of the structure is still not well defined, which is beyond the scope of this study. Another significant, yet undefined, geologic feature acts as a confining unit, creating a second, deep freshwater body below the primary basal aquifer (Thomas et al. 1996;Attias et al. 2020). Since this study focuses on the shallow groundwater, however, no attempt was made to simulate such a feature. To avoid such uncertainties, this study instead focuses on the Keauhou basal aquifer (150 km 2 ), with head levels less than 5 m and the eastern boundary coinciding with the geologic structure bordering the high-level aquifer (Fig. 1). The western boundary follows bathymetry and extends beyond the coastline to accurately account for processes at the land-ocean interface.

Water levels
Considering the scarcity of current measurements and necessity of spatially distributed data, initial pre-development water levels were used in the model calibration for a total of 54 measurements. Static, initial water levels (one per well) were obtained from the Commission on Water Resource Management (CWRM), where the majority of these levels were measured at the time of drilling, which occurred between 1944 and 2008 (Fig. 2d). While water levels change over time and available values at other wells do not necessarily accurately represent current conditions, initial water level data from four monitoring wells fell within the ranges of continuous measurements sampled on a quarterly basis (Fig. S1 of the electronic supplementary material (ESM). Further, the inclusion of all initial water level measurements yielded better results with a more robust dataset.

Nutrient and δ 15 N data
Median salinity, NO 3 − + NO 2 − and PO 4 3− concentrations from 20 wells (Tachera 2021) were used for model calibration (Table S1.2 of the ESM). Production well samples were collected three to four times over a year, thus median values were used to represent annual concentrations. Thirteen salinity concentrations and 19 NO 3 − + NO 2 − and PO 4 3− concentrations were collected at identified submarine springs by H. Dulai (University of Hawai'i at Mānoa, unpublished data, 2019) and utilized in model calibration ( Fig. 2d; Table S1.3 of the ESM). Hawaiian aquifers are well oxygenated, exhibiting dissolved oxygen saturations >90% (Richardson et al. 2017). It is therefore assumed that dissolved nitrogen species quickly oxidize to NO 3 − (Whittier and El-Kadi 2009). Because NO 3 − + NO 2 − and PO 4 3− (further referred to as N and P, respectively) were the major chemical species of nitrogen and phosphorus, these were used during model calibration. Concentrations of NO 3 − represent >70% of TN in the Keauhou aquifer, thus the error introduced by using NO 3 − + NO 2 − values for model calibration is assumed to be minimal.
All nutrient and stable nitrogen isotope (δ 15 N) samples were collected between November 2017 and March 2019 (Table S1.1 of the ESM). Well samples were collected from port heads in 10-L cubitainers and chilled in coolers when transported to the laboratory. Nutrient samples were transferred from the cubitainers to acid-washed 250-ml HDPE bottles and immediately chilled until analysis. The samples were analyzed for a suite of nutrients, which include ammonium (  Briefly, using a "reservoir" model, the recharge was calculated with rainfall (Giambelluca et al. 1986), fog interception (Juvik and Ekern 1978), and irrigation as inputs, and total evapotranspiration and runoff as outputs.

Current land use and land cover (LULC)
In areas of urban development, this model considered LULC and OSDS as primary nutrient sources. LULC can influence the groundwater concentrations of N and P due to wastewater and nutrient inputs to lawns and agricultural systems. In this model, LULC was divided into three main categories: natural (background), golf courses, and urban development (Fig. 2f). Natural N and P concentrations were assigned based on average measured concentrations from pumping wells across Keauhou (1.0 and 0.1 mg/L, respectively). Golf course locations were extracted from the Hawai'i Statewide GIS Program (State of Hawai'i 2021) and assigned predetermined N and P concentrations of 7.59 and 0.54 mg/L, respectively (Delevaux et al. 2018). Urban development was further divided into three levels of development intensity (low, medium, and high), as defined by the 2005 Hawai'i Coastal Change Analysis Program (National Oceanic and Atmospheric Administration 2006). In this study, urban development was defined by the land's green space. Structures associated with urban development were treated separately through calculation of nutrient inputs from OSDS. Nutrient loads associated with residential wastewater from OSDS in urban developed lands were added separately to the model as a point-source coverage-further explained in section 'Onsite sewage disposal systems (OSDS) and wastewater treatment plants (WWTP)'. Delevaux et al. (2018) determined that the N and P concentrations leaching from a typical lawn are 0.2 and 0.01 mg/L, respectively. The low, medium, and high intensity developments correspond to 65, 35, and 10% of green space per land parcel, respectively. Therefore, low intensity development areas, with a larger proportion of lawns, leach more N and P into the ground via lawn maintenance. In contrast, high intensity development areas, with a smaller proportion of lawns, leach less N and P into the ground. To account for the development intensity of each land parcel, typical lawn concentrations were scaled according to the proportion of lawns for each development intensity and added to the background concentrations (Table 1). These LULC concentrations were intersected with recharge and applied to the model.

On-site sewage disposal systems (OSDS) and wastewater treatment plants (WWTP)
The OSDS point-source coverage was obtained from a statewide OSDS census database (State of Hawai'i 2021;Whittier and El-Kadi 2014). Within the study area, there are 1,094 class I (soil) units, 77 class II (septic tank to seepage pit) units, 16 class III (aerobic treatment to seepage pit) units, and 6,251 class IV (cesspool) units (Table S2 of the ESM). Nutrient mass loads for each OSDS point are reported in units of kg/ day and were determined by applying typical nutrient concentrations to the estimated effluent rate (Table S2 of the ESM; Whittier and El-Kadi 2014). This OSDS coverage has a much finer resolution than the model grid, so to properly simulate the total (rather than average) mass load entering the groundwater system, all loads from OSDS points that fall within a single grid cell were integrated (Fig. 2g). That summed value was applied to the top layer of that respective cell in the model. With this method, 250 grid cells represent the OSDS compared to the 6,593 actual units, but the total mass load of 890 kg/day was maintained. The Kealakehe WWTP treats wastewater in aerated lagoons before disposing the effluent in a 930-m 2 disposal percolation basin (Wilson Okamoto Corporation 2019; Fig. 2g). Currently, the basin has an average flow rate of 6,435 m 3 /day with N and P mass loads of 144 and 45 kg/day, respectively (Wilson Okamoto Corporation 2019).

Conceptual model
Rather than modeling the entire Keauhou aquifer, only the basal area was selected as the focus for this study (Figs. 1 and 2a) and the high-level Keauhou aquifer was excluded due to the following factors: 1. Water level observations from the high-level aquifer needed for accurate calibration are inadequate. 2. In contrast to the basal aquifer, there is less concern for increased salinity under extreme conditions due to the thick freshwater layer created by the high-level water. 3. Similarly, nutrient contamination threats are less serious than in the basal aquifer. 4. There is uncertainty regarding the inland (eastern) boundaries of the Keauhou aquifer ( Fig. 1). Such boundaries are based on surficial features (Commission on Water Resource Management 2019), which do not necessarily coincide with subsurface divides. Therefore, potentially significant inter-aquifer groundwater flow may occur and thus cannot be treated as no-flow boundaries (Fackrell et al. 2020). As an additional benefit, modeling the smaller area reduces computational expenses associated with modeling density-dependent water flow and transport. With the exclusion of the high-level aquifer, the eastern model boundary of the basal aquifer was assigned a specified flow condition and the flow rate was estimated by creating a regional, freshwater-only MODFLOW model for West Hawai'i (Fig. 2b). Since the eastern boundary of this regional area generally follows the rift zones of Mauna Kea and Mauna Loa, groundwater is not believed to flow from the windward to leeward side of the island (Izuka et al. 2018), and thus the boundary can safely be assigned a no-flow condition. Such a regional model would eliminate the concern about interaquifer flow. The geologic structure was simulated as a horizontal flow barrier in the regional model and was reasonably calibrated to estimate hydraulic conductivity (K h ) and the barrier conductance. MODFLOW's zone-budget option was used to estimate the flow into the eastern boundary of the basal aquifer area. The flows through the northern and southern boundaries were negligible and therefore it was reasonable to treat them as no-flow boundaries.
Submarine springs are dominated by brackish spring flows along the coastline of the study site (Johnson et al. 2008;Duarte et al. 2010;Kelly et al. 2019). In contrast to diffuse SGD, these point-source springs are caused by preferential flow (Taniguchi et al. 2019). Explicitly modeling preferential flow is impractical due to the absence of detailed aquifer information, mostly related to identification of fractures or conduits contributing to preferential flow. As an approximation, springs were simulated as point drains, or head-dependent sinks, at their respective coordinate locations along the coastline. A similar approach was adopted for volcanic aquifers by Oki (2005) and El-Kadi et al. (2014). The approach is motivated by studies that indicated that discharge from springs is directly related to the aquifer's head gradient (Taniguchi et al. 2002;Duarte et al. 2010). In MODFLOW, the drain flow is estimated through Darcy's law, utilizing conductance and hydraulic head terms. In this study, with head levels roughly constant near the shoreline, it was reasonable to assume that conductance of each submarine spring was linearly proportional to the respective measured flow rate. Thus, each submarine spring point has a unique conductance value representing the varying degrees of preferential flow that can be encountered across this complex aquifer system. The bottom of the drain was assumed to be 0.5 m below msl, representing the rough depth from which the springs were sampled. The drains were thus located within the top layer of the model.

Numerical model
The study applied SEAWAT, a transient, multispecies density-dependent model (Langevin et al. 2008), which combines MODFLOW (Harbaugh et al. 2000) with MT3DMS (Zheng and Wang 1999). The software package interface GMS (Aquaveo 2021) was used in facilitating the modeling process. Quasi steady-state simulations were utilized by running the models for 30 years to reach near steady-state conditions. The model consists of 21 layers, 82 rows, and 36 columns (totaling 26,323 cells) and is rotated 10°counterclockwise. The model grid covers the Keauhou basal aquifer and extends out into the ocean. The grid's top elevation follows topography and bathymetry and the bottom of the grid has a flat elevation of 550 m below msl (Fig. 2c). With head levels less than 5 m above msl, and assuming validity of the Ghyben-Herzberg relation, the freshwater-saltwater interface is expected to be located approximately 200 m below msl. The additional model depth was included as a buffer to ensure the flow in the model area is not significantly affected by the bottom boundary condition. Because the bathymetry offshore of Keauhou is highly variable, the model grid extends perpendicularly offshore as far as 12,034 m and as close as 743 m. The bottom of the first layer is set to 1 m below msl to ensure dry cells are not produced if the water table drops below the cells' bottom. The remaining 20 layers vary gradually in thickness, where the upper layers are the thinnest (Fig. 2c). The horizontal grid spacing is a consistent 490 m × 490 m.

Calibration approach
Typical calibration approaches were used to estimate parameters of interest that gave the best match between measured and simulated head, salinity, nutrient concentrations, and submarine spring flux rates, with root-mean-square error (RMSE) as the criterion for assessing match quality. Available measurements were described earlier in section 'Available data'. The number of calibration parameters was minimized by utilizing independently estimated values to the extent of their availability. The calibrated parameters were horizontal hydraulic conductivity (K h ), vertical anisotropy (K h /K v ), porosity (φ), and longitudinal dispersivity (α L ), and their values were constrained based on available literature ranges. Parameters that were not calibrated were specific yield (S y ) and specific storage (S s ).
Two approaches were adopted for calibrating K h . In the first, a homogeneous K h was assumed, and nine trial-anderror iterations were performed for calibration, where only one parameter varied per iteration. In addition, sensitivity analyses were completed to quantify the effect of potential model parameter variations on the modeling results (Cacuci 2003) (see section S.3 of the ESM for details). In a second approach, a heterogeneous K h distribution was calibrated using the principal component geostatistical approach (PCGA) as a scalable inversion method for spatially distributed field characterization (Lee and Kitanidis 2014;Lee et al. 2016). Using both head levels and salinity measurements in wells for calibration, spatially distributed K h points were estimated. A prior statistical correlation structure, modeled by an exponential covariance function with the scale length of 2,000, 2,000, and 200 m in the x, y, and z directions, respectively, was determined by a Bayesian hyperparameter optimization approach (Kitanidis 1996) and used as a constraint for the calibration. An open-source Python package (Lee 2021) was linked to the USGS SEAWAT model and approximately 200 SEAWAT model runs were needed to obtain calibration results.

Parameter ranges
Few studies have computed K h values for West Hawai'i aquifers. Using specific capacity data, Rotzoll and El-Kadi (2008) computed K h values that range from 155 to 11,000 m/day within the Keauhou aquifer, and Oki (1999) simulated Kona's groundwater system by assigning a K h of 2,300 m/ day to the dike-free volcanic rocks. A K h /K v value of 200 was calculated by Souza and Voss (1987), which is consistent with lava flow orientation and fracture configuration, and is generally accepted for Hawai'i-based groundwater models (e.g., Voss and Souza 1987;Oki 1999). Values of φ for O'ahu's aquifers range from 0.05 to 0.5 due to the varying degrees of island weathering and age (Izuka 2011). Within that range, Oki (1999) and Izuka (2011) used a homogeneous value of 0.1 for Kona groundwater models. Estimates of α L were introduced by Oki (2005), Izuka (2011), and Glenn et al. (2013); however, the scale dependence of α L is emphasized by previous studies (e.g., Gelhar et al. 1992;Lee et al. 2018). Schulze-Makuch (2005) estimated that α L ranges from 0.1 to 100 m for basalts using a scale-dependent regression equation. Previous values of S y range from 0.04 to 0.055 (Gingerich and Voss 2005;Izuka 2011) and values of S s range from 2.3 × 10 −5 to 2.4 × 10 −5 L/m (Oki 1999;Izuka 2011).

Boundary conditions
The regional model estimated that approximately 470,000 m 3 of water enters the basal aquifer per day. That value was adjusted based on calibration to 388,399 m 3 /day and assigned to the top two layers of the eastern boundary (Fig. 2c). Because urban developments also exist within the high-level aquifer, the eastern boundary was assigned specified salinity, N, and P concentrations of 0.26 ppt, 1 mg/L, and 0.1 mg/L, respectively. These concentrations were the average of the median concentrations measured from pumping wells located within 1 km around the boundary. These include a combination of five high-level wells and five basal wells. The ocean floor was treated as a general head (or head dependent) boundary (Fig. 2c), with a conductance of 1 m 2 /day/m 2 and a salinity of 35 ppt. The conductance was estimated based on coastal leakance ranges proposed by Oki (1999). All other boundaries were assigned a no-flow condition. δ 15 N model validation There are numerous ways nutrients are introduced to the groundwater system in the study site, including pointsources such as OSDS and WWTP, and diffusive nonpointsources such as LULC (Kendall and Aravena 2000). Because multiple nutrient inputs are located within Keauhou and because their δ 15 N and δ 18 O of NO 3 − composition ranges overlap, it is difficult to identify the relative contribution of each nitrogen source. Relative N contributions from OSDS, WWTP, LULC, and eastern boundary condition (BC) were estimated by assigning different N species to these sources. The model is equipped to simultaneously simulate these as distinct species. Following the principle of superposition, the total N concentration was estimated by summing these individual concentrations (Reilly et al. 1984).
To validate the modeled source estimates, δ 15 N in dissolved NO 3 − were used to identify nitrogen sources in the aquifer; however, the analysis is complicated by the fact that each source exhibits overlapping ranges in isotopic composition. Abaya et al. (2018) measured δ 15 N from various N sources in Puakō, which is located on the coast north of the study area. The study estimated natural soil δ 15 N ranges from approximately 1-4‰. While most δ 15 N samples collected within the Keauhou basal aquifer fall within the natural soil δ 15 N range (e.g., Kendall and Aravena 2000), elevated NO 3 − measurements for some samples suggest N is added by other inputs as well. To test this, a simple isotope mixing equation was used to compute δ 15 N values for wells and submarine springs based on simulated relative contributions and δ 15 N endmember assumptions. Since mixing ratios are used for stable isotopes of dissolved components, δ 15 N needs to be weighted by the simulated concentrations in the mixing model (Hunt and Rosa 2009 where δ 15 N model is the computed δ 15 N at wells and submarine springs, f is the simulated fraction of each source component at a well or submarine spring location, C is the simulated concentration in the component, and dN is the δ 15 N endmember of each source component. The δ 1 5 N endmembers used for OSDS, WWTP, LULC, and BC were 7‰ (Wiegner et al. 2016), 1‰ (Abaya et al. 2018), 20‰ (Hunt 2014), and 2‰ (Hunt 2014), respectively. Single endmember values were selected to obtain the best match, but as previously mentioned, each source ranges in δ 15 N. Wiegner et al. (2016) summarized that δ 15 N for sewage in Hawai'i ranges from 7 to 20‰, while Abaya et al. (2018) narrowed the range to 9-11‰ in Puakō, Hawai'i Island. Abaya et al. (2018) reported that values for fertilizers range from −2 to 0‰ while soil ranged from 1 to 4‰. Since the modeled LULC coverage contains natural land and urbanized land, both ranges were considered. Hunt (2014) reported that the treated wastewater from Kealakehe had a δ 15 N measurement of 12‰ and a water sample from a well within the Kealakehe wastewater plume had a δ 15 N measurement of 20‰. Hunt (2014) reported that the Honokōhau well, which is a high-level pumping well in Keauhou, had a δ 15 N measurement of 2‰. This is similar to this study's high-level samples, which ranged from 1 to 2‰.

Freshwater-only model
As previously stated, models overlooking salinity changes in nearshore aquifers are commonly used to assess both water levels and quality. To test applicability of such an approach, a freshwater-only model was simulated to compare to the density-dependent results. The model was easily set by converting SEAWAT to run MODFLOW and MT3DMS in sequence, with MODFLOW set to run in a steady-state condition. The model retained the boundary conditions, coverages, and various input parameters previously described, except K h , which was recalibrated to obtain a reasonable match for hydraulic head under freshwater-only settings.

Model applications
To demonstrate the capability of the density-dependent model as a management tool addressing nutrient and salinity distribution in the basal aquifer, the model was tested with example applications. This study therefore focused on two different aspects that are important for water resource management in West Hawai'i: (1) aquifer salinity and storage changes as a response to future development and climate change, and (2) aquifer response to sea level rise.

Aquifer response to urban development and climate change
Accounting for future water demand is crucial in ensuring that the aquifer will not be over-pumped or subjected to water quality deterioration with anticipated urban expansion. As an example use of the model, a future urban development application was created by altering the LULC map described in section 'Current land use and land cover (LULC)' (Fig. 2f). The Kona Water Use Development Plan identified areas where new development is expected to occur based on allocated water permits (Fukunaga 2017). These new development areas were designated anticipated water demands, where a demand of 1.5 m 3 /day corresponds to one residential unit (Fukunaga 2017). The water demand, however, does not necessarily correlate with a particular development intensity level. In the majority of cases, new development occurred in areas characterized as natural background, where there is no current development. These areas were therefore assigned nutrient concentrations associated with medium intensity development, the most common development intensity in the area. In some areas, however, new development involved adding units to previously developed land. In this case, anticipated water demands for the new parcels (Fukunaga 2017) were used to determine whether the development intensity would significantly change with the additional development. Since one residential unit does not appear to be a significant addition to an already medium-or high-density area, development intensity was assumed to remain the same where anticipated water demand was 1.5 m 3 /day (adding one or less new home) but was assumed to increase from medium intensity to high intensity where anticipated water demand was >1.5 m 3 /day (adding more than one new home). Based on current legal regulations in Hawai'i (Act 125), all cesspools are required to be replaced by 2050, so future development will need to be connected to sewer lines or replaced with septic or aerobic treatment (ATU) systems. For simplification purposes, it was assumed that the first option would be adopted. It is also assumed that 50% of the total water demand from new development will be used outdoors (Board of Water Supply 2021), so the remaining 50% will be routed to the Kealakehe WWTP via sewer lines. With a total future water demand of approximately 64,080 m 3 /day across Keauhou (Fukunaga 2017), 32,040 m 3 /day will be directed to the Kealakehe WWTP, which results in 724 and 225 kg/day of additional N and P mass load, respectively. Combining this with the already existing mass load, the disposal percolation basin was simulated to receive total N and P mass loads of 868 and 270 kg/day, respectively. Concurrently, the total future water demand is assumed to be pumped from the high-level wells in the Keauhou aquifer, so 64,080 m 3 /day was removed from the eastern boundary condition, thus reducing the influx to 324,319 m 3 /day. . Based on interpolated maps, the Keauhou area may experience approximately 20% reduction in rainfall following the RCP 8.5 scenario (Elison-Timm et al. 2015). It is outside the scope of this project to create a new recharge coverage by computing a water budget under climate change conditions; therefore, to account for the overall reduction in rainfall expected at the end of the century, the current recharge computed by the USGS (Engott 2011) was simply reduced by 20% across the entire basal aquifer.
Aquifer response to sea level rise As a coastal aquifer, impacts of future sea level rise are of great concern. Global sea level is predicted to rise approximately 1 m by the end of the century, if greenhouse gas emissions continue to rise at the current rate (e.g., Pfeffer et al. 2008;Nichols and Cazenave 2010;Hawai'i Climate Change Mitigation and Adaptation Commission 2017). To simulate future sea level rise, the offshore head stage was set to 1 m above msl to assess how increasing levels could potentially impact submarine spring flux rate and salinity.

Density-dependent model calibration
Utilizing a heterogeneous K h distribution, as described in section 'Calibration approach', did not lead to a significant improvement in results (section S.4 in the ESM). Accuracy of salinity calibration in wells improved but was accompanied with a slightly reduced accuracy of hydraulic head calibration; additionally, the accuracy of N and P predictions worsened. It was concluded that such an effort is only useful if appropriate representation of various processes is included in the conceptual model. Under scarcity of data, utilizing a simplified K h distribution is preferred, considering the complications and computational expenses associated with developing a more involved K h distribution. The results presented from here forward are all for a homogeneous K h applied to the entire aquifer. The final calibrated parameter values used in the Keauhou basal aquifer model are provided in Table 2.
Head, salinity, N and P concentration calibration results are shown in Fig. 3. Head was reasonably calibrated with a relatively low RMSE (<0.5 m; Fig. 3a). Under density-dependent conditions, the calibration process is complicated by the nonlinear relationship between the hydraulic head and salinity. Attempts to improve on the accuracy of salinity calibrations resulted in less precise head calibrations. The reasonable accuracy of head prediction, nonetheless, implies validity of the continuum approach via Darcy's law (e.g., Bear 1972) for simulating the hydraulic head for this complex domain. However, lower accuracy and the scatter regarding salinity and nutrients are likely attributed to the complexities in geological characteristics of the field such as the potential existence of preferential flow, which would naturally affect various simulation processes. The simulated salinity in submarine springs exhibit lower accuracy due to the nature of sampling location. Submarine springs are simulated along the shoreline and are positioned at the interface between fresh groundwater exiting the aquifer and saline ocean water moving inland. Thus, most of the submarine springs have a simulated salinity confined within a narrow range in comparison to the wide range of salinities measured in the field (Fig. 3b). This mismatch is possibly caused by the model's inability to simulate a very precise spring point location. The field-based submarine spring salinity measurements were assessed from the freshest terrestrial water that could be detected at the coast; thus, they are expected to be biased towards lower values. In contrast, the modeled salinity was averaged over an entire cell, and therefore falls within a mid-salinity range.
The N simulations produce reasonable calibration results (Fig. 3c,d). In contrast, many of the observed P measurements fall within a small range that is not well represented in the model results (Fig. 3e,f). The imprecise results for P should be expected considering that the assumption of a conservative behavior is not valid (Delevaux et al. 2018); however, the results can be useful while applying the model for screening purposes. In contrast to salinity, the relatively better correlation between simulated and measured N and P for submarine springs is due to the fact that these species are primarily sourced from terrestrial sources located in freshwater zones and not drastically influenced by marine processes at the land-ocean interface.
Modeling submarine springs as drain points with varying conductance produced accurate results (RMSE = 1,237 m 3 / day, r 2 = 0.98). Under such conditions, the model simulated 107,673 m 3 /day of groundwater discharging at point-source submarine springs, closely matching the observed values. The  Voss and Souza (1987) b Oki (1999) c Izuka (2011) d Hunt (1996) model also simulated an additional 1,230,000 m 3 /day of fresh and saline water diffusively discharging along the coast and offshore. However, accurately simulated flux rates combined with lower accuracy of salinity, N, and P concentrations emphasize the fact that a continuum approach for flow modeling can be acceptable for preferential flow but fails to simulate solute transport.

Validity of the freshwater-only model
Although comparable head accuracy was obtained for the freshwater-only model and the density-dependent model (Fig. 3a), a significantly lower K h of 300 m/day was required to calibrate the hydraulic head. Compared to a value of 2,500 m/day for the density-dependent model, such a value is not realistic for a young basalt, as documented in previous modeling and other studies (Oki 1999;Izuka 2011). Hence, such a value can simply be considered as a fitting parameter that is not physically based. It is of interest to note that K h appears to be a process-based parameter, thus it depends on the nature of the conceptual model and processes simulated. Saltwater circulation clearly affects the estimated effective aquifer thickness, which cannot be accurately estimated by a freshwater-only model. RMSE of submarine spring flux also increased from 1,237 m/day to 3,273 m/day. Although the submarine spring drains retained the same conductance values, the lower K h assigned to the freshwater-only model reduced overall flow fields (Fig. 4b) and resulted in lower discharge at the submarine springs. The overall accuracy of N and P concentrations are comparable for the two models at wells (Fig. 3c,e), in contrast to the values in the submarine springs (Fig. 3d,f). The freshwater-only model overestimates nutrient concentrations in the springs due to the fact that water circulation is not simulated in the nearshore area, where in reality, nutrientpoor saltwater intrudes into the aquifer. As should be expected, the influence of circulation diminishes far from the shoreline, leading to better accuracy for wells. Flow fields for the two models are compared in Fig. 4, showing a significant contrast in the nearshore area. It should be noted that salinity does not directly affect nutrient concentrations, considering that they are treated as conservative species. The case can be different for reactive chemicals such as contaminants of emerging concern (CECs), where salinity can directly play a role (e.g., Patel et al. 2019;McKenzie et al. 2020). It can be concluded that, not only does the freshwater-only model overlook salinity, but it also ignores the significant water circulation in the nearshore environment, leading to unacceptable errors in water quality assessment. Salinity is an important factor regarding water sustainability concerning different water uses, such as reducing estimates of aquifer sustainable yield and causing damage to anchialine ponds. Fresh groundwater in Hawai'i is considered to have less than 1,000 mg/L of dissolved solids (Yee 1987), while the USGS also classifies Fig. 3 Scatter plots of density-dependent and freshwater-only model calibration results for a head, b salinity, c N in wells, d N in submarine springs, e P in wells, f P in submarine springs. As it is not possible to simulate salinity using a freshwater-only model setup, part b only contains density-dependent results. Root-mean-square errors (RMSE) are included and carry the same units as the respective measurements water as "slightly saline" with concentration less than 3,000 mg/L (Swenson and Baldwin 1965). If it is assumed that salinity measurements directly correspond to dissolved solid concentrations, then 45% of the 20 wells used in the Keauhou basal aquifer model had salinity measurements that exceed the "slightly saline" water standard, while 75% of the simulated salinity measurements exceed the threshold. Such an assessment clearly indicates the need for a densitydependent model to accurately assess water sustainability and estimate aquifer sustainable yield.
Relative N source contributions and δ 15 N of NO 3

− model validation
Based on simulated N mass loadings entering the Keauhou basal aquifer, the majority of N is sourced from OSDS (54%), while urban development LULC contributed 10%, the Kealakehe WWTP contributed 9%, and groundwater applied at the eastern boundary contributed 26% (Fig. 5). It is clear that OSDS have the greatest impact on nutrient contribution to the groundwater system. Areas with urban development, where there are significant clusters of OSDS, experience plumes of higher nutrient concentrations, which travel with groundwater flow to the coast (Fig. 5b). The results can be useful for managing water quality issues such as efforts to meet Federal regulations regarding the total maximum daily load and identifying major contributor sources where upgrade efforts should be directed. For this study, it is important to recognize the significant nutrient mass load that is sourced from the eastern boundary condition (Fig. 5e), especially while considering the boundary uncertainties previously discussed. The mass load that enters through the boundary was dependent on the assigned nutrient concentrations and water flux.
Based on simulated N concentrations in wells and submarine springs, relative N source contributions were also estimated at each sample location (Fig. 6). Overall, minimal amounts of LULC-sourced N were detected in simulated samples, with an average contribution of 0.13 mg/L (13% of total concentration). The N sourced from the WWTP contributed a significant amount of N to a few submarine springs. An estimated 57-71% of WWTP-sourced N was simulated in the two submarine springs (HHA and HHB) located directly downgradient of the WWTP. This was also validated by the enriched δ 15 N measured at these submarine springs (Fig. 6). Simulated N sourced from OSDS and the eastern boundary were both significantly detected in most wells and submarine springs. Of importance is to note that OSDS-sourced N Fig. 4 Vertical cross section i-i′ (Fig. 2a) of flow vectors for the a density-dependent model and b freshwater-only model. The vector lengths vary according to magnitude. The density-dependent model also contains salinity contours, displaying a relatively thin freshwater lens.
Flow vectors within the lens are directed towards the coast, while flow vectors within the underlying saltwater body are directed inland. As should be expected, freshwater-only vectors move towards the ocean appears in all but one sample across the aquifer, demonstrating the widespread effect of wastewater on groundwater quality. Eight of the 25 simulated samples had at least 50% of the simulated N concentration sourced from OSDS, while 11 samples had at least 50% of N sourced from the eastern boundary (Fig. 6). Conversely, the average OSDS-sourced N concentration was 0.6 mg/L, while the average boundarysourced N concentration was 0.4 mg/L. Calculated δ 15 N values were used to validate simulated N contributions. Results from the described isotope mixing model (section 'δ 15 N model validation') produced a reasonable fit between the measured and computed δ 15 N values (r 2 = 0.76), where 14 of the 25 samples matched within ±1‰ (Fig. 7). Considering that δ 15 N values vary over a range within each source, the mixing results are influenced by uncertainties associated with endmember designations. Since relative N contributions from LULC and WWTP were minimal, variations in δ 15 N LULC and δ 15 N WWTP endmembers did not drastically affect overall δ 15 N model results (Fig. 7a,b). Although the N contribution from the BC was relatively high, the δ 15 N BC range was not large enough to significantly affect δ 15 N model results (Fig. 7c). Since the BC is an ambiguous endmember with mixed water sources,  Corresponding δ 15 N measurements are represented by circle markers. Small grey squares represent OSDS locations and the red square represents the WWTP location the most reliable δ 15 N reference were water samples from highlevel pumping wells, which showed very little variation. Modeling results show that OSDS serves as the most significant N source, thus changes in the δ 15 N OSDS endmember exhibited the greatest variation (Fig. 7d). An enriched δ 15 N OSDS endmember would potentially enrich δ 15 N model results and thus skew relative contributions. There is a large range of OSDS δ 15 N measurements because of the different individual units, but this approach is limited to assigning one OSDS endmember. Therefore, differences between measured and modeled δ 15 N is to be expected; nonetheless, a reasonable fit was obtained with a δ 15 N OSDS value of 7‰. Although this endmember is slightly depleted compared to previously published values from other areas in Hawai'i (e.g., Wiegner et al. 2016), findings from previous studies (e.g., Abaya et al. 2018) already demonstrate that δ 15 N measurements in West Hawai'i are isotopically depleted compared to universal ranges typically found in literature (e.g., Kendall and Aravena 2000). These results therefore emphasize the importance of site-specific measurements, suggesting that more sampling from varying conditions is required and an approach acknowledging wider ranges should be considered. With these limitations in mind, the model captured overall δ 15 N ranges. These results validate the simulated N transport results and show that mixing between N from various sources can be constrained; hence, the model can be used as a tool to estimate relative contributions of multiple nitrogen sources.

Future development and climate change
Future urban development and climate change negatively impacted simulated head, salinity, N and P concentrations to various degrees. When only future development was considered, head was not affected (Fig. 8a) and salinity was minimally affected (Fig. 8b), due to the fact that water fluxes and input salinities were not altered. Under these development conditions, N and P concentrations primarily increased around the WWTP because the WWTP mass load increased to account for the additional wastewater generated with the additional development (Fig. 8c,d). When additional water demand was considered along with the future development, head levels expectedly declined, while salinity, N and P concentrations increased throughout the model due to the reduced water flux applied at the eastern boundary. Further, when climate change impacts were considered along with future development and additional water demand, head, salinity, N and P changes were exacerbated due to the additional water flux reduction. With less freshwater entering the simulated aquifer system, hydraulic head levels declined, with a maximum change of −0.3 m (Fig. 8a). Submarine spring flow also decreased, leading to more saltwater intrusion, thus increasing salinity with a maximum change of 2 ppt. The increased salinities are highest along the southern coast of the aquifer, which is likely due to the narrower aquifer model width (Fig.  8b). Under these future conditions, the increase in N concentrations generally remained under 1 mg/L, but increased notably up to 122 mg/L directly around the WWTP. It appears that reducing the water flux, whether due to increased water demand or climate change impacts, has the greatest impact on simulated results, while nutrient inputs from development green space does not appear to significantly alter simulated nutrient concentrations in the groundwater (Fig. 8c,d). The drastic reduction of water flux ultimately altered nutrient concentrations to a greater degree compared to the addition from urban development.

Modeling sea level rise
With 1 m of simulated sea level rise, the hydraulic head nearly increased by 1 m across the entire area of the modeled aquifer. The unaltered flux at the eastern boundary did not create a significant head gradient with the sea level rise to counter the increase in saltwater head gradient. This also resulted in approximately 1 m of rise of the freshwater-saltwater interface, thus ultimately lifting the entire freshwater lens. The change in salinity, however, did not occur in a simplified trend throughout the aquifer. The maximum change in salinity was 2.75 ppt, but most cells typically demonstrated less than 1 ppt change (Fig. 9). Very little change was observed in the relatively fresher upper layers (Fig. 9a), while the greatest change was observed in the middle of the freshwater lens (Fig. 9c). Additionally, under sea level rise conditions, there was an observed increase in simulated coastal discharge. Diffuse (ocean bottom) coastal discharge decreased by approximately 120,000 m 3 /day, but submarine spring discharge (via drains) increased by approximately 130,000 m 3 /day. Total discharge, therefore, increased by approximately 10,000 m 3 /day with 1 m of simulated sea level rise.
The results of this analysis are supported by Lee et al. (2013), who concluded that when longer time scales are considered, sea level rise could potentially result in an increase in submarine spring flow under certain conditions (Lee et al. 2013). Werner and Simmons (2009) suggest that sitespecific boundary conditions have a significant impact on simulation results when modeling sea level rise in coastal aquifers, and thus argue that groundwater discharge and saltwater intrusion under sea level rise conditions are sensitive to the assigned inland boundary conditions. Werner and Simmons (2009) conclude that with sea level rise, fluxcontrolled systems experience less saltwater intrusion and persistent groundwater discharge in comparison to headcontrolled systems. Results from this study support the concepts proposed by Lee et al. (2013) and Werner and Simmons

Study assumptions and limitations
The model developed for this study simulates groundwater flow and nutrient transport through the Keauhou basal aquifer, calibrated with measured head, salinity, N and P concentrations, and submarine spring fluxes. The major assumption in the approaches adopted concerns validity of the continuum approach (e.g., Bear 1972) for simulating density-dependent flow and transport in a domain that is highly complex with apparent preferential flow. Results suggest that such a concept, as the basis of Darcy's law, is applicable with good calibration accuracy for hydraulic head and submarine spring flows. Some scatter obtained for the results of salinity and nutrient calibration is mainly due to effects of local hydrogeologic medium variability. Use of models that account for small-scale variability such as accounting for preferential flow, is impractical due to the related extensive data needs. Another assumption concerns the treatment of N and P as conservative chemicals based on the justifications presented by Delevaux et al. (2018). Additionally, measurements of NO 3 − + NO 2 − and PO 4 3− were used as major chemical species of N and P, respectively, due to well oxygenated aquifer conditions (Richardson et al. 2017). Considering that a relatively high fraction of TN was represented by NO 3 − in this groundwater, it was assumed that these measurements could serve as representative values for model calibration.
Another model uncertainty is related to the estimation of flow through the eastern boundary of the model. Occurrence of such a flow is confirmed by Fackrell et al. (2020). A reasonable estimate was utilized in this study by a simplified regional model. Since calibration accuracy is sensitive to the eastern boundary condition (section S.3 of the ESM), such an estimate can significantly impact results. There is also an uncertainty regarding the flow condition at such a boundary regarding an assumed vertical flow distribution, considering that the actual hydrogeologic nature of this interflow aquifer boundary is not well understood. In addition, the model did not account for deep offshore freshwater discharge that exists off the coast of Keauhou (Attias et al. 2020).
The model has limitations due to the relatively coarse grid resolution that reduces accuracy of results regarding the spatial distribution of OSDS and LULC. These variables were averaged and assigned as a single lumped value to the corresponding modeling grid since it is certainly impractical to design the grid on an individual OSDS scale. However, numerical accuracy was not impacted by the spatial resolution and mass conservation for water flow and constituent transport was maintained by the model under all applications.
Finally, model calibration is affected by the lack of continuous head records, which is a serious drawback that hinders researchers' efforts in Hawai'i and other similar areas. Available data normally include initial head values measured during the drilling process. Most pumping wells are not continuously monitored, and the initial head measurements were therefore collected at different times. However, the assessment here indicated that the values used for model calibration are within the average fluctuations in water level measurements (section S.1 of the ESM). Fluctuations are not drastic due to the large K h values associated with volcanic rock materials.

Conclusions
As is the case for most coastal aquifers, groundwater from West Hawai'i Island is vulnerable to contaminants derived Positive change represents an increase in salinity with sea level rise, while negative change represents a decrease in salinity with sea level rise from point-and nonpoint-sources such as OSDS, LULC, and WWTP. To ultimately explore local water resource management strategies, this study developed a density-dependent groundwater model for the Keauhou basal aquifer, utilizing multiple measurements and species for calibration. By including head, salinity, N and P concentrations, and submarine spring fluxes in the calibration process, this model obtained overall reasonable result ranges. Comparisons between the results from the density-dependent and a freshwater-only model concluded that the former is a better fit to assess nearshore processes. The freshwater-only model overlooks salinity effects on water circulation, and consequently fails to properly simulate nutrient transport. These effects are critical when investigating submarine springs, which provide important freshwater and nutrient supplies to coastal ecosystems that are susceptible to environmental changes. Calibrating the freshwater-only model necessitates the use of an unrealistic K h due to the failure to account for nearshore flow circulation and the need to define an effective aquifer thickness. As such, in addition to depending on aquifer material-type, K h is a process-based parameter that depends on model conceptualization and the simulated processes. Despite the extensive data and computation expenses, utilization of a density-dependent model is critical for the accurate assessment of water flow and contaminant transport in a coastal aquifer.
Relative N source fluxes were estimated from the calibrated model and validated by using δ 15 N of NO 3 − , which is the dominant chemical species in this oxic aquifer. Simulated N was primarily sourced from OSDS and the eastern boundary. Further, N in submarine springs downgradient of the Kealakehe WWTP were overwhelmingly sourced from the WWTP, as validated by δ 15 N results. Hence, there is a need to improve on current practices for collecting and treating wastewater and releasing effluents from a single, seemingly inefficient WWTP.
As should be expected, future developments and climate change impact head, salinity, and nutrient concentrations. Specifically, the reduced water flux caused by higher water demands and reduced recharge induced by climate change resulted in declining hydraulic head levels and greater saltwater intrusion. Additionally, contaminant concentrations increased, most prominently around the WWTP. A 1 m of sea level rise increased hydraulic head by an equivalent amount aquifer wide, but effects on salinity distribution prove to be more complicated. Sea level rise can effectively increase submarine spring flux due to an increase in head gradient.
To further improve the model, a need exists for a better understanding of the eastern boundary structure and assigning its related aquifer interflow. There is also a critical need to identify the range of OSDS δ 15 N, which has proven to vary in isotopic composition. These details have a significant impact on the modeled results, which if used to inform groundwater management without acknowledgement of the various uncertainties, could result in unexpected or undesirable outcomes. Despite the assumptions and limitations, however, it is believed that this model is qualified for use in understanding various water flow and transport processes in the Keauhou coastal aquifer and as a management tool for assessing relative aquifer responses under various management scenarios. In addition to the examples discussed in the study, the model can address others, including OSDS upgrade efforts that include LULC and climate change, as well as the adverse effects in the nearshore environment in response to changes in both quantity and quality of submarine springs. Although this model was developed for a specific aquifer, the methodology and analyses are applicable to various coastal aquifers, where water quantity and quality are of major concern, especially with urban expansion and the increasing severity of climate change-related problems. statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.