This toy model explores whether deep reefs provide a refugia from storm damage. Deeper reefs have less wave damage, but also less light to grow. The code is available in GitHub.
Sets the x axis range for all plots. The model is designed for 2 to 40 m depths.
More light increases coral growth, up to a point (determined by half-saturation). Turbid waters decrease the available light (Kd turbid). We model two scenarios, clear and turbid for comparison. The growth is supplemented by heterotrophic growth. This model assumes annual average light conditions. See eReefs BGC model for typical Vertical attenuation and Light intensity values.
Meaning: Daily Light Integral at the surface, representing the total photosynthetically active radiation received per day averaged over the year. This value should account for typical cloud cover. GBR typical values range from 36–49 mol m⁻² d⁻¹ depending on region and season. Reduce this value for cloudier regions or years.
Meaning: Diffuse attenuation coefficient for relatively clear oceanic water. Expressed as the exponential rate of light loss per metre. See eReefs BGC for typical values. Typical 0.05 to 0.10 m⁻¹ offshore.
Meaning: Diffuse attenuation for turbid shelf water scenario. Typical 0.15 to 0.40 m⁻¹ on nearshore reefs.
Meaning: DLI where coral growth is half of maximum. Based on GBR Acropora studies showing strong growth at 12–16 mol m⁻² d⁻¹. Using ~14 as half-saturation means growth approaches maximum around this range.
Meaning: a normalised annual growth potential at full light. Units are arbitrary here but consistent across the model.
Meaning: non photosynthetic feeding and other energy sources that provide a baseline growth even at low light. Typical 0.02 to 0.10 of gmax.
This section models common storms (that mostly affect rubble) and typical cyclones that do not strip a reef.
Meaning: height of the waves during trade-wind and monsoonal storms. Typical wave height of 2.5 to 4.5 m offshore (medium confidence).
Meaning: period of the dominant waves. The longer the wave period the deeper the wave penetration. Typical values 6 - 9 s (medium confidence).
Meaning: how frequently the moderate storms occur. This uses a Poisson distribution so a setting of 1 results in a probability of 63% per year. The more severe the storm settings the less frequent they should be.
Meaning: height of the waves during cyclone events. Typical wave height of 4 to 7 m offshore. Can reach 10 m in deep water (medium confidence).
Meaning: period of the dominant waves. The longer the wave period the deeper the wave penetration. Typical values 10 - 16 s (medium confidence).
Meaning: how frequently the severe cyclones occur. This uses a Poisson distribution. The more severe the storm settings the less frequent they should be.
This section determines how much water movement is needed to move rubble (which breaks coral), and break coral directly.
Meaning: Critical shear current from waves needed to mobilise rubble. Estimated to be 0.3 to 0.5 m s^-1 (low confidence).
Meaning: How sharp is the rubble threshold. 0 - represents hard threshold, 0.2 represents a broad threshold.
Meaning: What fraction of the coral is removed when the rubble moves. Note: The rubble areas only represent a small fraction of reefs. This represents the average over the whole reef (low confidence in the typical value).
Meaning: The shear velocity, caused by waves, needed to break the coral into rubble. Acropora would have a much lower breakage velocity than a non-branching species.
Meaning: How sharp the threshold is? A value of 0, means values below the threshold have no effect. A value of 1 means that the breakage occurs over a wide range of values.
Meaning: Fraction of coral lost when breakage occurs. This is effectively the amount of coral lost during a cyclone.
This section models extreme cyclones (Cat 5) that can strip a reef to the rock. Under these conditions recovery is slower and has a lag. We include this category because while these are very infrequent these cyclones will produce larger waves that penetrate deep, undermining potential slow growing deep refugia.
Meaning: Water velocity from waves that results in stripping back to the substrate, i.e. not just some breakage, but very high levels of damage. This should be set to a higher velocity than the typical storm threshold.
Meaning: How sharp is the threshold. 0 - very sharp, 1 - occurs over wide range of values.
Meaning: Fraction of coral that is lost when a catastrophic event occurs.
Meaning: Number of years before the reef returns to normal recovery growth rate.
Meaning: What fraction of the normal growth rate do we see during a lag period.
Light decays exponentially with depth and growth follows a saturating curve, with a modest heterotrophic floor that supports some growth at low light.
Wave driven near bed orbital velocity decays with depth and intersects rubble and breakage thresholds at characteristic depths that depend on wave height and period.
Combining event rates with threshold exceedance gives annual probabilities that fall off with depth, with breakage dropping faster than rubble motion.
The difference between growth and expected storm loss peaks at an intermediate depth rather than cancelling out across all depths.
What this model shows. This interactive tool explores how coral reef resilience varies with depth. The core idea is that there may be an optimal depth range—a "refugia"—where corals are deep enough to avoid the worst storm damage, but shallow enough to receive sufficient light for growth. The four plots illustrate this trade-off:
How to use this tool. Adjust the sliders in the left panel to explore different scenarios. Try changing water clarity (Kd), storm frequency, or wave height to see how the optimal depth shifts. The CSV export button lets you download the current profiles for further analysis.
Technical details: Variables and units. Depth z in metres. Surface daily light integral DLI0 in mol photons m⁻² d⁻¹ (total PAR received per day, averaged over the year including typical cloud cover). Diffuse attenuation coefficient Kd in m⁻¹. DLI at depth DLI(z) in the same units as DLI0. Half-saturation light Ik in mol photons m⁻² d⁻¹—the DLI at which growth is half of maximum. Maximum annual growth potential gmax in normalised units per year. Heterotrophic growth floor g_floor as a fraction of gmax. Wave height H in metres and peak period T in seconds. Angular frequency omega = 2π/T in s⁻¹. Wavenumber k in m⁻¹ from the dispersion relation. Gravitational acceleration g in m s⁻². Near bed orbital velocity ub(z) in m s⁻¹. Damage thresholds ucrit_rubble, ucrit_break, ucrit_cat in m s⁻¹. Transition steepness parameters sigma_r, sigma_b, sigma_c in m s⁻¹. Event rates lambda_i in events per year. Annual probabilities P_rubble(z), P_break(z), P_cat(z) are unitless. Loss severities Mr, Mb, Mc are unitless fractions of cover. Post-catastrophe recovery lag lag in years and growth during lag g_lag as a fraction of normal growth. Net recovery margin R(z) in the same normalised units as gmax.
Light decay with depth. Daily light integral at depth is modelled as DLI(z) = DLI0 × exp(-Kd × z). DLI0 is the surface DLI representing the average daily PAR received, accounting for typical cloud cover and seasonal variation. Kd is the diffuse attenuation coefficient that captures water clarity. Lower Kd means clearer water, deeper light penetration. This exponential form (Beer-Lambert law) is well-established for broadband PAR in the ocean. Typical GBR surface DLI values range from 36–49 mol m⁻² d⁻¹ depending on region, with central GBR around 44–48 mol m⁻² d⁻¹.
Growth response. Potential growth follows a saturating relationship with light: g_photo(z) = gmax × DLI(z) / (DLI(z) + Ik). Ik is the half-saturation DLI—the light level where growth is half of gmax. GBR studies show strong Acropora growth at 12–16 mol m⁻² d⁻¹, with little additional benefit above this range. Using Ik ≈ 14 mol m⁻² d⁻¹ means growth approaches its maximum around this light level. To include non-photosynthetic feeding (heterotrophy) we add a floor, so actual potential growth becomes g(z) = g_floor × gmax + (gmax - g_floor × gmax) × DLI(z) / (DLI(z) + Ik). This captures appreciable growth at moderate depths and saturation in shallow water.
Storm kinematics and damage. Linear wave theory gives the dispersion relation omega² = g × k × tanh(k × z), which we solve for k at each depth and period. The near bed orbital velocity magnitude is approximated by ub(z) = (π × H / T) / sinh(k × z), which decays with depth. We compare ub to velocity thresholds. The probability that a given event causes damage is mapped with a logistic curve p = 1 / (1 + exp(-(ub - ucrit) / sigma)), where sigma sets the sharpness of the transition. With event rates lambda_i we compute annual probabilities using a Poisson process: P(z) = 1 - exp(-sum_i lambda_i × p_i(z)). Breakage is split into a non-catastrophic portion and a catastrophic portion by subtracting the catastrophic exceedance share.
Catastrophic damage and recovery. A high threshold ucrit_cat with steepness sigma_c and severity Mc represents events that remove most cover and impose a recovery delay. If the rate of catastrophic events at a depth is lambda_cat(z) and the recovery lag is lag years, the expected fraction of time spent in the lag state is phi(z) = (lambda_cat × lag) / (1 + lambda_cat × lag). During lag, potential growth is reduced to g_lag times its normal value. The effective growth is g_eff(z) = g(z) × (1 - phi) + g(z) × g_lag × phi. Expected loss combines categories as loss(z) = Mr × P_rubble + Mb × P_break + Mc × P_cat. The net margin is R(z) = g_eff(z) - loss(z). As catastrophic frequency or lag grows, phi approaches one, suppressing recovery even if a severe event does not occur in the current year.
Limitations. This is a toy model for exploring concepts, not a predictive tool. It does not include: spectral or directional waves, reef crest wave transformation, infragravity energy, wave-current interactions, topographic sheltering, internal waves, temperature/bleaching effects, sediment dynamics, recruitment supply or survival, larval connectivity, species-specific mechanics, or successional changes. The model also assumes constant conditions over time and does not capture interannual variability. Use this tool to explore sensitivity and to shape hypotheses, then move to more detailed models when you need predictions.
Event frequency uses a Poisson process. If the mean rate is lambda = 1 event per year and an event, when it happens, certainly causes damage at a given depth, then the probability of at least one event in a year is 1 - exp(-lambda), which is approximately 0.632. To get close to certainty within a year you need lambda much larger than 1 or a deterministic annual event. The catastrophic lag captures that even without a new severe event, a recent one can keep the system in a suppressed state.
Running...