Predicting Coral Reef Carbonate Chemistry Through Statistical Modeling: Constraining Nearshore Residence Time Around Guam

To accurately predict the impacts of ocean acidification on shallow-water ecosystems, we must account for the biogeochemical impact of local benthic communities, as well as the connectivity between offshore and onshore water masses. Estimation of residence time can help quantify this connectivity and determine the degree to which the benthos can influence the chemistry of the overlying water column. We present estimates of nearshore residence time for Guam and utilize these estimates to model the effects of benthic ecosystem metabolism on the coral reef carbonate system. Control volume and particle tracking approaches were used to estimate nearshore residence time. These estimates were paired with observed patterns in the reef carbonate system around Guam using water samples collected by NOAA’s National Coral Reef Monitoring Program. Model performance results suggest that when considering the effects of benthic metabolism on the carbonate system, it is paramount to represent the contact time of the water volume with the benthos. Even coarse estimates of residence time significantly increase model skill. We observed the highest predictive skill in models including control volume derived estimates of residence time, but only when those estimates were included as an interaction with benthic composition. This work shows that not only is residence time critically important to better predict biogeochemical variability in coral reef environments, but that even coarse hydrodynamic models can provide useful residence time estimates at management relevant, whole-ecosystem scales.


Introduction
Ocean acidification (OA) has already reduced global average surface ocean pH by at least 0.1 since the beginning of the industrial era (Sabine et al. 2004;Jiang et al. 2019), and models forecast an additional decline of 0.3-0.4 pH units by the end of the century (Brewer 1997;Wickett 2003, 2005;Jiang et al. 2019;Orr et al. 2005). Coral reef ecosystems are particularly vulnerable to OA (Hoegh-Guldberg 2007;Andersson and Gledhill 2013). Scleractinian corals, which form new reef framework habitat through the process of calcification, experience impaired skeletal growth under low pH conditions (Andersson et al. 2009). Simultaneously, acidification causes carbonate dissolution and higher bioerosion (Wisshack et al. 2012;Webb et al. 2017;Stoltenberg et al. 2021), tipping the balance toward net habitat erosion (Enochs et al. 2015a(Enochs et al. , 2016 and contributing to ecosystem shifts (Enochs et al. 2015b). Understanding the relative importance of different drivers of biogeochemical variability on coral reefs-especially characteristics or functions that can be restored, protected, or amplified-may provide a means to locally manage the global challenge posed by OA.
Global OA projections typically focus on changes to the carbonate system in openocean surface seawater. Open ocean models can widely predict oceanic carbonate chemistry (Ricke et al. 2013;Jiang et al. 2019), especially where ample data is available (Lauvset et al. 2016 (GLODAP v2); Bakker et al. 2016 (SOCAT v3)). The complexity of coastal systems, however, has made it difficult to make similar OA predictions in shallow-water coral reef environments. The vulnerability of coral reefs to OA is complicated by natural diel and seasonal cycles in the overlying seawater carbonate chemistry, which can structure the physiological responses of keystone species (Enochs et al. 2018) and modify their exposure . Primary production, including photosynthesis from corals, algae, and seagrass beds, can buffer acidification stress when daytime photosynthesis elevates local pH as dissolved inorganic carbon (DIC) is removed (Unsworth et al. 2012;Manzello et al. 2012). Respiration can exacerbate the effects of global OA by driving an influx in DIC and contributing to a further decline in pH. Calcification and dissolution also influence DIC, as well as total alkalinity (TA; Smith and Key 1975). The diurnal metabolic cycle, combined with long-term OA-driven decreases in buffering capacity, may lead to an amplification in the variability experienced on reefs, exposing the coral reef community to exacerbated low pH conditions as well as opportunities for buffered conditions (Shaw et al. 2012).
Nearshore environments are spatially and temporally dynamic with respect to carbonate chemistry due to a combination of biological and physical processes. While alterations to seawater carbonate chemistry are driven by benthic community metabolism, they are further influenced by environmental factors including light intensity, nutrients, flow speed, and depth (Langdon and Atkinson 2005;Page et al. 2016Page et al. , 2019Monismith et al. 2007;Falter et al. 2013;Cyronak et al. 2020). Flow speed and depth are especially important because they determine the residence time and the volume of water that interacts with the benthos, therefore altering the magnitude of the metabolic signal. A longer residence time provides greater contact time between the benthos driving the metabolic change and the seawater, resulting in greater chemical modification. Residence time should be properly accounted for as a driver of carbonate chemistry variability-especially in the context of the benthic community structure on the reef.
Understanding the local drivers of biogeochemical variability on coral reefs is critical for effective management of OA at local scales. Widespread acidification threatens to reduce calcification and shift coral reefs to a state of net dissolution by the end of the century (Silverman et al. 2009), but experimental (Albright et al. 2016) and modeling (Mongin and Baird 2016) studies have demonstrated that a reversal of this trend is possible at the local scale through interventions to buffer OA. Albright et al. (2016) observed increased net community calcification in a natural coral reef community in response to enriched alkalinity (accomplished by pumping a sodium hydroxide [NaOH] solution onto the reef flat). And Mongin and Baird (2016) demonstrated the feasibility of reef-scale OA mitigation by using seaweed farming to remove carbon and locally increase aragonite saturation state. These studies show how local mitigation and reversal of OA trends may be possible at the reef-scale.
Some nearshore environments already act as natural coral reef refugia due to local processes (Shamberger et al. 2014;Barkley et al. 2015). To manage reef persistence into the future, we must better understand and predict the factors that contribute to ecosystem vulnerability and/or resilience to climate change stressors. By understanding natural drivers of resilience, managers can take actions to enhance the features that contribute to resilience. For instance, planning coral restoration projects in proximity to communities of primary producers that can naturally buffer daytime pH (Manzello et al. 2012;Mongin and Baird 2016), slowing the reduction in aragonite saturation, and supporting calcification.
To predict spatiotemporal gradients in carbonate chemistry on coral reefs, models need to incorporate data about the benthic community driving metabolic changes in the carbonate system, as well as the local hydrodynamics modulating the magnitude of those changes. While the former is available at multiple resolutions/scales (e.g., diver surveys-satellite remote sensing), there is a distinct paucity of applicable flow data. Instrument-based measurements, for example, typically represent small footprints of the reef environment, and cannot be used for evaluating the flux of large water masses over multiple habitats. Readily available, large-scale hydrodynamic models, however, usually describe basin-scale circulation, and rarely capture the intricacies of nearshore environments. While a few detailed hydrodynamic models have been built for specific coral reefs (e.g., Palmyra (Rogers et al. 2016;Rogers et al. 2017) and Kaneohe Bay (Lowe et al. 2009)), these models are limited in their scope and may represent specific time domains.
Here we develop a modeling approach that incorporates both benthic data and residence time estimates at large, whole-ecosystem scales to better predict nearshore changes in the carbonate system. We apply this approach to the nearshore waters of Guam, an area with diverse benthic communities and complex flow regimes. We evaluate two residence time estimation techniques (control volume and particle tracking) to assess the importance of including residence time in the model, as well as other drivers including benthic composition, hour of the day, and light availability. An improved capacity to model reef metabolism-driven changes in the carbonate system will enable us to better locally manage coral reefs by predicting the features that contribute to vulnerability and resilience.

Methods
This model requires four main sources of data: carbonate chemistry, benthic composition, sunlight exposure, and estimates of residence time from local hydrodynamics information. Below we describe each of these components in detail, ending with our assessment of model performance.

General Site Description
Guam is the southernmost island in the Mariana Archipelago in the western Pacific and the largest island in Micronesia. The island is surrounded by a diversity of reef types: fringing reefs border most nearshore areas, barrier reefs surround the patch reef and lagoonal areas of the southern Cocos Lagoon and Apra Harbor, and the western and eastern coasts are dotted with several sheltered bays (including Tumon Bay, Agana Bay, and Pago Bay). There are distinct biological and hydrodynamic environments on the windward (east) and leeward (west) sides of the island driven by easterly trade winds. The predominant benthic cover on Guam's reefs is turf algae (59-65%) with relatively low coral cover (the highest mean island coral cover was 14.3% in 2011 before decreasing to 13.0% in 2014 and 11.6% in 2017) and a small percent cover of crustose calcifying algae (CCA, 4.1-6.8%; see Barkley et al. 2022).

Carbonate Chemistry
A total of 64 discrete seawater samples (  (Dickson et al. 2007). Fifty discrete seawater samples were collected around Guam during daylight hours (08:00-15:00 local time) at depths ranging from 0.9 to 16.5 m. Fourteen additional samples were collected at a single site on the south coast to capture natural diurnal fluctuations in carbonate chemistry (every 4 h for 24 h on April 6-7, 2017). Additionally, three offshore surface samples (over 10 km from land) were collected in 2017. Temperature, salinity, and pressure were measured concurrently with a Seabird 19 + conductivity-temperature-depth (CTD) profiler. Seawater samples were shipped to the NOAA Pacific Marine Environmental Laboratory (PMEL) in Seattle, WA and analyzed for TA (open-cell titration using an instrument custom-built by the Dickson Lab, Scripps Institution of Oceanography) and DIC (coulometric titration using two Single Operator Multiparameter Metabolic Analyzer [SOMMA] systems). TA and DIC were used in conjunction with temperature, salinity, and pressure to calculate the full carbonate system parameters in the R package seacarb (Gattuso et al. 2020), using dissociation constants from Lueker et al. (2000). Additional information about carbonate chemistry data collection and analysis for the Marianas region is available in Oliver et al. (2018), Barkley et al. (2017), and Barkley et al. (2021).

Determining Offshore-Onshore Changes in Carbonate Chemistry
Nearshore changes in DIC and TA (values normalized to a salinity of 35) were calculated using the discrete measurements on the reef minus an offshore endmember value. The offshore endmember DIC and TA was the average of the three discrete sample values collected. Since all three offshore samples were collected in 2017 but onshore samples span 2011-2017, we corrected all offshore reference values for global progression of DIC using the linear temporal trend in DIC (Supplementary Material, Figure S1) from the Hawaii Ocean Time-series data (https:// hahana. soest. hawaii. edu/ hot/ metho ds/ dicalk. html) and the difference in time from collection of the onshore samples and its paired offshore reference.

Benthic Composition
Benthic cover and community composition data were collected using a stratified random benthic survey design according to standard NCRMP sampling protocols (Ayotte et al. 2015;Winston et al. 2020;NOAA Coral Program 2021). In total, 346 unique sites were surveyed around Guam in 2011, 2014, and 2017: May 5-9 and June 6-16, 2011 (132 sites), March 25-April 4, 2014 (117 sites), and May 3-14, 2017 (103 sites). Percent benthic cover was determined from photoquadrats using Coral Point Count with Excel extensions (CPCe; Kohler and Gill 2006) in 2011 and 2014 and using CoralNet (Beijbom et al. 2015) in 2017. The organism (genus/morphology for corals, genus/functional group for algae) or type of substrate was identified for ten randomly overlaid points in each image (Lozada-Misa et al. 2017). For the purposes of our model, we simplify benthic composition as the percent cover of calcifiers (scleractinian coral; Halimeda, crustose coralline algae and Peyssonnelia) and non-calcifying algae (macroalgae and turf).

Sunlight Exposure
We paired 8-day composite Photosynthetically Active Radiation (PAR) data from Moderate Resolution Imaging Spectrometer/Aqua level 3 binned imagery with our sampling dates of interest, matching the sampling date with the composite's centered time. This data 1 3 can be downloaded from the ERDDAP server (https:// upwell. pfeg. noaa. gov/ erddap/ gridd ap/ erdMH 1par0 8day. graph).

Nearshore Sector Designation
Distinct nearshore environments around Guam were determined using a statistical downscaling approach developed by Oliver et al. (2020) based on spatially contiguous clustering and mixed model analysis. Sector boundaries were informed by clustering of carbonate chemistry samples (DIC and TA) and were achieved by linking neighboring sampling points, estimating the similarity between linked points, and finally calculating a minimum spanning tree to define natural divisions between groups of neighbors (Assunção et al. 2006). The resulting groups underwent mixed model analysis (Zuur et al. 2009) to define the optimal balance between fine spatial resolution and statistical robustness, thereby informing the total number of distinct sectors around the island (Fig. 1b). Sectors were translated into the 2 km gridded environment of the hydrodynamic model data (Fig. 1c) and sector IDs were assigned to each grid cell based on the sector number associated with it as well as its identity as a "land" cell or an "ocean" cell with velocity data. Cells on the border between sectors were assigned the numbers of both sectors. The resulting sector size is a function of three drivers: providing enough samples per sector to be statistically robust, capturing sample clusters that have low within-group variation in their measured carbonate systems, and resulting in a sector size large enough to accommodate multiple 2 km grid cells.

Estimation of Residence Time
Nearshore hydrodynamic information within predefined sectors ( Fig. 1b) was used to estimate residence time from both a control volume (CV) approach using volume fluxes calculated over the boundaries of each sector and a particle tracking (PT) approach that traced the path of water parcels through each sector. Estimates of residence time were made using Regional Ocean Modeling System (ROMS) model 3-hourly flows for the area of interest over model years 2015-2021. Control volume residence times were estimated from March 17 to April 11 (2016-2021) and from April 25 to May 21 (2015-2021). Particle tracking residence times were estimated over the same time periods, but only for 2018-2021 due to the inconsistency in landmass extent prior to 2018 (see Supplementary Material for details about land mask, Figure S2).

Hydrodynamic Information Availability
ROMS is a free-surface, terrain-following, primitive equations ocean model (http:// myroms. org). We used a regional implementation of ROMS that covers an area around Guam and parts of the CNMI (Powell 2013 April 11 (2016April 11 ( -2021 and from April 15 to May 21 (2015-2021), encompassing the days that carbonate chemistry was sampled. While a 2 km grid is still coarse for studying nearshore dynamics, it is the highest resolution hydrodynamics data available for this region.

Control Volume Approach to Residence Time
The 3 h volume flux into and out of each sector was calculated using the ROMS velocity field. Both in and out fluxes were estimated as the sum of 3-hourly current flows through the grid cell area over the boundaries between sectors as well as over the sector-ocean boundary of each sector to the mixed layer depth (MLD). In areas where the MLD differed between two adjacent points composing a cell boundary the deeper depth was used. The control volume approach was not used to estimate residence time if data were missing in a sector. For instance, on most dates prior to 2018, in sectors 1 and 3, the land mask was larger for Guam in the ROMS model (see Supplementary Material Figure S2 for land mask discrepancy). This resulted in missing data for the nearshore cells at the northern (sector 1) and southern (sector 3) ends of the island, so residence time was not estimated in those sectors.
Boundary fluxes were converted to residence time (hours) by dividing by sector volumes, which were calculated by multiplying the surface area of each cell in a given sector by the MLD or seafloor depth (whichever was shallower), and sea surface height of that cell. This guarantees that we are focusing our calculations on a well-mixed, upper water mass. Finally, a hindcast running mean (12 and 24 h) was calculated for the reported residence time to capture the relevant benthic effect on water mass inside each sector. From this point, the two control volume (CV) residence time estimations will be described as RT CV12 and RT CV24 , to distinguish between the running means applied.

Particle Tracking Approach to Residence Time
Residence time was also estimated using a Lagrangian particle tracking approach. Particle trajectories were modeled using the Parcels v.  Figure S3) and released every two hours for the duration of the period specified (March 17-April 11 [2018][2019][2020][2021] and April 15-May 21 [2018][2019][2020][2021]). Particles were advected using the ROMS surface velocity field with eddy diffusivity added (10 m 2 s −1 , adopted from Lindo-Atichati et al. 2020) using a five-minute model timestep. After the first round of advection and diffusion, each particle was evaluated to determine if it had been "beached" on a land cell. Beached particles were returned to their previous location and re-advected using an adjusted longshore velocity field (further described in Supplementary Material, Text S2 and Fig. S4) to simulate longshore flow and prevent advection onto the island. Following this process (or if the particle was never beached), a new round of advection and diffusion was conducted to continue the particle trajectory.

3
To avoid the velocity components going to zero at the coast, the default no-slip boundary condition in Parcels was replaced with the free-slip interpolation condition that separately considers the cross-shore and along-shore velocity components and corrects the velocity field based on the direction of the land boundary. Model output of particle locations was saved at one hour intervals.
Residence time was determined for two domains using the modeled trajectory output from Parcels. For every particle inside a nearshore sector at each time point, we determined the total uninterrupted time the particle had spent in that focal sector so far. We also calculated the uninterrupted time spent in any nearshore sector to capture the time when the benthic environment of another sector could influence the water before it arrived in its current sector. We followed the age tracer approach described in Mongin and Baird (2014; see also Monsen et al. 2002 andHall andHaine 2002) to account for the accumulation of nearshore time as each particle was advected and diffused. This approach also incorporates the expected decay in metabolic signal outside of the nearshore sectors when water column processes and air-sea gas exchange likely counteracted reef-driven changes in carbonate chemistry. While we based our expected decay on the treatment described in Mongin and Baird (2014), we applied a more aggressive damping rate of 0.8 h −1 for each hour spent outside a nearshore sector. From this point on, the two particle tracking (PT) estimates of residence time (RT) will be referred to according to the portion of their trajectory considered, focal sector only or all nearshore sectors, as RT PTsector and RT PTnear, respectively.

Climatological Residence Time
Because of the availability of ROMS data, our residence time calculations start on April 15, 2015. Thus, to pair estimates of sector residence time with all 64 carbonate chemistry samples collected around Guam in 2011, 2014 and 2017, we calculated three-hourly climatological residence times and paired them with the carbonate chemistry samples by fractional Julian day (i.e., day 1-366 of a given year, with fractional 3 h intervals). We generated the annual 3-hourly climatology by summarizing the residence time by sector for each fractional Julian day. For the control volume estimates (CV), we generated a climatological mean residence time for two periods: April 15-May 21, spanning the years 2015-2021, and the earlier period, March 17-April 11, for which data were only available in 2016-2021. For the particle tracking estimates (PT), we spanned the years 2018-2021 for the same two periods. Sectors varied in their temporal coverage, with sectors 2, 4, 5, and 6 having all seven years of data present, but due to changes in the ROMS model's land mask (Supplementary Material, Figure S2), sectors 1 and 3 lacked estimates for 2015 and 2016, and had only partial coverage in 2017.
The fractional Julian day was also calculated for each carbonate chemistry sample so that the relevant sector residence time could be matched to each discrete sample for modeling. Thus, the carbonate chemistry samples collected in 2011, 2014, and 2017 were matched with the appropriate Julian day climatological residence time.

Modeling Nearshore Carbonate Chemistry Deltas
We built linear models of increased complexity using the lme4 package for R (Bates et al. 2015) to describe the nearshore changes in DIC and TA (the difference between endmember chemistry and chemistry on the reef) in 64 discrete samples for Guam. The null model included only benthic community structure (the percent cover of calcifiers and percent cover of non-calcifying algae). We then successively added parameters to build more complex models, starting with hour of the day, then satellite-derived photosynthetically available radiation (PAR), before including estimates of climatological residence time using each of the four methods described above (RT CV12 , RT CV24 , RT PTsector , RT PTnear ). Finally, we included each climatological residence time estimate as an interaction with the benthic composition.

Model Performance
The predictive skill of each model was evaluated by comparing adjusted R 2 values and the difference in Akaike's Information Criterion (ΔAIC). Adjusted R 2 values allowed comparison of model accuracy, or the goodness-of-fit for each linear model. ΔAIC provides a way to compare the skill and simplicity of each model, allowing us to evaluate if increasing complexity in the model improves the fit without overfitting. An increase in model skill is indicated by a higher adjusted R squared value and a lower ΔAIC.

Results
Our models utilized 64 discrete carbonate chemistry samples collected around Guam in 2011, 2014, and 2017 as well as benthic survey data to describe the communities present at the sites where those samples were collected. These two foundational data sources are summarized in Fig. 2.   Fig. 2 a Benthic composition by sector. The x axis is the percent cover of non-calcifying algae (macroalgae, other encrusting macroalgae, and turf) and the y axis is the percent cover of calcifiers (corals, crustose coralline algae (CCA), Halimeda algae, and calcifying Peyssonnelia sp. algae). Points represent the mean benthic composition for each sector in each sampling year. Error bars depict standard error; b Nearshore DIC (µmol kg −1 ) and TA (µmol kg −1 ) by sector. Triangles represent offshore endmember chemistry. Lines depict the expected changes in chemistry driven by photosynthesis/respiration (horizontal left/right) and calcification/dissolution (2:1 slope up/down). All points are colored by sector number and the inset map (top right) shows the sector locations around Guam. Point shapes indicate the sample year 1 3

Benthic Composition
The highest percent cover of non-calcifying algae was observed at the southern end of Guam, in sectors 3-5 (Fig. 2a). Sector 6, on the western side of Guam, had the lowest noncalcifying algal cover and the highest percent calcifier cover.

Nearshore Carbonate Chemistry
The largest changes in both onshore TA and DIC relative to the offshore endmember (TA: 2313.6 µmol kg −1 , DIC: 1980.0 µmol kg −1 ) were observed in sector 6 on the west side of Guam, where TA was as low as 2244.6 µmol kg −1 and DIC was as low as 1874.9 µmol kg −1 (Fig. 2b). Almost all nearshore sites had lower TA relative to the offshore endmember sample. Typically, DIC was also lower nearshore, with the exception of sectors 1 and 3, where DIC was as high as 1994.9 µmol kg −1 and 1991.6 µmol kg −1, respectively. Sector 3 includes data from the diel suite experiment in 2017 and therefore, includes samples collected late in the day and night, when we expect respiration to elevate DIC in the absence of photosynthesis.

Estimated Residence Time by Sector
Among the four residence time estimates compared for each sector (Fig. 3), the control volume estimates for sector 6 deviated the most from sector 6's particle tracking estimates and from residence estimates from any methods for the other sectors. In sector 6, RT CV12 and RT CV24 showed median values of 98.28 and 90.37 h, respectively, while the next highest median residence times were observed in sector 2 and 3's particle tracking estimates, specifically PT near (Fig. 3). The particle tracking estimates for the southern coast and leeward (west) coast of Guam (sectors 3-6) were relatively low Fig. 3 Boxplots of residence time estimates (hours) by sector. Each sector includes residence time estimated using each of the four methods. Control volume (RT CV12 and RT CV24 ) estimates are shown in blue colors and particle tracking (RT PTsector and RT PTnear ) estimates are shown in green colors. Y axis is a log scale (minimum of 8.5 h in sector 5), and not significantly different between trajectory summing methods (nearshore versus focal sector) with the notable exception of sector 3, where the RT PTnear was significantly higher (ranging from 45.3 to 74.2 h). RT PTnear was similarly higher relative to RT PTsector in sector 2, on the windward east coast. In sectors 1 and 2 (wrapping around the northern coast and down the east), both particle tracking estimates yielded relatively higher residence times compared to the control volume estimates. This was also true in Sector 3 (southern coast) for the nearshore particle tracking estimate, but the focal sector particle tracking estimate was not significantly different from the control volume estimates.

Model Skill
We evaluated the performance of models including benthic composition, hour of the day, satellite-derived PAR, four residence time estimates, and residence time as an interaction with benthic composition (Fig. 4). Increasing model complexity beyond basic benthic composition, by including hour of the day and then PAR, did not substantially increase model skill. Inclusion of residence time in the models usually increased model skill, though not in all cases. The greatest increase in model skill resulted from including the control volume residence time estimates (RT CV12 and RT CV24 ) as an interaction with the benthic composition. The best model skill in predicting ΔDIC was observed using RT CV12 (Fig. 4a), whereas the best skill for predicting ΔTA was seen using RT CV24 (Fig. 4b). Including control volume estimates without the benthic interaction did not improve model skill relative to models lacking residence time.
Models including the residence time estimates derived from nearshore particle tracking (RT PTnear ) also outperformed models lacking residence time, with or without the Fig. 4 Model performance comparison for models with and without residence time included. Models without residence time are plotted in black. Models including residence time (not as an interaction) are unfilled with shape and color indicating the residence time estimate used. Filled markers indicate models that include residence time as an interaction with benthic composition. The legend indicates with parameters were included in each model (Benthos = benthic composition, HR = hour of the day, PAR = satellite-derived Photosynthetically Active Radiation, RT_CV12 = control volume (12 h running mean) residence time estimate, RT_CV24 = control volume (24 h running mean) residence time estimate, RT_PTsector = particle tracking (focal sector) residence time estimate, RT_PTnear = particle tracking (nearshore) residence time estimate, RT:Benthos = residence time as an interaction with benthic composition) benthic interaction included (Fig. 4). Among the models including RT PTnear , the inclusion of residence time without the interaction with benthic composition outperformed the models that included RT PTnear as an interaction. The models including RT PTsector did not demonstrate a notable increase in model skill.

Overview of Model Performance
Coral reef metabolic processes can counteract or exacerbate the local effects of OA (Cyronak et al. 2018). Therefore, improving quantitative predictions of how these processes affect the local carbonate system may help us better manage reefs in the face of climate change. Small-scale and experimental studies have recognized residence time as an important driver of local biogeochemical variability (Page et al. 2019;Cyronak et al. 2020), since the degree of metabolic alteration of the carbonate system is a function of the time that a given parcel of water spends in contact with the local biological community. However, residence time estimates are often unavailable at ecosystem scales relevant to management, and therefore there is often a disconnect between nearshore carbonate system modeling and direct management applications.
Here we demonstrate the utility of even coarse estimations of residence time in predicting the contribution of benthic metabolism to changes in the nearshore carbonate system on the reefs around Guam. Our statistical modeling results confirm that, while models that lack residence time as a correlate have effectively no power to explain reef carbonate chemistry dynamics, the inclusion of residence time elevated predictive skill, even for very simple estimation methods, i.e., those based solely on the current flux over nearshore sector boundaries (CV approach- Fig. 4). Indeed, using residence time calculated with the control volume approach yielded more predictive power than the particle tracking approach, which was only effective when considering residence time in all nearshore sectors along the particle trajectory.

Strengths of Our Approach in the Context of Other Methods used to Understand Carbonate Chemistry Variability on Coral Reefs
The modeling method presented here is a valuable addition to existing approaches commonly used to describe and quantify benthic metabolic interactions with the seawater carbonate system. The method (a) estimates hydrodynamics at a scale relevant to reef interactions, (b) provides repeatable, high-coverage sector-scale estimates of benthic metabolic effects, and (c) incorporates both benthic cover and exposure time, thus disaggregating physiology and local hydrodynamics. Of the common methods utilized, no single approach meets all these functions. Most previous descriptions of reef biogeochemical variability describe simple flow environments or capture limited footprints of benthic impact. Methods widely used to describe and quantify metabolic processes on reefs (e.g., net ecosystem calcification, TA-DIC regressions, etc.) are limited in the ability to capture both the direction and magnitude of metabolic changes, while also accurately attributing these changes to the relevant combination of drivers.

3
The importance of residence time in our models corroborates previous studies that have examined drivers of change in the carbonate system on coral reefs and have highlighted time spent on the reef as an important constraint to metabolic change (Falter et al. 2013;Cyronak et al. 2020). While the influence of residence time on regulating metabolismdriven changes in water column chemistry is well established in the literature with previous studies pointing to the importance of depth and currents as drivers of coral reef biogeochemical variability (Falter et al. 2013;Page et al. 2019;Kekuewa et al. 2021 among many others), few studies have described detailed local hydrodynamics outside of very simple flow environments (slack water method for negligible flow [e.g., Ohde and van Woesik 1999;Shaw et al. 2012;Kekuewa et al. 2021] or unidirectional flow [e.g., Koweek et al. 2015a, b;Albright et al. 2013;Falter et al. 2012;DeCarlo et al. 2017]) where general behavior and the effects of currents and/or residence time can be better assumed. Predicting residence time in complex coral reef environments is critical for appropriately describing the physical mechanisms that limit or amplify metabolic effects on carbonate chemistry. Some reef hydrodynamics have been modeled in detail (e.g., Palmyra Atoll, Rogers et al. 2016;Rogers et al. 2017), but hydrodynamic data and models are lacking for most shallow coral reef ecosystems. Despite limitations in spatial resolution, our estimates of residence time on the nearshore coral reefs of Guam represent an important step toward accounting for appropriate timescales for metabolic change in the carbonate system on coral reefs and predicting future carbonate chemistry variability at spatial and temporal scales that are relevant for reef management.
Previous studies have used net ecosystem calcification (NEC) to describe the balance between calcification and dissolution in coral reef environments (see DeCarlo et al. 2017 for a recent synthesis). NEC is predicted to decline with open ocean OA, ultimately resulting in a shift on reefs from net accretion to net erosion (Silverman et al. 2009). It is typically calculated using the alkalinity anomaly technique and applied over small-scale reef environments. The sign change of the alkalinity anomaly conveys whether a reef is undergoing net dissolution or net calcification; however, determining the exact rate of change in CaCO 3 driven by precipitation and dissolution requires knowledge about flow conditions (residence time) and volume. Usually, residence time is determined by tracking a parcel of water across the reef using a Lagrangian (e.g., Albright et al. 2013;Koweek et al. 2015a, b) or Eulerian (e.g., Falter et al. 2012;DeCarlo et al. 2017) framework. In other cases, NEC is determined under conditions of negligible flow (slack water method, e.g., Ohde and van Woesik 1999;Shaw et al. 2012;Kekuewa et al. 2021). Less conventional alkalinity anomaly techniques have been applied to slightly more complex environments; for instance, Shamberger et al. (2018) estimated NEC in a semi-enclosed lagoon in Palau using a combination of TA, salinity, and volume budgets.
As the slopes of TA-DIC plots change across a reef in space and time, they reflect the primary drivers of metabolic alteration, i.e., NEC and net ecosystem production (NEP), a balance between photosynthesis and respiration (see synthesis in Cyronak et al. 2018). In communities dominated by NEP, the slope will approach 0 and in communities dominated by NCC, the slope will approach 2. Changes in TA and DIC, as visualized in propertyproperty plots, can indicate the sign and magnitude of change in chemistry driven by metabolic processes especially when combined with data describing proportions of calcifying and non-calcifying organisms in the area. For instance, in this study, reductions in TA and DIC relative to offshore indicate a combination of photosynthesis and calcification across all sectors (Fig. 2) around Guam. While property-property plots and TA-DIC slopes can be useful tools for understanding changes in the reef carbonate system over space or time, including breaking down the dominant processes driving seawater chemistry modification 1 3 in different reef habitats (Andersson and Gledhill 2013;Anthony et al. 2011), these data do not provide a measure of the time scale over which changes occur. It is difficult to determine whether the magnitude of change should be attributed to changes in metabolic rate or if it is the result of exposure time to the benthic community.
While we can use TA-DIC slopes and property-property plots to infer the metabolic processes driving change in the carbonate system, we cannot directly assess the metabolic rates responsible for changes in chemistry. Even if the changes are likely linked to benthic habitats, we do not know how long water is spending in each environment for the signal to manifest. A notable strength of our approach here is that our residence time estimates provide a timescale for local change that we can pair with the sign changes in TA and DIC as well as with information about benthic composition. Together, this information provides us with a more complete description of the local drivers of change. Furthermore, we are building a modeling framework where predictions of nearshore reef carbonate chemistry can be made in areas without direct in situ measurements of the carbonate system (a requirement for both the NEC and TA-DIC regression methods of quantifying coral reef metabolism).
Our approach incorporates hydrodynamics at a scale relevant to reef interactions with the water column, but we also generalize by dividing Guam's nearshore environments into only six distinct sectors. This limitation is imposed by the spatial and temporal resolution of the hydrodynamic model data readily available as well as the spatial distribution of the discrete seawater samples. The benthic environment is variable inside each sector, and residence time varies over space and time within these sectors, too. Nonetheless, capturing even a fraction of the physical variability at this scale is valuable as it allows us to assess the independent and interactive impacts of biological and physical drivers of carbonate system variability.

Relative Performance of Control Volume Versus Particle Tracking Estimates
The interaction of residence time and benthic composition captures the mechanisms behind benthic community metabolism driving changes in chemistry (ΔDIC and ΔTA) and residence time determining the magnitude of that change. Inclusion of the control volumederived residence time estimates (RT CV12 and RT CV24 ) resulted in the greatest model skill when those residence times were included as interactions with benthic composition, but model skill was low when CV residence times were included independently. This suggests that this relatively simple estimate of the time a water mass spent in each sector captures the relevant timescale for the benthic environment to interact with the water column, and therefore has predictive utility. The more computationally complex particle tracking estimates of residence time also improved model skill, but the importance of the interaction with benthic composition was less clear and the two distinct particle tracking estimates varied considerably in their utility, with RT PTnear estimates outperforming the focal-sector-focused RT PTsector . The RT PTsector interaction model was among the worst for predicting both ΔDIC and ΔTA. This may be the result of RT PTsector not fully capturing the time that water spent in contact with a reef environment since this approach was limited to the time spent in only the focal sector. If water was previously exposed to a metabolic influence, that exposure, and any modification, would be underrepresented by this estimation.
The nearshore method (RT PTnear ) exhibited increased model skill both with and without an interaction with the benthic composition; however, the model lacking the interaction term performed better, unlike the pattern seen with CV residence time estimates. These results suggest that RT PTnear represents the time that water has spent in contact with the reef better than RT PTsector . The poor performance of models including the benthic interaction may be caused by RT PTnear including residence time on reefs outside the focal sector yet still pairing those estimates with benthic cover data from only the focal sector. Therefore, for RT PTnear , all change in chemistry is credited to the focal sector's composition even if a different benthic environment is responsible for part of the change.

Complexity of Nearshore Flow Limits Utility of Particle Tracking Estimates
The relatively poor performance of particle tracking estimates of residence time may be due to their reliance on interpolated velocity fields within the nearshore boundary conditions in our relatively coarse model descriptions of these complex coastal hydrodynamic environments. Despite providing the best available spatial and temporal resolution for the region, the 2 km ROMS model for Guam does not resolve nearshore flow. This limitation is particularly relevant for a method that is more reliant on detailed flow, like particle tracking, since ROMS does not capture the hydrodynamic complexity of the reef that we expect to dictate flow trajectories (for instance, wave action over shallow reefs). Interpolation of current velocity across the 2 km grid cells likely resulted in the particle tracking residence time estimates overestimating the time spent in each sector on the windward, upstream side of Guam, while underestimating the time spent on the leeward, downstream side of Guam (Fig. 1a).
In Parcels, the Python framework used to simulate Lagrangian particle tracking, particles are advected then diffused in each timestep (Supplementary Material, Figure S4). Particles slow down as they approach the zero velocity of "land". Depending on the timestep, this could bring a particle very close to the shoreline, but that would mean the next interpolated velocity used for advection is extremely small. We used the free-slip interpolation method to minimize this behavior by redirecting the velocity vector in the direction perpendicular to the land cell (not allowing velocity to go all the way to zero). We also utilized an "unbeaching" process described in the methods. Despite these interventions, particles that were advected toward Guam from the east followed slower trajectories down the coast before being pushed offshore (usually around the southern end of the island from east to west). Longer residence times were especially evident in RT PTnear since this method summed the time a trajectory spent in any nearshore sector, not just the focal sector. Thus, many particles in sector 3 (south) had previously followed trajectories along the east coast through sectors 1 and 2 resulting in very long residence times for RT PTnear , whereas RT PTsector only used a fraction of that path time.
In contrast, on the west side of the island, particles seeded further from Guam were less likely to reach the western nearshore sectors as the offshore velocity field swept them westward, away from the island. The particles that did travel inside the sectors on the west coast of Guam (sectors 4, 5, and 6) had shorter residence times relative to their eastern counterparts. They also had shorter residence times relative to the control volume estimates. This was especially evident in sector 6 where the control volume-derived residence times were especially long due to relatively low fluxes across the sector boundaries. Another important distinction between our CV and PT methods-and a possible explanation for some of the estimate discrepancies-is that the PT analysis was restricted to surface flow, while the CV integrated fluxes from the surface down to the mixed layer. Vertical velocities are not available in the ROMS data for the Marianas region, limiting our ability to model particle trajectories in three dimensions.
Together, these four estimates of residence time demonstrate the hydrodynamic variability along Guam's coast, and thus the complexity of physical forcing on the reef. Though we lack ground-truthed residence times in these areas for comparison, the improvements in model skill with the inclusion of these ROMS-based residence time estimates convey that we are capturing important variation in timescales for metabolic modification to occur.

Application to Management
To implement successful OA mitigation strategies on coral reefs, we first need to understand future exposure. In this study we have demonstrated how residence time can be estimated from readily available hydrodynamic information, and how these estimates can increase the predictive skill of benthic ecosystem metabolism models. Using this modeling framework, we can predict how benthic community structure and flow conditions may modulate carbonate chemistry on future reefs, as they are exposed to more acidic oceanic waters and experience altered benthic community composition.
For instance, studies have shown that primary production may be able to offset low pH conditions on coral reefs (Manzello et al. 2012;Unsworth et al. 2012). We can assess this capacity at target reefs by modeling the changes in the carbonate system resulting from changing the benthic composition (for instance, we might expect that increasing photosynthetic community biomass would drive down daytime DIC and buffer pH given the appropriate flow conditions for the signal to manifest). By modeling the effects on the carbonate system of different changes to the benthic community, we can better predict whether management strategies will exacerbate or alleviate OA stress on reefs.
While in situ carbonate system data on coral reefs is undoubtedly useful for interpreting changes over time or differences among habitats, this data is difficult to collect. Field campaigns tend to be limited in spatial and/or temporal scope (Shaw et al. 2012;Koweek et al. 2015a, b;Albright et al. 2015;DeCarlo et al. 2017), and findings are difficult to extrapolate beyond the study domain. Given the importance of the interaction between residence time and benthic composition in driving changes to the carbonate system, we can likely predict areas of favorable carbonate chemistry given information about these drivers. High resolution benthic data are available for many coral reefs worldwide (e.g., National Centers for Coastal Ocean Science [https:// coast alsci ence. noaa. gov/ proje ct/ benth ic-habit at-mappi ng-coral-reefs-flori da-carib bean-pacifi c/], Living Oceans Foundation [https:// www. livin gocea nsfou ndati on. org/ maps/]) and hydrodynamic models can be developed at scales more appropriate for reef scale analyses. The available 2 km ROMS output facilitates estimation of nearshore residence times for a relatively large island like Guam, but this resolution would be inappropriate for smaller islands and atolls. Hydrodynamics models will require downscaling to be appropriately utilized in most reef areas. With information about both benthic composition and local residence time, we can begin to identify probable locations of beneficial changes in carbonate chemistry and coral reef resilience.
Our findings from Guam can be used to demonstrate how information about local residence time paired with benthic composition can be utilized to assess resilience. In Guam, the most extreme changes in TA and DIC relative to offshore conditions were observed in sector 6 (Fig. 2). Sector 6 was characterized by high non-calcifying algal cover and higher calcifier cover relative to the other sectors (Fig. 2). Sector 6 also had the highest residence times (CV methods, Fig. 3). Prolonged exposure to the calcifying and non-calcifying communities in that area likely magnified the signals of calcification and primary production, decreasing TA and DIC, respectively. By contrast, we observe an opposite trend when residence time is lower, minimizing exposure time to benthic metabolism. For example, sector 3 had some of the highest non-calcifying algal cover (Fig. 2a) and therefore we might predict that sector 3 would be a viable location for OA buffering via photosynthesis. However, we observed a lower DIC drawdown in sector 3 than in sector 6, despite the higher cover of photosynthesizing organisms (Fig. 2a), likely due to the relatively short residence times in sector 3 (Fig. 3). Therefore, while sector 3 might initially present a promising area for OA refugia planning due to its abundance of photosynthesizing biomass, in actuality it is likely a poor choice for refugia-based management of reefs since short residence times would prevent the metabolic modification (DIC drawdown) required to combat acidification stress.
Given the importance of residence time in disentangling carbonate chemistry variability and the utility of model-derived residence time estimates, we should invest in the development of hydrodynamic models in complex coastal environments, including coral reefs. We are currently limited by the scale of the hydrodynamic models available, but our statistical model results demonstrate that, despite this limitation, we can make coarse estimates of nearshore residence time that increase predictive skill. To expand this approach, we require hydrodynamic models where we can resolve island mass at the very least and ideally where we can appropriately constrain flow at the scale of nearshore coral reefs. Pairing appropriately scaled hydrodynamic models with benthic habitat maps and oceanic source water chemistry will increase our capacity to improve and expand predictions of nearshore coral reef carbonate chemistry.

Supplementary Information
The online version contains supplementary material available at https:// doi. org/ 10. 1007/ s10498-023-09411-6. as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by 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:// creat iveco mmons. org/ licen ses/ by/4. 0/.