Holocene Closure of Lib Pond, Marshall Islands

Well-preserved sediment from closed water bodies of atolls such as Lib Pond are rare opportunities to reconstruct the past regional climate, which pieced together across a latitude and longitude range identify the range of movement patterns of wider scale climate phenomena such as the Intertropical Convergence Zone (ITCZ) and El Niño Southern Oscillation (ENSO). We conducted the first physico-chemical survey of Lib Pond, a shallow, closed-water saline lake located on remote and difficult to access Lib Island in the Marshall Islands at 8° 18′ 48.99″ N, 167 22′ 51.90″ E in the Pacific Ocean, in July 2009. We performed a bathymetric survey, recorded salinity, dissolved oxygen, pH, and temperature profiles, monitored the tidal variability, and conducted a vegetation survey surrounding the lake. From bathymetric data we calculated the lake volume, which we used to estimate the lake's salt budget, and ultimately the residence time of water in the lake basin. We took a series of sediment cores from the lake, cores which indicate Lib Island's changing environment and climate. Radiocarbon measurements determined sediment age, and reveal significant mixing over the last 2 ka of deposition. We conclude that prior to 3 ka, Lib Island was an atoll with a central lagoon connected to the open ocean, which was then closed off from the open ocean to form the brackish system that exists today. We predict that the sediment accumulation in Lib Pond evident today will continue. As seawater is inhibited from exchanging with fresh water, Lib Pond will become a shallower lake with increasingly fresh water.

. Inflatable row boats, used for the pond survey and coring, are shown with the lake shoreline in the foreground. The shoreline was highly vegetated and not easily accessible on foot.
The view is looking toward the SE, from near the north end of Lib Pond.

Prior Mentions of Lib in the Scientific Literature
Lib Island was briefly summarized by renowned 20 th century tropics researcher Francis Raymond Fosberg in June 1988 [7]: Apart from that, we know of no description of Lib Island in the literature despite an extensive search. With this geochemical study, we have aimed to undertake the study Fosberg suggested.

Kwajalein Atoll Annual Rainfall vs. SOI
The total annual rainfall over Kwajalein Atoll , dotted red line) superimposed on the Southern Oscillation Index , solid black line). Given that the precipitation appears to closely match SOI oscillations Kwajalein, and therefore Lib, has local characteristics influenced by larger, regional trends in the equatorial Pacific Ocean. Figure S7. Sources: [8] and the Southern Oscillation Index, 2009. Tables   Kwajalein Atoll Tide Tables, July 2009 Date

Lib Pond Bathymetry Data
The GPS coordinates were imported from the Garmin GPS units to GoogleEarth.
However, they needed to be offset given the GPS error, which placed several pond depth coordinates beyond the body of Lib Pond, and past the shoreline on land. Since the original pond orientation top was not true north, the lake satellite image from Google Earth was reoriented and rotated in Mathematica using the Lambert Azmuthal projection. There was only one good fit within the pond dimensions, given that we traversed the pond edge and hugged the shore at a section on the SE and NE shore. The outline of Lib Pond was drawn in Mathematica and exported back into GoogleEarth, along with the original depth pairings. For images showing this, please contact the corresponding author.

Satellite Image Alignment
The satellite image has a crosshair reticle at a specified point in the lake (Latitude, Longitude). It was aligned such that the crosshair was centered over the coordinates (0,0) in a Lambert Azmuthal projection in Mathematica, which has a built in function for displaying that particular projection. Thus, georeferencing error is negligible. Below, the red points indicate points at which the depth was taken, and the pink line indicates the extent of the floodplain from the lower left quadrant and out. The blue line are the perimeters of Lib Pond.

Lib Pond GPS Coordinates and Depth Pairings
This original dataset was used to construct the bathymetric map of Lib Pond (positions are before the offset correction from GPS error.) The transects were composed of 263 GPS coordinate and depth pairs. Data points were taken on July 21 st and 23 rd at comparable tidal periods.

Location of some of Sediment Cores and Geochemistry Profiles
For a zoomed-in satellite image of Lib pond, labeling core sampling and water profiling sites with available GPS data, please contact the corresponding author (image is a GoogleEarth screenshot).

