Summary of NCRMP carbonate chemistry data and reef metabolic processes analyses
Hannah Barkley
August 2018
1 Carbonate chemistry samples
To start, below are several maps and summary figures showing the spatial and temporal scales and variability of the NCRMP carbonate chemistry data.
1.1 Map of Pacific and Atlantic water samples
Interactive map of water samples collected between 2010 and 2017. Points are colored by \(\Omega_{ar}\) (values can be seen by hovering over individual points).
1.2 Table of locations
Table of regions, region codes, location codes and island names.
| RegionCode | Region | LocationCode | Location |
|---|---|---|---|
| AMSM | American Samoa | OFU | Ofu and Olosega Islands |
| AMSM | American Samoa | ROS | Rose Atoll |
| AMSM | American Samoa | SWA | Swains Island |
| AMSM | American Samoa | TAU | Ta’u |
| AMSM | American Samoa | TUT | Tutuila |
| CNMI | Commonwealth of the Northern Mariana Islands | AGR | Agrihan |
| CNMI | Commonwealth of the Northern Mariana Islands | AGU | Aguijan (Goat Is.) |
| CNMI | Commonwealth of the Northern Mariana Islands | ALA | Alamagan |
| CNMI | Commonwealth of the Northern Mariana Islands | ASC | Asuncion Island |
| CNMI | Commonwealth of the Northern Mariana Islands | FDP | Farallon de Pajaros |
| CNMI | Commonwealth of the Northern Mariana Islands | GUA | Guam |
| CNMI | Commonwealth of the Northern Mariana Islands | GUG | Guguan |
| CNMI | Commonwealth of the Northern Mariana Islands | MAU | Maug |
| CNMI | Commonwealth of the Northern Mariana Islands | PAG | Pagan |
| CNMI | Commonwealth of the Northern Mariana Islands | ROT | Rota |
| CNMI | Commonwealth of the Northern Mariana Islands | SAI | Saipan |
| CNMI | Commonwealth of the Northern Mariana Islands | SAR | Sarigan |
| CNMI | Commonwealth of the Northern Mariana Islands | TIN | Tinian |
| MHI | Main Hawaiian Islands | HAW | Hawai’i (Big Island) |
| MHI | Main Hawaiian Islands | KAH | Kaho’olawe |
| MHI | Main Hawaiian Islands | KAU | Kaua’i |
| MHI | Main Hawaiian Islands | LAN | Lana’i |
| MHI | Main Hawaiian Islands | MAI | Maui |
| MHI | Main Hawaiian Islands | MOL | Molokai |
| MHI | Main Hawaiian Islands | NII | Niihau |
| MHI | Main Hawaiian Islands | OAH | Oahu |
| NWHI | Northwestern Hawaiian Islands (Papahanaumokuakea Marine National Monument) | FFS | French Frigate Shoals |
| NWHI | Northwestern Hawaiian Islands (Papahanaumokuakea Marine National Monument) | KUR | Kure Atoll |
| NWHI | Northwestern Hawaiian Islands (Papahanaumokuakea Marine National Monument) | LIS | Lisianski Island |
| NWHI | Northwestern Hawaiian Islands (Papahanaumokuakea Marine National Monument) | MAR | Maro Reef |
| NWHI | Northwestern Hawaiian Islands (Papahanaumokuakea Marine National Monument) | PHR | Pearl and Hermes Atoll |
| PRIA | Pacific Remote Island Areas | BAK | Baker Island |
| PRIA | Pacific Remote Island Areas | HOW | Howland Island |
| PRIA | Pacific Remote Island Areas | JAR | Jarvis Island |
| PRIA | Pacific Remote Island Areas | JOH | Johnston Atoll |
| PRIA | Pacific Remote Island Areas | KIN | Kingman Reef |
| PRIA | Pacific Remote Island Areas | PAL | Palmyra Atoll |
| PRIA | Pacific Remote Island Areas | WAK | Wake Atoll |
1.3 Summary of island-level carbonate chemistry
Cross-Pacific \(\Omega_{ar}\) data for 38 islands collected during three sampling periods (2010-2011, 2012-2014, and 2015-2017). Boxplots show \(\Omega_{ar}\) values nearshore (benthic and surface) samples, and triangles indicate \(\Omega_{ar}\) of paired offshore samples (where collected). 
2 Offshore references
In order to conduct reef metabolic analysis, it is essential that nearshore reef carbonate chemistry be compared to an offshore oceanic end member. Paired offshore-onshore data exist for 25 island-years (where an “island-year” represents a set of water samples collected at an island during a particular year; e.g. Jarvis-2015). However, the majority of reef samples do not have paired offshore samples, which means that a large portion of the carbonate chemistry data set cannot be used for process inference (or I need a more creative solution for obtaining offshore values!).
2.1 Criteria for offshore reference samples
First, to ensure that offshore water samples used for process analysis represented true oceanic end member values, I selected reference samples using three criteria:
- the sample was collected >1km from the coast (for American Samoa, Main Hawaiian Islands, Mariana Islands, and Pacific Remote Islands) or the 20m depth contour (Northwest Hawaiian Islands, where the reef flat extends far beyond the extent of land),
- the sample was collected at a location where seafloor depth was >100m
- the calculated \(\Omega_{ar}\) for the sample was within two times the median absolute deviation of offshore \(\Omega_{ar}\) values for that island and year.
Mean values of salinity-normalized (to S = 35) offshore TA and DIC were paired with each onshore sample from the same island and sampling period.
2.2 Discrete offshore samples vs. GLODAP climatology
One approach I’ve tried using to address missing offshore data is using GLODAP climatological TA/DIC/\(\Omega_{ar}\) values as an offshore reference. This superficially seemed like a great solution because GLODAP climatology is gridded, global, and based on spatial interpolation of actual observations.
However, when I compared offshore TA/DIC/\(\Omega_{ar}\) data from discrete samples to values predicted by GLODAP, GLODAP data differed from observations. The top plot below shows offshore \(\Omega_{ar}\) vs. GLODAP climatology for islands and sampling years for which offshore data exist. The bottom plot shows the difference between measured offshore \(\Omega_{ar}\) and GLODAP \(\Omega_{ar}\). The mean difference in \(\Omega_{ar}\) is 0.11 with a range of -0.18 to 0.36. Because the magnitude and direction of the bias is not consistent, it is not straightforward to apply a simple bias correction. As a result, I don’t believe GLODAP climatology is a viable alternative for offshore reference values.
2.3 Neighbor island offshore references
The other solution that I’ve come up with to address missing offshore data is to pair offshore values from islands and years where samples exist to neighboring islands with no offshore data. I assigned neighbor island offshore values to islands where the distance between islands was <300km and where the percent change in GLODAP \(\Omega_{ar}\) was <3% between the location of the offshore reference sample and that of the at onshore samples at the neighboring island. This increases the number of island-years from 25 to 47.
3 Inferring reef processes
Because I don’t have any hydrodynamic data or estimates of residence time, I can’t directly calculate NEC or NEP.
Instead, I calculated the “calcification-dissolution signal” (CDS) and “photosynthesis-respiration signal” (PRS), which are essentially rate-less analogs of NEC and NEP, from offshore-onshore gradients in salinity-normalized TA and DIC to infer the net sign and significance of local reef processes. CDS and PRS were derived for each onshore water sample using the following equations:
CDS = -\(\frac{1}{2}\)(\(nTA_{onshore}\) - \(nTA_{offshore}\))
PRS = -(\(nDIC_{onshore}\) - \(nDIC_{offshore}\)) - \(\frac{1}{2}\)(\(nTA_{onshore}\) - \(nTA_{offshore}\))
The major disadvantage of this approach is that, because CDS and PRS are not scaled to residence time, the magnitude of the signal cannot be treated as a rate (and therefore cannot be use to directly compare islands/sectors). Only the sign of the change in TA/DIC on the reef (i.e. negative or positive relative to offshore) can be interpreted.
The other major constraint with this approach is that I’m assuming that all changes in TA and DIC are entirely due to biological processes. This may be a reasonable assumption in many of the uninhabited small atolls, but could be a poor assumption in places where other oceanographic processes are occurring (e.g. upwelling in the equatorial PRIAs) or where coastal runoff is a factor.
3.1 Sector-year calcification-dissolution signal from bulk chemistry offsets
I tested the significance of the CDS and PRS for each sector-year by conducting one-sample Wilcoxon sign rank tests against the null hypothesis of a zero median difference between offshore and onshore samples. I adjusted p values using a false discovery rate correction to control for multiple comparisons.
The boxplot below shows CDS calculated between an offshore reference and nearshore samples for all sector-years with paired offshore-reef water sample data. Boxplots are shaded by the sign and significance of each signal, with stars indicating a significant difference from zero. Sector years with a significant calcification signal are shaded pink, those with a significant net dissolution signal are shaded blue, and those not significantly different than zero are shaded gray.
3.2 Map of sector-level calcificaton-dissolution bulk signals
Here’s a map of what the CDS sector data look like. Sectors shaded pink have a net positive CDS, gray is non-significant, and blue is net negative. Sectors with insufficient data have been removed from this plot. The east side of Guam (“GUA_EAST_OPEN”) is the only sector with a “net dissolving” signal (however, I’m not convinced that this isn’t just due to river output on the east side).
3.3 Sector-year photosynthesis-respiration signal from bulk chemistry offsets
Boxplots below show PRS data plotted by sector-year. Boxplots are shaded as significantly positive (net photosynthesis; green), significantly negative (net respiration; yellow), or not significantly different than zero (gray). The number of samples is too low and the variability too high for this to really tell us anything. 
3.4 Slopes
Another method for evaluating reef metabolic processes that is residence-time independent is calculating the slope of the salinity normalized TA-DIC property-property plots to evalute the balance of NEC and NEP. Below are nTA-nDIC property plots for the sector-years where there were more than 5 samples and for which there was a significant linear relationship between nTA and nDIC (p < 0.05). Points are colored by sample type: red = benthic, blue = surface, and green = offshore.
4 Drivers of spatial and temporal carbonate chemistry patterns
4.1 CDS and PRS as a function of environmental variables and benthic composition
The most tempting and, unfortunately, least responsible way to compare carbonate chemistry to benthic cover is to relate the magnitude of the TA/DIC drawdown to percent cover of benthic groups. However, this is an “illegal manuever” because it ignores residence time as a driver of the bulk offshore-inshore change in TA and DIC (and the relationship between benthic cover and CDS/PRS isn’t incredibly compelling anyway). Points are colored by the significance of that sector-year: blue points are from a sector-year with a significant net positive CDS, green are non-significant, and red are net negative.
4.2 Categorical CDS significance
The more responsible way to compare CDS to environmental and benthic data is to categorize the data into groups based on sign and significance (as shown in the plots below).
I’ve evaluated the significance of the relationship between benthic cover and calcification-dissolution signal significance by building a GAM with calcification-dissolution signal treated as an ordered categorical variable. The resulting model does a poor job explaining variance in the data, and only sediment emerges as a potential driver. This doesn’t seem to be the way to go.
##
## Family: Ordered Categorical(-1,3.98)
## Link function: identity
##
## Formula:
## CalcDiss.rel ~ s(CORAL, k = 3) + s(TURF, k = 3) + s(CCA, k = 3) +
## s(OmegaAragonite, k = 3)
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.5778 0.2748 13.02 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(CORAL) 1.819 1.967 4.278 0.145
## s(TURF) 1.000 1.000 1.178 0.278
## s(CCA) 1.778 1.950 3.058 0.165
## s(OmegaAragonite) 1.000 1.000 1.979 0.159
##
## Deviance explained = 10.2%
## -REML = 44.462 Scale est. = 1 n = 614.3 Relationships between TA-DIC slope and benthic cover
I’ve also examined the relationship between nTA-nDIC slope and benthic cover. The figures below show the slope plotted against percent cover of benthic groups. 
4.4 GAM analysis of normalized slopes driven by individual benthic groups
GAMs highlight coral and halimeda percent cover as the most important drivers (with a surprisingly high 84% deviance explained)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## N.Slope ~ s(CORAL, k = 3) + s(TURF, k = 3) + s(MA, k = 3) + s(CCA,
## k = 3) + s(HAL, k = 3) + s(SED, k = 3) + s(OmegaA, k = 3) +
## s(temp, k = 3)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.02882 0.06021 17.09 5.19e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(CORAL) 1.974 1.997 6.497 0.0065 **
## s(TURF) 1.645 1.872 0.763 0.4289
## s(MA) 1.000 1.000 0.007 0.9348
## s(CCA) 1.538 1.785 0.549 0.5538
## s(HAL) 1.000 1.000 0.571 0.4591
## s(SED) 1.790 1.955 1.667 0.1987
## s(OmegaA) 1.000 1.000 6.482 0.0195 *
## s(temp) 1.000 1.000 0.572 0.4586
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.504 Deviance explained = 68.5%
## GCV = 0.18287 Scale est. = 0.1124 n = 31## Global model call: gam(formula = N.Slope ~ s(CORAL, k = 3) + s(TURF, k = 3) + s(MA,
## k = 3) + s(CCA, k = 3) + s(HAL, k = 3) + s(SED, k = 3) +
## s(OmegaA, k = 3) + s(temp, k = 3), data = gSl.nona, na.action = "na.fail")
## ---
## Model selection table
## (Int) s(CCA,3) s(COR,3) s(HAL,3) s(OmA,3) s(SED,3) s(tmp,3) s(TUR,3)
## 19 1.029 + +
## 83 1.029 + + +
## 147 1.029 + + +
## 23 1.029 + + +
## 20 1.029 + + +
## 55 1.029 + + + +
## df logLik AICc delta weight
## 19 5 -9.871 33.6 0.00 0.226
## 83 5 -9.324 33.7 0.13 0.212
## 147 7 -7.134 34.2 0.59 0.168
## 23 6 -9.293 34.7 1.04 0.134
## 20 6 -8.886 34.7 1.06 0.133
## 55 7 -6.556 34.8 1.16 0.127
## Models ranked by AICc(x)4.5 GAM analysis of normalized slopes driven by functional benthic groups
I’ve grouped the benthic categories into functional bins - coral, calcifying algae, and macroalgae - and rerun the GAM. Here, calcifying algae and macroalgae emerge as the strongest drivers.
##
## Family: gaussian
## Link function: identity
##
## Formula:
## N.Slope ~ s(CORAL, k = 3) + s(CAL.ALG, k = 3) + s(ALG, k = 3) +
## s(OmegaA, k = 3) + s(temp, k = 3)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.02882 0.06501 15.82 3.26e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(CORAL) 1.842 1.973 5.298 0.008712 **
## s(CAL.ALG) 1.000 1.000 0.379 0.543673
## s(ALG) 1.000 1.000 1.131 0.298018
## s(OmegaA) 1.120 1.223 12.209 0.000713 ***
## s(temp) 1.000 1.000 2.708 0.112687
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.422 Deviance explained = 53.7%
## GCV = 0.16896 Scale est. = 0.13102 n = 31## Global model call: gam(formula = N.Slope ~ s(CORAL, k = 3) + s(CAL.ALG, k = 3) +
## s(ALG, k = 3) + s(OmegaA, k = 3) + s(temp, k = 3), data = gSl.nona,
## na.action = "na.fail")
## ---
## Model selection table
## (Int) s(ALG,3) s(CAL.ALG,3) s(COR,3) s(OmA,3) s(tmp,3) df logLik AICc
## 13 1.029 + + 5 -9.871 33.6
## 29 1.029 + + + 5 -9.324 33.7
## 14 1.029 + + + 7 -7.968 35.2
## 30 1.029 + + + + 6 -8.915 36.2
## 15 1.029 + + + 6 -9.801 36.5
## 31 1.029 + + + + 6 -9.328 37.1
## delta weight
## 13 0.00 0.323
## 29 0.13 0.303
## 14 1.53 0.150
## 30 2.61 0.088
## 15 2.84 0.078
## 31 3.45 0.058
## Models ranked by AICc(x)5 Pacific vs. Atlantic carbonate chemistry
I’m in the beginning stages of integrating Atlantic data into these analyses. Incorporating Atlantic data adds the additional factor of seasonality, which is less of a concern for Pacific data because water samples are collected during the same season for each region. The boxplot below shows all Atlantic and Pacific data plotted by season (all are northern hemisphere seasons except for American Samoa).
6 Additional analyses
Here are a few other things I’ve done that indirectly relate to this project.
6.1 CO2 buoy data
Ridgeline plots showing the hourly, monthly, yearly, and year-month (daily time series data binned by month) distribution of pCO\(_{2}\) data from CO\(_{2}\) buoys.
6.1.1 Ala Wai
6.1.2 CRIMP
6.1.3 CRIMP2
6.1.4 Kaneohe
6.1.5 Kilo Nalu
6.1.6 WHOTS
6.1.7 Cheeca
6.2 Tutuila diurnal suite
Instrument data from diurnal suite deployments at Fagatele and Aunu’u (Feb-Mar 2015). (A) Map of instrument deployment locations at Tutuila. (B) SeaFET pH time series. (C and D) Diurnal suite instrument deployments at Fagatele (C) and Aunu’u (D) include pressure (ADCP), pH (SeaFET - blue; PUC discrete samples - red), temperature, and salinity (CTD). 