Marine lakes as biogeographical islands: a physical model for ecological dynamics in an insular marine lake, Palau

Marine lakes are emerging ecological and evolutionary natural experimental systems, with genetically isolated resident populations that exhibit extreme population dynamics and rapid phenotypic change. Marine lakes are posited to be marine islands, however, unlike terrestrial islands for which rich models have been developed over the past half-century, we know little of the mechanisms driving changes in marine lakes. This is a critical knowledge gap in efforts to reconcile theory on, or distinguish differences among, island and island-like systems. To reduce this critical knowledge gap, we present a mathematical model describing marine lakes based on a case study of Jellyfish Lake (Ongeim’l Tketau, Mecherchar: OTM), Palau. Empirical data show that marine lakes exhibit delayed and reduced tidal motions, suggesting exchange of a limited amount of water with the neighboring (‘mainland’) ocean. Our model tracks changes in lake level, allowing determination of an exchange rate that is a physical null model for biological colonization and a proxy for colonization distance in island biogeography theory. In addition, we track horizontally averaged in-lake quantities such as salinity and temperature (i.e., marine weather, climate) and stratification (i.e., habitat) — that are known to influence resident species’ distributions and population dynamics — by solving an advection-diffusion equation. We find that weather, ocean conditions, groundwater, and exchanges through tunnels determine the abiotic environment in OTM. By comparing simulations and data, we estimate the difficult-to-measure properties of the surrounding groundwater — the ‘matrix’ in the vernacular of habitat islands — and give a range of realistic values for the effective diffusion coefficient. This coefficient is found to increase in a environment and, therefore, likely influence the intensity of ecological filtering, the diversity of lake habitats, and community structure. - Our paper highlights the potential of marine lakes as model systems in island biogeography.