Figure S10. pH and DO % Correlation from Geochemistry Profiles
Graphic representation of the positive correlation between measured pH values and DO% saturation, with DO% reported on a log scale.

Filter System
We used a filter to check and collect plankton and other minute organic material in Lib Pond. Below is one such filter: Figure S11.

Algal Mat Example
The below example shows an algal mat which was often at the surface of the cores, including several at Lib. Here, a picture of a similar core from Mejit Island, also in the Marshalls, is shown as a visual example (taken on July 28 th , 2009, evidenced by the 090728 date.) It is approximately 0.5 cm thick. It is held up by a plastic sectioner. Figure S14.

Conventional Bathymetry Background
The conventional approach to bathymetry is a statistical fitting of a smooth surface through a set of depth points. Numerous depth data are collected in a grid or by crisscrossing transects [9]. A bathymetric map is then created from statistical techniques commonly available in GIS software such as kriging [10] or laplacian smoothing [11]. However, these statistical methods only work if sufficient depth data exists because they do not assume any knowledge about the lake other than the depth values and locations.
In this case, field conditions prevented the collection of a sufficient large data set (and depth transects) to use those conventional methods. A standard interpolation technique by itself produces a bathymetric map with unrealistic features: knife-edge ridges and sharp angles which are inconsistent with field measurements and experiences with similar types of lakes.
A prior methodology measures the central deep point of a body of water and just several other points in all four directions, orthogonal to each adjacent direction, forming a cross [12].
"Synthetic" radial transects are modeled upon the observed transects to create a smoothed transect with a realistic lake bottom (e.g. no rough or jagged edges). However, Lib Pond contains smaller features and is more irregular than the wetland basins and we did not collect perpendicular linear transect data (a prerequisite to creating radial transects). Additionally, the deepest part of Lib Pond is unknown and the required laser surveying equipment was unavailable.
Other conventional methods that do not employ the above statistical techniques include multispectral satellite imagery which produces a bathymetric map with centimeter-resolution [13] [14] and radial transects [12]. Yet satellite images of Lib Pond at multiple frequencies are not available.

Boundary Slope Fitting Methodology
Boundary slope fitting creates a model of near-shore lake regions using information nearby regions with similar characteristics and can be used for any shallow body of water. It is considerably more cost effective than other conventional bathymetric methods in terms of equipment and in terms of the number of required depth measurements. Original bottom data is sufficient to determine lake bathymetry if supplemented with data from the lake edges. Hence, boundary slope fitting adds modeled data points around the lake perimeter, which in turn are Thus lake regions without depth data can be indirectly determined from adjacent, known lake slopes. Field observations and unrecorded depth data were used to make the assumption that the slope characteristics are relatively homogenous in a given section where no direct depth data was taken: similar shore types and bottom conditions imply similar slope characteristics.
We used exponential functions since polynomials are not flat at the center of the lake as you extend from the shoreline (most lake bottoms, including Lib, are asymptotically flat). The unnavigable shallow floodplain is located to the bottom left and marked by a black line.
Edge types indicate similar slopes.

Created Contours from Boundary Slope Fitting, Before and After
For boundary slope fitting contours of Lib Pond shown against the satellite image of the pond, please contact the corresponding author. The channel on the lower left is accurately reflected in the bathymetry to the right of the channel opening. These contours are compared against contours from a standard interpolation of the original data set.

Black and White Luminosity vs. Depth (Remote Sensing)
As another possible way to gauge the depths of Lib Pond, the satellite image was converted from color to monochrome. The idea is that the original depth data were taken across pixels with different luminosities, and that a correlation between depth and pixel value might exist.
The color satellite image of the pond was converted to gray scale to see if there was a correlation between depth and color. But in this case there was no apparent correlation when pixel luminosity values are plotted versus the depth data. There was no correspondence between the color and the depth. As you can see from the graph on the right, there is little correlation, indicated by the wide variation of depths associated with each value of pixel luminosity. This is likely due to the different bottom types and specular reflection of sunlight off the surface of the pond.
Therefore, remote sensing is not helpful in further determining the depth of Lib Pond. This is not surprising since color differences in the water column are likely due to bottom type change or particulate and organic matter, rather than by depth change. Figure S56. Determine similar shoreline and bottom conditions of pond using satellite photo. Assemble linear transects of depth data.
Curve-fit transects that go up to shore using a rational exponential function. You need transects that go up to shore for each shoreline type present in the pond.
Create modeled datapoints along normal vectors from the pond edge, using the curve fits from transects. Length of vectors depends on transect data and shoreline curvature.
Now you have combined modeled data with measured data. The modeled data gives appropriate boundary conditions for interpolating/triangulating the measured data from deeper parts of the pond.
Instead of an interpolation, can use a new application of an existing computational geometry method: Delaunay Triangulation to get volume estimate. In essence you create a triangular mesh.
Conventional techniques to get bathymetry can still be applied to hybrid data set (modeled/observed) if desired. These interpolations can be used to compute the volume.

Volume Calculation Discussion
There are additionally several ways to calculate the lake volume after a bathymetric map is created. One could take modeled data and use contour plotting, kriging or other statistical techniques to calculate the volume but we used Delaunay triangulation. It and other computational geometric methods are widely used in finite element modeling and computer graphics, both of which require modeling a surface with a triangulated mesh at a high level of accuracy, and are a particular way of taking a set of data points in 2D and dividing them up into a set of triangles [15] [16].
Results from Delaunay triangulation were compared to a conventional kriging technique as implemented in ArcGIS. To better visualize the difference between the original depth data and the added shoreline points (i.e. modeled data)a difference map was created from a kriging interpolation done in ArcGIS after boundary slope fitting was applied. The volume difference between a Delaunay triangulated surface and kriging is negligible for the original data point set and the modeled one.  Figure S58.
The boundary slope fitting variation compared against the volume variation. This is done by putting the for the volume and the curve fitting, and plotting the ratios against each other on different axes, for 5,10, and 20% standard deviations. Note that curve fitting variation is the parameter variation, since the parameters compose, and can be derived from, the curve fits.
This graph shows that there is roughly 4:1 ratio of curve fitting error to volume error. In other words, if the curve fitting error is known to be 10%, then the resulting volume error from the curve error would only be 2.5%. So even if the curve fits were off (highly unlikely; see prior analysis) the pond volume would not be off by much as the curve fitting is off. The fit to the line is: .
Bathymetric maps comparing kriging using the original dataset (right) and the dataset including the modeled points with boundary slope fitting (left). The depth transitions are more realistic with the addition of modeling, as seen from the color display.

Figure S59.
A comparison of the bathymetric map of the original data (right) versus including the modeled data points (left), shown using a Triangulated Irregular Network (TIN) created in ArcGIS. The region beyond the channel is visible on the left but not the right, demonstrating that boundary slope fitting makes the resulting depths adhere to the physical properties of the pond.

Figure S60.
A difference map comparing depths using kriging with their corresponding depths using the bathymetric map from Delaunay Triangulation. White areas indicate no change. (Ignore the color outside the boundary).
The map is shown on the next page: Figure S61.