Introduction
Marine lakes -pieces of seawater entirely surrounded by land -include nearly isolated marine ecosystems that display unique biological and physical characteristics Hauri 1981, Hamner andHamner 1998) akin to islands (Dawson 2016). Over 200 marine lakes exist globally (see Fig. 1), constituting a novel and valuable, but generally poorly known, class of marine habitat. Of the known lakes, approximately 35 have been described in part, primarily with a focus on their geomorphology Hamner 1998, Becking et al. 2011) and genetic (e.g., Dawson and Hamner 2005, Gotoh et al. 2009, Maas et al. 2018 or species diversity (e.g., Fautin and Fitt 1991, Tomascik and Mah 1994, Meyerhof et al. 2016, Rapacciuolo et al. 2019. In contrast, the physical dynamics of marine lakes have been largely ignored, with few exceptions Hamner 1998, Martin et al. 2006), and remain understudied relative to other tropical lakes (e.g., Fazi et al. 2018) and lakes at higher-latitudes (e.g., Mantzouki 2018) and higher elevations (e.g., Kraemer et al. 2017, Hanson et al. 2018). Thus, current understanding is limited, not only of this class of lake but also of how they compare with other lakes and might illuminate novel or rare aspects of limnological systems (Hamner and Hamner 1998). This knowledge gap, in turn, limits understanding of how physical drivers may impact highly valued endemic biological resources, particularly as climate changes (e.g., Kraemer et al. 2017).
Marine lakes, recognized as evolutionary 'natural experiments' for their potential to illustrate patterns and processes influencing diversification of marine populations (Dawson and Hamner 2005) and species (Hanzawa et al. 2012, Dawson 2016, have similar potential to illuminate ecological dynamics. They are a convenient place -akin to closed coastal lagoons, fjords, and human-made lake systems -for studying in situ population dynamics (Dawson et al. 2015) in absence of the high emigration and immigration typical of many open coastal marine systems. Marine lakes thus offer an opportunity to overcome the long-standing challenge in dynamic coastal and oceanic settings of identifying physical boundaries and also delimiting distinct populations and communities (Cowen et al. 2000, Kaplan et al. 2017, Nolasco et al. 2018. Marine lakes, like 'true' oceanic islands before them (Darwin 1859, Wallace 1881, MacArthur and Wilson 1967, Whittaker and Fernández-Palacios 2007 thus offer the potential to clarify the relationships between physical and biological variation, both spatially (e.g., Rapacciuolo et al. 2019) and through time, because of our ability to more accurately identify natural discontinuities in physical and biological processes (Dawson et al. 2015, Dawson 2016. The most substantial limitation on progress is that, unlike terrestrial islands for which rich models have been developed over the past half-century, the physical setting and internal dynamics of these island-like marine systems remain only coarsely described.
Marine lakes are found mostly in sea-level karst landscapes (Hamner and Hamner 1998, Cerrano et al. 2006, Vilibić et al. 2010, Becking et al. 2011, Žic et al. 2012, 2013, though not always (Strelkov et al. 2014). Marine lakes usually occur in clusters (i.e., are archipelagic) and within each archipelagic cluster is a diversity of lakes of varying shapes, sizes, orientations, and distances from the sea. These basic geomorphological differences are associated with differences in gross limnological structure: usually, a lake is either holomictic (and of nearly uniform density at all depths) or meromictic (and strongly density stratified; Hamner and Hamner 1998), which is reflected by gross differences in biological assemblages (Meyerhof et al. 2016, Rapacciuolo et al. 2019. Despite these long-term 'average' characteristics of a lake, a suite of spatial and temporal variations exists within and between archipelagoes Hamner 1998, Becking et al. 2011). Thus, meromictic lakes harbor genetically distinct populations (Dawson and Hamner 2005, Gotoh et al. 2009, Maas et al. 2018) whose dynamics appear to respond to local basin-scale ecosystem variation despite presumably similar weather patterns driven by common regional climatic drivers (Dawson et al. 2001(Dawson et al. , 2015. However, we currently understand little about the processes driving these responses, which include dramatic warming, salinification, and population crashes associated with large-scale climate perturbations, including El Niño and La Niña (Dawson et al. 2001, Martin et al. 2006. The spatial variation and temporal dynamics of marine lakes, illustrated schematically in Fig. 2, are multifaceted. They are influenced by weather patterns, tidal motion, exchanges with the ocean and/or ground-water (Hamner and Hamner 1998), as well as in-lake currents and possibly 'biomixing' attributable (1) Ocean height, tidal range and timing, temperature, and salinity; (2) lake height and tidal range and timing; (3) precipitation; (4) wind velocity, incoming solar radiation, air temperature measured using weather stations (WS), and a vertical profile of water temperature and salinity; (5) exchange or mixing between the lake and surrounding ocean, (6) mixing of the lake at intermediate depths by floodtide waters, and (7) exchanges between the hypolimnion, or monimolimnion, and the mixolimnion.
to large populations of plankton, including medusae (Katija and Dabiri 2009). Marine lakes thus offer insight into biological responses to climate change (Martin et al. 2006) and also complement studies of other somewhat isolated, or bounded, marine systems (Dawson 2016, Wolanski 2017. Here, we develop a model that is a first step towards understanding the impact of exchanges with the ocean and the conditions that determine the stratification, or lack thereof, of a marine lake. We focus on the short-term physical dynamics of "Jellyfish Lake" (Ongeim'l Tketau, Mecherchar (OTM)), the most famous marine lake due to a perennial millions-strong population of golden Mastigias papua medusae that supports a substantial tourist industry, which has been the subject of long-term study since extreme dynamics in 1998-2000 (Dawson et al. 2001, Martin et al. 2006 overturned prior indications of long-term stability (Hamner and Hamner 1998). OTM has been studied monthly for almost two decades, including measurement of various physical limnological parameters and a decade-long synchronous record of local weather patterns from which we draw in this study. We present a model that incorporates weather variations such as rainfall, solar heating, and potentially wind and evaporation, which shape the abiotic environment of the lake. We also focus on novel portions of our model -exchanges with the ocean through porous rock, tunnels, or channels (non-permanent surface inflows) -which tend to be lacking from most extant limnological models and which provide a physical null model and a proxy for biological colonization distance in island biogeography theory.

Methods
To fully describe the physical dynamics of marine lakes, we distinguish three major categories of processes that are modeled and measured. First, tidal processes drive lake level dynamics (#s 1, 2 in Fig. 2). Second, the, largely unknown, subterranean exchanges provide specific sources and sinks that affect the composition of lake water (#5 in Fig. 2). Third, the interior dynamics of the lake redistribute properties such as temperature and salinity and thus influence vertical stratification and the distribution of habitat (#s 3, 4, 6, 7 in Fig. 2). We introduce and integrate each of these effects into our model below.

Lake level dynamics
Exchanges between marine lakes and the neighboring ocean can take place through a combination of porous rock, tunnels, and fissures. The exact geometry and flow rates of those exchanges may not be realistically measured directly. However, the hydrostatic pressure difference between the lake and ocean is easily measured in terms of the time-dependent water level of both as ∆P = g(ρ o H o − ρ l H l ), where g is the gravitational constant, ρ o and ρ l are the ocean and lake densities, respectively, and H o and H l are the ocean and lake surface heights, respectively. All the notation used in our model is summarized in Box 1. We aim to model how the lake level changes for given tidal variations in the ocean (given H o (t)) and supposing known ocean and lake densities. We assume that, on average, the velocity of the flow between the ocean and the lake is proportional to the pressure difference, ∆P, between the two liquid masses, in a manner similar to the celebrated Darcy's law (Darcy, 1856) in porous media: where L is the horizontal distance between the lake and the ocean (estimated as the average straight-line distance between ocean and lake shorelines using the ruler function in Google Earth), κ is the effective permeability of the material separating them, μ is the water's viscosity, and u ex is the average horizontal fluid velocity of the exchanges between the ocean and the lake. While equation (1) is typically applied in porous media, it is also applicable to small tunnels, where the average velocity is also proportional to the pressure gradient provided the flow is laminar, i.e., with a Reynolds number less than approximately 4000. This approach is thus valid for any laminar flow and has also been used to model hydrodynamics of Clipperton Atoll (Jean-Baptiste et al. 2009). The total volumetric flow rate into the lake may then be expressed as: where A v is the area over which water seeps between the ocean and the lake. Given the surface area of the lake, A s , the water level in the lake, H l (t), is then subject to the differential equation For convenience, we define a transport parameter ω as v l s A g A L ρ κ ω µ = with units of 1/time, which encompasses the rate at which exchanges take place between the lake and the ocean. One may therefore think of ω as the inverse of the time taken by the lake to adjust to ocean level changes. The dynamics of the lake level therefore satisfy: We note that as the lake level changes, both the area over which exchanges take place, A v , and the effective permeability, κ, may vary. As a consequence, in general, ω may also be a function of H l , typically increasing slowly as H l increases and more avenues for exchanges become available. However, in certain cases, ω may increase rather abruptly if the lake level reaches a new tunnel or even a channel (surface flow) that serves to suddenly increase exchanges with the ocean.
Although for general time-dependent ocean levels H o (t) equation (2) cannot be solved analytically, the lake level will generally tend to oscillate in a manner Box 1: Summary of notation and parameters used in the model Independent variables: -z: vertical coordinate, in meters. -x: horizontal coordinate, within the tunnels, in meters. -t: time, in seconds. Dependent variables computed by the model: -T: temperature (lake water (no index), within tunnels (tun)).
-C: generic fluid property, here either salinity or temperature.
-C : horizontally averaged property C within the lake.
w : horizontally averaged vertical velocity within the lake.
-u tun : horizontal velocity of the water between the lake and the ocean in a given tunnel.
-H: water level (o for ocean, l for lake).
-φ: lag time between lake and ocean tides, in seconds. Easily observable physical parameters: -μ: viscosity of the water.
-ρ: density of water (o for ocean, l for lake).
-R(t): rate of rainfall, in m/s. -Sun: intensity of solar radiation.
-L: tunnel length, or underground distance from the lake to the ocean.
-κ: effective permeability of the porous rock separating lake and ocean.
u  : velocity vector of the lake water. -u ex : mean velocity of the water exchanged between the lake and the ocean.
-A v : area over which water seeps between lake and ocean. -A s : horizontal area of the lake. -So C,ex : source of C due to exchanges with the ocean. -So C,surf : source of C due to near surface inputs.
-Q: Flux between lake and ocean. -L w : maximum distance traveled by lake water within a tunnel.
-k: tidal frequency, in 1/s. Model parameters to specify: -ω: lake-ocean exchange transport parameter, in 1/s. -D: diffusivity (e: effective, top: at the lake surface, bot at the lake bottom).
-d type : distribution of sources of various types (in, out, rain, rad for radiation).
-β: ratio of lake surface area to tunnel surface area, for each tunnel.
-α: rate at which contents of tunnel water equilibrates with surrounding groundwater.
-T gw temperature of the groundwater around tunnels.
-S gw salinity of the groundwater around tunnels.
resembling that of the ocean, with a smaller amplitude and a time delay. For illustration purposes, we consider the special case of sinusoidal oceanic oscillations, where the ocean level is given by H o (t) = H sin(kt), with k the tidal frequency. The lake level then also behaves in a sinusoidal manner, after a brief transition period, with the long-term lake response given by: where the lag time φ satisfies (ω 2 +k 2 ) cos 2 kφ = ω 2 .
To obtain the results shown in the next section, the differential equation (2) was solved by integrating forward in time using a second order Runge-Kutta method. Ocean and lake levels, used respectively as inputs and outputs to compare against, were measured contemporaneously using Onset HOBO water level data loggers (#U20-001-04) deployed in each lake and at a long-term reference ocean location at Coral Reef Research Foundation, Malakal. The pressure loggers, which have resolution of at least 0.014 KPa (or approximately 0.14 cm), were deployed at 0.9-1.6 m depth from March 1 to March 26, 2012.

Tunnel modeling
We now describe how the subterranean water exchanges between the lake and the ocean are modeled and how they affect properties such as temperature and salinity. We consider systems where a significant fraction of the exchanges takes place through well-defined tunnels. Within each tunnel, tides will cause water to move back and forth, toward, or away from the lake. The temperature, salinity, and velocity of water were measured at the mouth of tunnels using a Seabird Conductivity-Temperature sonde coupled to a low-profile upward-facing Nortek acoustic doppler current profiler (AQUADOPP 2 MHz). The instruments were deployed at each of the five major tunnels in OTM. The data were collected during a 24 hours deployment on October 15, 2014. Similar results were previously obtained in a two-week deployment started on October 23, 2008.
The water within each tunnel has a composition that varies in time. Around ocean high tide, when water flows into lake OTM, it was empirically observed to always have properties that differ from those of ocean water. Around ocean low tide, when water flows toward the ocean, tunnel water tends to be identical to lake water at that same level, after a short delay or recovery time. To capture this time variation, we use a model that assumes that each tunnel extends horizontally, in the x-direction, over a distance L. The contents of each tunnel move away or toward the lake with a velocity u tun (t), taken to be positive when water moves toward the ocean. In addition, the water within the tunnel simultaneously equilibrates its temperature and salinity with that of the groundwater within the porous rock, denoted with the subscript gw. The rate at which this equilibration takes place is given by α, in a manner similar to Newton's cooling law. As a result, the temperature or salinity within the tunnels satisfies a one-dimensional advection equation: Equation (3) was solved using finite differences to approximate spatial derivatives and an upwinding method. Integration in time was performed via a second order Runge-Kutta method. When solving this equation numerically, we initially also considered horizontal diffusion, but for realistic parameters we found that diffusion along the length of the tunnels plays no dynamically significant role, and so we elected not to include it in the final model. A boundary condition must be specified at one end of the tunnel, and the selection of that end is determined by the direction of the flow, which is given by the sign of u tun . If water is flowing toward the ocean, we specify a value of C at the lake, C tun = C tun,l , and if it is flowing toward the lake, we specify it at the ocean end, C tun = C tun,o . Near the lake-end of the tunnel, where the data were recorded, the boundary value should be constrained to be between the value of C tun at the last tide change, when water was flowing into the lake, and the value of C within the lake at the tunnel's depth.
In light of the observation that the water within the tunnel has the same composition as that in the lake, approximately 10 minutes (600 seconds) after the beginning of the ocean-bound flow, we propose to use a boundary condition that varies linearly between the tunnel water when flow stopped and the lake water for the first 600 seconds of the flow, as illustrated in Fig. 3. After that time, the boundary value is the same as that of the lake. We apply a similar approach at the ocean-end of the tunnel.
To calculate the value of C within the entire tunnel, we thus need to specify the values of u tun , L, α and C gw (x). The velocity within the tunnel is assumed to be at this level within the lake, w (z), in such a way that the local volume flux is conserved: where A tun is the cross-sectional area of the tunnel, and β is a dimensionless coefficient. We determine the value of β for various tunnels based on observations of the velocity within a tunnel, with typical values of the order of 1000. Similarly, the appropriate values of α and the properties of the groundwater are based on comparing the model results with observed measurements of salinity and temperature at the mouth of the tunnels and are discussed in the Results section (below). We discuss in the following subsection how to obtain w (z).

Lake interior dynamics
To obtain a useful description of the interior of the lake, we model the time dynamics of the temperature, T, and salinity, S, within the lake. Note that other quantities, such as dissolved oxygen level, chlorophyll, or pollutant concentration, etc., may be treated in a similar manner. However, those properties would also generally require complex biological modeling. To focus on the novel aspects of our model, exchanges between lake and ocean, we restrict our attention to temperature and salinity. These properties generally satisfy the advection-diffusion equation, potentially supplemented by sources and sinks (referred to as sources from here on, for conciseness) accounting for exchanges with the ocean: where C can stand for either T or S. Here the divergence is taken in three dimensions, u  is the velocity of the water, and D e is an effective diffusivity that includes turbulent water motion and may depend on both time and depth. To reduce the dimensionality of the model, we consider horizontally averaged properties, denoted with a bar, C , defined as where A s (z) is the area of the lake at a given height z, which can be accurately estimated from satellite data.
Here we set the level z = 0 at the bottom of the lake, and z increases upwards, so that z = H l at the surface.
To obtain an equation of the averaged quantities, we integrate equation (4) horizontally. To allow further progress, certain assumptions must be made. First, we assume that over the scale of the entire lake, the water velocity is not correlated with the quantities we study. This implies that there is no large-scale vertical circulation in the lake. Rather, all fluid motions occur on a small scale so that uC uC =   . Second, we assume that the properties at the edges of the lake are not systematically different than they are near the center, so that the average of vertical derivatives can be approximated by the derivative of averages: Our model may then be expressed solely in terms of horizontally averaged quantities: We note that this same equation may also be obtained via a careful mass balance (Montroy, 2013). We solved equation (7) numerically using an explicit, upwinding scheme for the advection term, and an implicit Crank-Nicholson method for the diffusive term. Here, w is the averaged vertical velocity, positive upwards. This velocity is positive when the lake level increases and negative when it decreases. To describe it more precisely, we introduce the vertical distribution with depth of incoming flow into the lake, d in , and that of outgoing flow d out . Each distribution accounts for all inflow or outflow, respectively, due to tunnels and seepage, and so satisfies We note that these two distributions may not be identical, as incoming flow can quickly adjust its depth once it reaches the lake to reach its level of neutral buoyancy. The distribution function therefore refers to an averaged distribution after this readjustment has taken place. The averaged vertical velocity at any given level therefore satisfies, by mass conservation, where the choice of which distribution to use depends on the direction of the flow, as given by the sign of ( ) In addition to influencing the vertical fluid velocity, exchanges with the ocean also act as a source of heat and salinity. At any given height z, these averaged sources, So C,ex , take the form where the inward or outward distribution is chosen depending on the sign of w, and C tun stands for the salinity or temperature within the tunnels connecting the lake and the ocean, as described above. Another type of source is due to the vicinity of the upper surface of the lake, where inflow of fresh water (rain) and heat (solar radiation) may take place. These inflows are spread over a few meters near the surface. The surface source (or sink) of salinity due to rain is where R(t) is the rate of rainfall, assumed to contain no salt, in length per unit time, and d rain is the depthdistribution of rainwater, taking into account run-off from surface groundwater. This distribution is chosen to be uniform over the top two meters of the lake. We note here that the rate of rainfall is at most of the order of 10 -7 m/s, which is small compared to the rate of change of the depth of the lake due to tidal effects, which is of the order 10 -5 m/s. The surface source of heat, which affects the lake temperature, includes the heat influx due to rainfall, and that due to solar radiation where T rain is the temperature of the rain, d rad is the distribution of solar radiation, taken to decay exponentially over the first four meters, and Sun(t) is the intensity of solar radiation. Our simulations have shown that over time scales exceeding a day, the exact choice of the rain temperature and intensity of solar heating are negligible compared to the influence of the surface air temperature. For conciseness, we therefore only present results where the rain temperature is set to the air temperature and where solar heating is neglected. These effects could be included in a study focused on time scales of the order of hours. Finally, boundary conditions at the top and bottom of the lake must be specified for each property tracked. We note that the effect of rainfall on salinity is accounted for by the source term So S,surf rather than as a boundary condition, which would be equivalent to assuming that all incoming rainfall enters the lake at the very surface. The temperature is assumed to match that of the surroundings at both the surface and bottom of the lake where T bot (t) and T air (t) are, respectively, the temperature of the rock at the bottom of the lake and the air temperature above the lake. The air temperature can be obtained from weather data and changes over the course of each day. The rock temperature, however, will be taken as constant, matching the observed temperature at the bottom of the lake. We note that the stated boundary conditions neglect evaporation, which can affect both salinity and temperature. However, estimates of evaporation rates give an order of magnitude of approximately 10 -8 m/s, which renders its effects negligible compared to those of other factors included here. Given the ocean level over time and weather data, such as air temperature and rate of rainfall, all the parameters derived in the previous section may be completely specified except for one (as will be described in the Results section). The exception is the effective diffusion (or mixing) coefficient as a function of depth. While there exist empirical relations between wind velocity and mixing (e.g., Riley andStefan 1988, Botte andKay 2002) and descriptions of the mixing rate decay with depth (Kaufman et al. 1991), these were mostly obtained in deeper lakes or in the oceans and for a significantly thicker upper layer (above any pycnocline). These relationships also involve complicated dependencies on weather conditions. To focus on the characteristics that are specific to marine lakes, we follow here a simpler approach.
To estimate the diffusion coefficient, we consider that turbulent fluid motion due to waves, currents, and swimming organisms account for most of the mixing within the lake, while molecular diffusion is present but negligible. We also anticipate that mixing due to in-lake motion (excluding the influence of tunnels which are treated separately) is smallest near the bottom of the lake, since it is where water is the most stagnant, and largest near the surface where the wind acts as the main source of fluid motion. Finally, we expect that density variations, in particular the presence of a pycnocline, will greatly reduce mixing.
Over the scale considered here, diffusivity is expected to remain nearly constant above the pycnocline and to remain nearly constant below it, though with a different value. Below the pycnocline, the overall temperature and salinity profiles of the lake are expected to vary rather slowly, and the model should be able to capture this base state. We propose a diffusivity based on variations in the lake density and write the logarithm of the diffusivity as a linear interpolation between the logarithm of a prescribed surface diffusivity, D top , and that of a prescribed deep-water diffusivity, D bot : This form has the advantage of being sensitive to abrupt density variations such as those found near the pycnocline. We thus need to specify only two diffusion coefficients, D top , D bot , as the density may be obtained as a function of the salinity and temperature (Millero and Poisson 1981). A more detailed discussion of the form of the diffusion coefficient may be found in Montroy (2013).
To determine the values of D top , and D bot , we compare simulation results with vertical profiles of salinity and temperature measured at meter intervals from the lake surface to lake bottom. These profiles were obtained using an Hach DS5 sonde with conductivity, temperature, dissolved Oxygen, pH, Chlorophyll a, and depth (i.e., pressure) probes. These profiles also serve as initial conditions for the model. The initial conditions for the simulations were taken on August 16, 2010.

Full model with simulation of tropical storm Bopha
As a short-term 'ecological' example, we considered the response of our model to short-lived, vigorous changes. We focus on a tropical storm (Bopha) that took place in November 2012. For three days, significant rainfall and violent winds battered the Palau area. To model this event, we therefore imposed a rate of rainfall of 4.5 x 10 -7 m/s, corresponding to a little less than 4cm of rain per day. We also decreased the air temperature at the surface from 31.1 °C to 29.7 °C, matching the measured surface lake temperatures before and after the storm. The ocean properties were left unchanged, as they did not vary much over such a short time scale. We observed the response of the model after three days of identical conditions and compared it with measurements taken after the storm.

Impact of tunnels on mixing
Our model allows us to determine the effect of the presence of tunnels on the dynamics within the lake. For lake OTM, four main tunnels have been identified, at depths of 1 m, 6 m, 9 m, and 11 m. As a first approximation, and based on direct observations, we estimate that each tunnel accounts for an equal fraction of exchanges with the ocean. We further assume that together all four tunnels account for half of the total exchanges with the ocean. The remaining half of the lake-ocean exchanges are evenly distributed over all the rest of the lake, which goes down to a depth of approximately 33 m. Using this information, we may build the distribution of inflow into the lake and outflow from the lake. For comparison, and to characterize the impact of the presence of the tunnels, we may simply run our model with the same rate of exchanges with the ocean, but replacing the inflow/outflow distribution concentrated at the level of the observed tunnels with a uniform distribution of outflow and inflow. We discuss the impact of the presence of tunnels in the next section.

Results
Following the structure of the Methods section, we present our results on the lake level dynamics, tunnel modeling, and lake interior dynamics.

Lake level dynamics
We show in Fig. 4 numerically computed solutions of the dynamic level of OTM using measurements of the ocean level and of an initial lake level. To illustrate the impact of the transport parameter ω, we show computed results for a range of values for this parameter.
Larger ω clearly results in a greater lake amplitude and in a shorter lag. If we assume an exchange rate that is independent of lake level, the best choice for lake OTM appears to be ω ≈ 5.5 x 10 -5 /s (Appendix S1). Fig. 5 shows measurements of the velocity, salinity, and temperature made at the entrance of one of the largest tunnels of lake OTM (at a depth of 6 m at the lakeside). We show in Fig. 5a the measured computed velocity for a single tunnel (black curve). Based on this data, we find that choosing β = 1500 yields good agreement between model (red curve) and tunnel velocity for velocities toward the lake (negative values). However, the measurements of velocities away from the lake are significantly lower. In fact, the measurements recorded far more water coming into the lake than exiting the lake. Over time, this would lead to a sustained rise of the lake level, which does not happen. We therefore postulate that this is likely due to the position of the measurement device, slightly outside the tunnel itself, which leads to fairly accurate measurements of the concentrated jet of water coming from the tunnel into the lake but captures poorly the drainage of the lake into the tunnel which takes places  in a much less concentrated manner. We therefore choose the value of β based on measurements of the inflow only. In the case of lake OTM, measurements were taken at the mouths of four tunnels. For all other depths that do not have well-defined underground tunnels, we estimated the value of u tun to be ten times smaller than within well-defined tunnels. This value was chosen because the velocity of water flowing into the lake through porous rock should be significantly less than that through a hollow tunnel. As a result, exchanges at those levels are less significant as their contents always closely resemble lake water.

Tunnel modeling
We note from the measurements shown in Fig. 5b and 5c that, for OTM, tunnel water is never identical in content to ocean water. The highest salinity reached is 30.84 kg/m 3 , while the salinity of the ocean is about 33.5 kg/m 3 . Moreover, at high tide the temperature of the incoming flow drops below 29 °C, while both the lake and the ocean temperatures approach 31 °C. We therefore conclude that ocean water does not directly reach the lake within a single tidal cycle. In contrast, when water flows toward the ocean at low tide, the tunnel water is seen to be nearly identical to lake water. There is a short but clearly visible delay of 10 to 20 minutes (600 to 1200 seconds) between the time when the tide turns and the time at which the contents of tunnel water begins to change. This delay is presumably due to the low velocities as the tide turns and to the ebb and flow of the same water moving back and forth. More specifically, water that just entered from the lake into the tunnel is likely to reenter the lake upon a tide change, and vice versa. Our boundary condition is chosen to capture this delay and is illustrated in Fig. 3.
The tunnel length L can be estimated from geographical observations and must be compared to L w, the maximum distance traveled by water within the tunnels: If L w > L, water coming from the ocean would reach the lake at every tidal cycle, and the contents of tunnel water would remain the same after a finite transition period. Measurements of salinity (Fig. 5b) and temperature (Fig. 5c) of tunnel water do not show such a plateau, indicating that the water does not travel from the lake to the ocean within a single tidal cycle, so that L w < L. The exact value of L is thus relatively insignificant so long as it remains greater than L w . For lake OTM, we therefore set L = 200 m, a value consistent with geographical measurements and greater than L w for any of the tidal cycles considered here.
The rate of exchanges between tunnel and groundwater, α, is more difficult to evaluate, but comparison with tunnel data, as shown in Fig. 5, indicates that α . 1.17 x 10 -4 /s yields good agreement. These simulations also showed that the results are not very sensitive to the choice of α. More critically, we need to estimate the temperature and salinity in the groundwater surrounding the tunnel, C gw , for which we have very little data. For the salinity, a simple linear interpolation between the lake and ocean, consistent with salt slowly diffusing in from the ocean, yields good agreement with tunnel measurements. However, the temperature of the water exiting the tunnel is less than that of both the lake and the ocean. We thus selected to use a groundwater temperature profile that varies across a relatively thin boundary layer at both the lake and ocean ends of the tunnel to a cooler rock temperature of 28.5 °C in the center, with a boundary layer thickness of 1/10th that of the entire tunnel. The temperature and salinity profiles used in our model are shown in Fig. 6.
Based on the preceding general considerations and specific observations, all the parameters required to simulate tunnel flow have been set. Using the ocean level over time, we may thus compute the lake level over time and obtain l dH w dt = . This in turn yields the tunnel velocity, u tun = w. From the data collected on tunnel inflow velocity, we determine β (for example β = 1500 for the data shown in Fig. 5a), and we may therefore compute the salinity and temperature of the tunnel water. The computed salinity and temperature are compared to the measurements made at the mouth of the tunnel in Fig. 5b and 5c. We first note that at low tide, when water flows into the ocean, between 400 and 800 minutes (24000 and 48000 seconds) in Fig. 5, the water in the tunnel is the same as in the lake, which our model automatically recovers. The onset of the variations associated with inflow into the lake is also well captured by our model. The salinity of the tunnel water is actually predicted very well by our model, with a similar time profile and maximum reached. The temperature of tunnel water is also generally well predicted by our model, although the exact shape profile is somewhat different. This is most likely due to our ignorance of the true temperature distribution of the groundwater, which in turns affects Figure 6. Salinity and temperature profiles used to model the properties of groundwater. The profile for salinity assumes that diffusion leads to a linear transition between ocean and lake values. The profile for temperature assumes that the groundwater is cooler than either lake or ocean, and it is constant outside a relatively small boundary layer near either water mass. the temperature of the incoming water. We note a small discrepancy in the duration of the inflow, with our model predicting a longer inflow duration, which could be due to differences between the tidal height in the immediate vicinity of the lake compared to where measurements were taken. However, overall, our model predicts the tunnel velocity, salinity, and temperature variations to a satisfactory degree of accuracy, which will allow us to use it as source terms in the equations governing the dynamics of the lake's interior.

Lake interior dynamics
The only parameters left to determine in our model are the values of the diffusion coefficient at the top and bottom of the lake. The value of D bot is determined by noting that within the data available to us, the temperature and salinity of the lower layer remained effectively constant over a time scale of weeks. For the model to maintain nearly constant bottom temperature and salinity over time scales of up to a month, we found the requirement that D bot remains below a threshold value. Since small values correspond to virtually no mixing, we use as a bottom diffusivity the upper bound of that range and set D bot = 8.33 x 10 -7 m 2 /s. For comparison, the molecular thermal diffusivity of water is approximately 1.5 x 10 -7 m 2 /s, and that of salt is considerably smaller.
The value, or range of realistic values, of D top was determined by examining weekslong simulations. In Fig. 7a and 7b we show the salinity and temperature profile over time, assuming no rainfall and a constant surface air temperature, for D top = 8.33 x 10 -5 m 2 /s. This choice clearly results in a significant erosion of the density profile. The temperature being fixed at the surface, most of the changes in temperature occur around a depth of 10 m. The salinity profile, subject to a zero-derivative at the top and bottom, flattens at the surface first, and the extent of the region of constant salinity increases over time. In the absence of any restoring mechanism, the pycnocline is thus eroded in a matter of weeks for this choice of D top , which thus appears to be larger than the diffusion constant in the upper layer corresponding to regular conditions. Here the temperature is kept fixed at the top and bottom. The lower most curve is the initial data and modelled temperature curves become progressively higher as days pass. (b) In-lake salinity. Here the derivative is set to zero at the top and bottom. The left-most curve is the initial data and modelled salinity curves move progressively to the right as days pass. (c) and (d) modelled in-lake temperature and salinity of lake OTM after one week, for various surface diffusion coefficients. (c) In-lake temperature vs depth. (d) In-lake salinity vs depth.
To obtain a range of realistic upper diffusivity, we proceed to determine the largest diffusion coefficient that does not result in a significant change of temperature and salinity profile over a one-week period. Fig. 7c and 7d show the temperature and salinity profiles after one week for various values of the upper diffusion coefficient, D top . As discussed in the next section, the range 5-10 m deep is strongly affected by the presence of tunnels, and so it is not the most indicative of the effects of the diffusion constant. On the other hand, significant differences are visible near the top of the lake. For D top = 1.67 x 10 -6 m 2 /s, and smaller values of D top (not shown), only very minimal changes are present near the surface over a period of one week. However, choosing D top = 8.33 x 10 -5 m 2 /s shows clear differences within a week, and so does choosing D top = 1.67 x 10 -5 m 2 /s, though to a lesser extent. We thus conclude that a realistic range of diffusion coefficients is that 10 -4 ≤ D top ≤ 8.33 x 10 -5 m 2 /s. Local conditions such as the wind speed and the extent of bio-mixing induced by swimming organisms are all expected to influence the diffusivity. This parameter appears to be the most determinant of all the parameters appearing in this model, and more precise quantification of its dependence on various factors based on direct, targeted measurements will help improve the accuracy of the model. Fig. 8a and 8b show the model response for three different choices of the upper layer diffusion coefficient in storm conditions, including heavy rainfall. A relatively small diffusion coefficient of D top = 8.33 x 10 -6 m 2 /s is seen to be insufficient to capture the changes near the surface, indicating that the strong winds increased mixing in the upper layer. A better fit to the measured data is achieved using a diffusion coefficient of D top =5.0 x 10 -5 m 2 /s, for which the salinity profile agrees with post-storm measurements. The temperature profile also agrees better, but the significant drop in temperature in the 5-10 m depth is not captured by the model. Even a larger diffusion coefficient does not reduce the temperature significantly in this range. However, this range is where tunnels are most influential, as discussed in the prior section. We therefore postulate that the groundwater temperature during this storm was probably also reduced slightly, by drainage of rain through the karst, and that this lower temperature water would have been transported into the lake via the tunnels.

Effect of tunnels on mixing
In Fig. 8c and 8d we show the temperature and salinity profiles after one week with, and without the presence of tunnels, for two extreme values of the diffusion constant. The tunnels periodically bring saltier, colder water to the lake and this is most clearly seen at depths ranging from 5-10 m. The temperature at this depth is clearly reduced by the inflow of cold water from the tunnels, which thus contribute to increase the heat exchanges between the groundwater and the lake. The salinity level is also increased by the presence of tunnels, as they enhance the transport of saltier water into the lake. We note that the tunnels maintain the presence of the pycnocline where the density changes abruptly. In fact, they may have been determinant in setting the depth of the pycnocline as the water they provide to the lake has properties that approach those of the bottom of the lake.

Discussion
Island biogeography theory has since its earliest incarnations drawn heavily from the physical aspects of terrestrial islands, most notably their geographical isolation and size Wilson 1963, 1967). Subsequent developments integrated additional physical categories such as small versus not-small islands (Niering 1963), mainland-island geometry (Taylor 1987), the physical stages of oceanic island formation, erosion, and subsidence (Whittaker et al. 2008), different geological types of island (e.g., Novosolov andMeiri 2013, Ali 2017), and whether islands are archipelagic, meta-archipelagic, or neither (Whittaker et al. 2018). Our goal here, therefore, is likewise to begin to describe the basic physical properties that vary among and may define the island-like qualities of insular marine lakes.
Through a first-principles approach, we have obtained a complete and manageable model of a marine lake that includes exchanges with the ocean through a karst matrix including tunnels. Our model predicts the variations of the lake level over time and can therefore describe the lag of lake oscillations relative to the tidal cycle of the neighboring ocean. We also compute the dynamics of the horizontally averaged temperature and salinity over the depth of the lake, as well as the flow rate, temperature, and salinity of the water within well-defined tunnels connecting the lake to the ocean. These physical characteristics determine the physical connectivity of marine lakes and also their internal structure, which ultimately manifests in the coarse categorization of lakes as holomictic (mixed) or meromictic (stratified) Hamner 1998, Colin 2009) and in turn strongly influences their biological characteristics (Meyerhof et al. 2016, Rapacciuolo et al. 2019. Our model uses as input some well-known physical quantities (water density and viscosity, gravitational acceleration; see Box 1 for a full list of parameters). In addition, some easily measured inputs are required, such as the distance to the ocean, weather data, and oceanic level over time. Some inputs are more difficult to determine but have been estimated for certain marine lakes, such as the distribution and importance of tunnels (Hamner and Hamner 1998). We also inferred that a realistic range of diffusion coefficients is 1.67 x 10 -6 ≤ D top ≤ 8.33 x 10 -5 m 2 /s. Such factors are nominally easily described or applied for other marine lakes (e.g., Becking et al. 2011, Žic et al. 2012, Strelkov et al. 2014) and related ecosystems (e.g., Vilibić et al. 2010). Using such inputs, our model is thus broadly applicable to additional marine lakes in Palau and also other karstic marine lakes globally (e.g., Dawson et al. 2009).
Nonetheless, some parameters used in our current model may not be estimated directly. For these we demonstrated the utility of an inverse-problem approach -consisting of estimating parameters based on their effects on the model computations -which may also be applied broadly. In Palau, observations of changes in the surface temperature and salinity may be used to estimate the effective diffusion coefficient within the lake, and we used this approach to obtain a range of realistic diffusion coefficients near the surface. Our present model does not predict the diffusivity in varying weather conditions, but Fig. 8 shows that more precise predictions would require better estimates of the local diffusivity as a function of conditions such as wind strength. The lack of information regarding the groundwater composition and temperature also limits the accuracy of the predictions that can be made by our model. On the other hand, by examining the results of the model, it is possible to gain information on both the groundwater and the diffusion coefficient, two quantities difficult to measure directly. Since both factors act most strongly in different locations, it is possible to determine a realistic range for each. In the case of tropical storm Bopha, the groundwater appears to have had a temporarily reduced temperature, at least in the ground in contact with the upper tunnels, and the diffusion coefficient appears to have increased in the upper portion of the lake.
We also used inverse reasoning to estimate the rate of exchanges between the lake and the ocean, ω, by comparing the predicted and recorded lake levels. This rate may then be used to quantify the volumetric rate at which water flows between the ocean and the lake and to estimate parameters such as the effective permeability of the material separating them. Notably, the parameter ω is an important estimator of the physical connectivity of marine lakes and, thus, may be an important predictor of biological connectivity sensu Cowen et al. (2000). Although averaged estimators have ambiguous biological implications and inferences of marine dispersal are more commonly based on Lagrangian models (Siegel et al. 2003), ω has potential as a first estimate of relative connectivity of populations in lakes relative to the ocean given a species' dispersal kernel. This rate of exchange between the lake and the ocean, ω, may provide a better estimation of the oceanographic colonization distance and insularity of these marine systems than does the overland distance from ocean to lake.
The ability of the inverse-problem approach to better delimit 'hidden' characteristics that influence marine lakes is of great value. The question of connectivity depends not only on the distance and duration of transport ('isolation by distance') but also on the environmental characteristics that differ between two points ('isolation by resistance'; Sexton et al. 2014). Predictions of the flow rate and contents of tunnels can also be used to infer the composition of the surrounding groundwater, or water caught in pockets of some tunnels, through which marine larvae and other propagules presumably must pass to enter a lake. By comparing model predictions and data recorded near the mouth of tunnels, we deduced that the groundwater salinity profile was approximately linear between the lake and the ocean, and that groundwater temperature, away from the immediate vicinity of the lake and ocean, quickly reached a value lower than that of either water mass. Osmotic stress due to differences between ocean and lake salinity, therefore, likely is gradual rather than stepped but, nonetheless, increases the relative isolation of a lake from the ocean. In contrast, the cooling of larvae when subterranean may slow development and imbue relatively higher dispersal duration, thus decreasing the relative isolation of a lake from the ocean, sensu Brown (2014). These physical estimates (and biological inferences) may be compared to future measurements of the effective diffusivity in various circumstances and of groundwater properties (though the latter measurements are considerably more difficult to obtain). The range of systems that fit the definition of marine lakes is, however, broad, and some may equally be suited to simulation preliminarily using models already applied to freshwater lakes. By exploring the range of lakes and models and their intersections we will surely learn much about the physical dynamics of marine lakes, a strategy we encourage.

Conclusion and outlook
Our long-term purpose in creating this model is to understand how changes in weather and climate affect the coupled physical and biological dynamics of marine lakes. Through the simulation of a storm, we have shown that short-term severe weather can influence lake structure (Fig. 8) but may not have a long-lasting effect (Hamner and Hauri 1981). We expect that the model will inform the dynamics of analogous semi-isolated and bounded marine systems, such as the lagoonal basins of Palmyra and Clipperton Atoll (Gardner et al. 2011, Sanvicente-Añorve et al. 2017. Moreover, our model provides a foundation for any user to investigate the importance of additional factors, such as the potential for 'biomixing' by large jellyfish populations (Katija and Dabiri 2009), the diurnal migration of micro-and mesozooplankton (Houghton et al. 2018), or changes in the number of human visitors. The model can also be extended in a straightforward manner to include components other than salt, such as the concentration of dissolved oxygen or pollutants, and thus, for example, allow for estimates of the accumulation, or flushing, of sunscreen brought in by human visitors (Bell et al. 2017). The model is sensitive enough to account for short-lived storms, yet it is also sufficiently stable as to allow investigations of long-lasting effects (years).
We intend to use the approach described here to study lake dynamics over longer time scales, including asynchronous dynamics of different lakes (Dawson et al. 2015) and how density stratifications may be established, or eroded, over year-long events. For example, lakes that are usually meromictic may occasionally over-turn and temporarily become holomictic (e.g., Goby Lake, Koror), moderately deep lakes may variously be oxic or anoxic at the bottom (e.g., T-Lake (Lake 10)), the pycnocline may shoal and epilimnion warm excessively in association with severe El Niño-La Niña events (e.g., OTM (Dawson et al. 2001)). Such events have potential to extirpate endemic populations and, in a period of increasingly extreme climate, demand our understanding and, if possible, mitigation. The isolation and stratification of the lakes has been responsible for generating novel populations and subspecies but may also put them at risk. To better understand the long-term dynamics of marine lakes, we also may study their in-lake physical dynamics as sea-levels rise, as they did during formation after the Last Glacial Maximum. Stratification and mixing are apparent properties of marine lakes that are related to population dynamics and patterns of biodiversity (Martin et al. 2006, Rapacciuolo et al. 2019. Understanding the formation and dynamics of marine lakes is key for understanding the origins and fate of their endemic biodiversity.