Error Analysis of Volume Estimation
An error analysis of the data and methods is presented to demonstrate the validity the boundary slope fitting method and the volume and bathymetry of Lib Pond. First, there is error in the raw data used to compute the lake bathymetry (GPS coordinates and associated depths, Each GPS coordinate was randomly scattered according to an error distribution. The Monte Carlo simulation of 500 trials, which produced 500 different lake volumes, had a mean volume value of 170,205 m 3 and a standard deviation of 594.474 m 3 . Thus, lake volume did not vary substantially due to GPS error. A larger error source comes from the method used to estimate the volume (in this case boundary slope fitting). To test the effect of variance on volume, Gaussian distributed random parameters for the exponential function curve fits were added. (Gaussian noise with a mean of zero and a standard deviation of 10% of the value found in the curve fit for each parameter a,b,c, and d). This results in 10% in the slope parameter, in the depth parameter, etc., quite a large cumulative error to apply, as seen in the resulting slope differences.
The modeled error is far higher than the actual residuals from boundary slope fitting and is intended to simulate lake regions not well-approximated by nearby transect data. Again, a Monte Carlo simulation of 500 trials produced 500 lake volume values with a mean of 170,040 m 3 and a standard deviation of 4648.33 (substantially higher than the GPS error but nevertheless relatively low.) Expressed as a percentage, the ratio of the standard deviation divided by the mean lake volume is 2.7%. Despite a simultaneous 10% difference in each parametersuch as curve slope and lengththe end lake volume was essentially identical. Monte Carlo results using parameter standard deviations of 5% and 20% had similar effects to the 10%.
The effect of boundary slope fitting on the resulting volume can be quantified by a ≈ 4:1 ratio when the boundary slope fitting errors are plotted against the volume error. In other words, one quarter of the error from the curve fit gets passed onto the volume calculation. This is consistent with the fact that the boundary slope fitting method only added about 25% of the volume of the lake. The volume calculation is insensitive to boundary slope fitting due to modeling shallow regions of the lake as opposed to deep regions.
Finally, any echo sounder depth reading will have penetration error. Hence, the definition of the bottom using sonar is dependent on sediment porosity, and varies depending on echo sounder strength and frequency. Thus knowing the specifications of the echo sounder, including the measurement resolution, is more important for shallow bodies of water as than the open ocean. A typical GPS gamma error distribution. A similar figure is seen in [17].  The three graphs below show that the boundary slope fitting error is relatively small. The histograms fits a Gaussian "bell" curve. Volume is in cubic meters (m 3 ). (Figure S66) Using 10% standard deviations. The same boundary slope fitting variation compared against the volume variation, done for 5% deviations ( Figure S67) and 20% deviations (Figure S68) of the parameters. The R 2 value of the 5% SD is 0.999934, for example. Figure S66.

Water Volume Estimation
The calculated Lib Pond volume increases by ~25% modeling the bathymetry using boundary slope fitting. Computing the volume after kriging versus Delaunay triangulation makes little difference to the resulting volume. However, the lake bottom consists of extremely porous sediment so the volume estimate only accounts for water outside of the sediment, what can be called the open water column.
Therefore, taking the sediment porosity into account when calculating the water volume the system contains is a high remaining source of uncertainty. Since 8.5+ m of sediment were recovered at the deepest point, of which the first 4 m were very porous (> 0.5, or 50% water), the total lake sediment volume is likely greater than the open water volume. A single value of porosity is likely appropriate for the entire sediment column, since the depth is on the order of meters (8 m) instead of a km (i.e. 1-2 km, which would be modeled linearly) or many kms (i.e. 5-10 kms, which would be modeled exponentially) [18], but permutations involving an exponential drop off in porosity should also be tested.
An extremely porous orange gelatinous sediment layer that was 90%+ water was likely represented as part of the water column since the echo sounder "read" through it. Field experiments demonstrated that the depth sounder was inaccurate and overestimated depth on the order of 0.5 m when tested in shallow water (< 1 m, where the bottom could be visualized.) Thus, the current estimate for lake water volume is almost certainly an upper bound. Multiplying the surface area of the lake by 0.5 m and subtracting it from the volume estimate subtracts 60,346 m 3 (ignoring significant figures, for the sake of the example). For boundary slope fitting using Delaunay triangulation, the volume would change from 169,900 m 3 to 109,554 m 3 , a baseline volume that is almost certainly a lower bound. Therefore a range of Lib Pond water volumes is ~110,000 m 3 to ~170,000 m 3 .
Thus the current estimated lake volume is really an open water volume calculation; there could be a substantial amount of water in the sediment. Thus the current state of Lib Pond is consistent with a lake in a late stage of sediment accumulation that is filling in. We believe the shallow floodplain region will eventually expand to cover the entire lake based on the high sediment accumulation rates inferred from the Lib Pond sediment cores, which would reduce the volume estimation of open water in Lib Pond.