""" ocean-darkening.py -- the Science Journaling Club's own simplified light-and- volume model for the field note "Twenty One Percent of the Ocean Got Darker and Almost Nobody Noticed" (Data Story, Oceanography). THIS IS NOT THE PAPER'S ANALYSIS. Davies, T. W. & Smyth, T. (2025). Darkening of the Global Ocean. Global Change Biology 31(5), e70227. DOI: 10.1111/gcb.70227 Davies and Smyth did the satellite work: twenty annual composites of MODIS Aqua's diffuse attenuation coefficient at 490 nm, Kd(490), on a 9 km grid, 2003 to 2022, ARIMA-filtered and then fitted with quantile regression to the median. They report the areas. Everything below is the club rebuilding the optics from the textbook equation and then doing arithmetic on THEIR reported areas. Where our numbers look like theirs, that is a sanity check. Where they differ, assume we are the ones who are wrong. ============================================================================= THE MODEL ============================================================================= PART A -- BEER-LAMBERT ATTENUATION Downwelling irradiance in water falls off geometrically with depth: I(z) = I0 * exp(-Kd * z) I0 is irradiance just below the surface, z is depth in metres, and Kd is the diffuse attenuation coefficient in units of inverse metres. Kd is the single number that says how murky the water is. Kd = 0.02 is the clearest open gyre. Kd = 0.5 is a river plume. The euphotic (sunlit) depth is conventionally the 1% light level. Set I(z)/I0 = 0.01 and solve: z_1% = ln(100) / Kd = 4.60517 / Kd PART B -- THE NONLINEARITY Differentiate: dz/dKd = -ln(100) / Kd^2 The depth lost for a fixed increase in Kd goes as the inverse SQUARE of Kd. Clear water is fragile: the same absolute amount of added murk costs a clear gyre hundreds of times more metres than it costs an estuary. The FRACTIONAL loss, by contrast, is scale free. A 10% rise in Kd always removes 9.09% of the euphotic depth, whatever Kd you started from. Both facts are in the same equation and they pull in opposite directions. PART C -- SECCHI DISK CONVERSION A Secchi disk is a white plate on a rope. You lower it until you lose sight of it and write down the depth, Z_sd. The standard conversion is Kd = 1.7 / Z_sd from Poole & Atkins (1929). The constant is empirical and genuinely varies: published values run from about 1.4 to 2.0 in marine water and can reach 3 in turbid estuaries. We use 1.7 throughout and show the spread the other constants would give. Substituting gives a result worth memorising: z_1% = ln(100) * Z_sd / 1.7 = 2.709 * Z_sd Euphotic depth is about 2.7 times your Secchi reading. PART D -- LOST LIT VOLUME The paper reports shoaling by AREA in nested bands. We un-nest them into disjoint bands, assign a representative shoaling depth inside each band, and multiply out: lost volume = sum over bands of ( area_band * depth_band ) Representative depths are chosen at or below the midpoint so the estimate leans low. The open-ended top band (>100 m) is charged at exactly 100 m, which is certainly an underestimate. A seeded Monte Carlo then samples the representative depth uniformly inside each band to put an interval on it. ============================================================================= ASSUMPTIONS, STATED PLAINLY (A1-A8) ============================================================================= A1 Kd is constant with depth. It is not. The real water column stratifies, and a subsurface chlorophyll maximum can sit below a clear surface layer. Davies and Smyth make the same assumption and say so. A2 Kd(490) at one wavelength stands in for the whole photosynthetic band. Blue-green 490 nm penetrates furthest in clear water, so this flatters clear water slightly. A3 The 1% level is a convention, not a biological boundary. Banse (2004) argued it should have been retired decades ago. The paper's headline depths use a different and much dimmer threshold (the light level that makes a Calanus copepod migrate, 0.027 microwatts per square metre at 490 nm), which is why the paper's photic depths are far deeper than any euphotic depth we compute here. A4 CONSEQUENCE OF A3: the shoaling areas in Part D come from the paper's Calanus-threshold photic depth. Our Beer-Lambert work in Parts A to C uses the 1% convention. These are two different depths. We do not mix them inside a single calculation, and the volume result inherits the paper's definition, not ours. A5 Every pixel inside a reported band is charged the same representative shoaling. Real distributions inside a band are not uniform. A6 Shoaling is treated as loss of a slab of water that used to be lit and now is not. Habitat is not really a slab and organisms are not evenly spread through it. This is a volume bookkeeping exercise. A7 Trends are treated as linear over 2003 to 2022 for the per-decade rate. That is 19 years of elapsed change, not 20. A8 Global ocean area 361.9 million km2, ocean volume 1.335 billion km3. ============================================================================= LIMITS ============================================================================= This script cannot tell you whether the satellite trend is real. It assumes the reported areas are correct and works out what they imply if they are. It has no error model for the satellite retrieval, no treatment of sensor calibration drift, and no way to separate a genuine twenty-year trend from the low-frequency tail of natural variability. Those are the arguments that decide whether the result stands, and none of them are in here. Python 3.12. numpy used if importable, pure-Python fallback otherwise. All randomness seeded with 20250527 (the paper's publication date). """ import math import random SEED = 20250527 try: import numpy as _np HAVE_NUMPY = True except ImportError: # pragma: no cover _np = None HAVE_NUMPY = False LN100 = math.log(100.0) # --------------------------------------------------------------------------- # Constants: the paper's reported areas, and standard ocean geometry # --------------------------------------------------------------------------- # Davies & Smyth (2025), reported to the square kilometre. AREA_KD_UP = 75_341_181 # Kd(490) increased (21%) AREA_KD_DOWN = 37_269_515 # Kd(490) decreased (10%) AREA_SHOAL_10 = 68_402_842 # photic depth reduced by >10 m (19%) AREA_SHOAL_50 = 32_449_129 # photic depth reduced by >50 m (9%) AREA_SHOAL_100 = 9_392_219 # photic depth reduced by >100 m (2.6%) OCEAN_AREA_KM2 = 361_900_000.0 # km2 OCEAN_VOL_KM3 = 1_335_000_000.0 # km3 EPIPELAGIC_M = 200.0 # m, conventional top of the twilight zone YEAR0, YEAR1 = 2003, 2022 ELAPSED_YEARS = YEAR1 - YEAR0 # 19 SECCHI_CONST = 1.7 # Poole & Atkins (1929) SECCHI_ALTS = {"Devlin 2008": 1.4, "Holmes 1970": 1.44, "Poole & Atkins 1929": 1.7, "Murray & Markager 2011": 2.0} def rule(ch="-", n=74): return ch * n def z_euphotic(kd, frac=0.01): """Depth at which irradiance falls to `frac` of the surface value.""" return -math.log(frac) / kd def light_fraction(kd, z): """I(z)/I0 for a given Kd and depth.""" return math.exp(-kd * z) def main(): random.seed(SEED) if HAVE_NUMPY: rng = _np.random.default_rng(SEED) else: rng = None out = [] P = out.append P(rule("=")) P("OCEAN DARKENING -- CLUB LIGHT-ATTENUATION AND LOST-VOLUME MODEL") P("Anchor: Davies, T. W. & Smyth, T. (2025).") P(" Darkening of the Global Ocean.") P(" Global Change Biology 31(5), e70227.") P(" DOI: 10.1111/gcb.70227") P("This is the club's own simplified model, NOT the paper's analysis.") P("numpy: %s | seed: %d" % (("yes, v" + _np.__version__) if HAVE_NUMPY else "no (pure-Python fallback)", SEED)) P(rule("=")) P("") # ----------------------------------------------------------------------- P("PART 0 -- WHAT THE PAPER REPORTS, IN OUR UNITS") P(rule()) P(" Kd(490) increased over %14s km2 (%.2f%% of ocean)" % (f"{AREA_KD_UP:,}", 100 * AREA_KD_UP / OCEAN_AREA_KM2)) P(" Kd(490) decreased over %14s km2 (%.2f%% of ocean)" % (f"{AREA_KD_DOWN:,}", 100 * AREA_KD_DOWN / OCEAN_AREA_KM2)) P(" net area darkened minus lit %14s km2" % f"{AREA_KD_UP - AREA_KD_DOWN:,}") P(" ratio darkened : brightened %14.2f : 1" % (AREA_KD_UP / AREA_KD_DOWN)) P("") P(" photic depth shoaled >10 m %14s km2 (%.2f%%)" % (f"{AREA_SHOAL_10:,}", 100 * AREA_SHOAL_10 / OCEAN_AREA_KM2)) P(" photic depth shoaled >50 m %14s km2 (%.2f%%)" % (f"{AREA_SHOAL_50:,}", 100 * AREA_SHOAL_50 / OCEAN_AREA_KM2)) P(" photic depth shoaled >100 m %14s km2 (%.2f%%)" % (f"{AREA_SHOAL_100:,}", 100 * AREA_SHOAL_100 / OCEAN_AREA_KM2)) P("") P(" For scale: the darkened area is %.2f times the land area of Russia" % (AREA_KD_UP / 17_098_246)) P(" and %.2f times the total land area of Earth (148.9 Mkm2)." % (AREA_KD_UP / 148_900_000)) P("") # ----------------------------------------------------------------------- P("PART A -- BEER-LAMBERT: HOW DEEP IS THE 1% LIGHT LEVEL?") P(rule()) P(" I(z) = I0 * exp(-Kd * z) -> z_1%% = ln(100)/Kd = %.5f/Kd" % LN100) P("") P(" Kd(490) water type z_1% z_10% z_0.1%") P(" (1/m) (m) (m) (m)") P(" " + rule("-", 66)) kd_ladder = [ (0.020, "clearest subtropical gyre"), (0.030, "clear oligotrophic open ocean"), (0.040, "open ocean, low chlorophyll"), (0.060, "open ocean, global average-ish"), (0.080, "productive open ocean"), (0.120, "shelf edge"), (0.200, "coastal, moderate"), (0.350, "coastal bloom"), (0.600, "turbid nearshore"), (1.000, "river plume"), ] for kd, label in kd_ladder: P(" %7.3f %-28s %7.1f %7.1f %8.1f" % (kd, label, z_euphotic(kd, 0.01), z_euphotic(kd, 0.10), z_euphotic(kd, 0.001))) P("") P(" Light remaining at fixed depths, as a percentage of the surface:") P("") depths = [10, 25, 50, 75, 100, 150, 200] P(" Kd(490) " + "".join("%9s" % ("%d m" % d) for d in depths)) P(" " + rule("-", 9 + 9 * len(depths))) for kd, _label in kd_ladder: row = "".join("%8.3f%%" % (100 * light_fraction(kd, d)) for d in depths) P(" %7.3f %s" % (kd, row)) P("") # ----------------------------------------------------------------------- P("PART B -- THE NONLINEARITY: METRES LOST PER UNIT OF ADDED MURK") P(rule()) P(" dz/dKd = -ln(100)/Kd^2. The cost of added murk scales as 1/Kd^2.") P("") P(" Kd start z_1% start +0.005 1/m +0.010 1/m +10% relative") P(" (1/m) (m) lost (m) lost (m) lost (m)") P(" " + rule("-", 68)) kd_probe = [0.020, 0.030, 0.040, 0.060, 0.080, 0.120, 0.200, 0.350, 0.600, 1.000] for kd in kd_probe: z0 = z_euphotic(kd) d005 = z0 - z_euphotic(kd + 0.005) d010 = z0 - z_euphotic(kd + 0.010) drel = z0 - z_euphotic(kd * 1.10) P(" %7.3f %9.1f %10.2f %10.2f %10.2f" % (kd, z0, d005, d010, drel)) P("") ratio = ((z_euphotic(0.02) - z_euphotic(0.025)) / (z_euphotic(1.0) - z_euphotic(1.005))) P(" The same +0.005 1/m of added murk costs the clearest gyre") P(" %.1f metres and the river plume %.4f metres." % (z_euphotic(0.02) - z_euphotic(0.025), z_euphotic(1.0) - z_euphotic(1.005))) P(" Ratio: %.0f to 1." % ratio) P("") frac_loss = 100 * (1 - 1 / 1.10) P(" But in RELATIVE terms the nonlinearity vanishes completely.") P(" A 10%% rise in Kd removes exactly %.4f%% of the euphotic depth" % frac_loss) P(" at every single value of Kd. Check, across five orders of magnitude:") for kd in [0.001, 0.01, 0.1, 1.0, 10.0]: z0 = z_euphotic(kd) z1 = z_euphotic(kd * 1.10) P(" Kd = %8.3f z_1%% %10.3f m -> %10.3f m loss %.4f%%" % (kd, z0, z1, 100 * (z0 - z1) / z0)) P("") P(" So 'the ocean lost 9% of its lit depth' and 'the ocean lost 50 m") P(" of lit depth' are statements about very different water.") P("") # ----------------------------------------------------------------------- P("PART C -- SECCHI DISK LOOKUP TABLE (bring your own rope)") P(rule()) P(" Kd = %.1f / Z_secchi -> z_1%% = ln(100)*Z_secchi/%.1f = %.4f * Z_secchi" % (SECCHI_CONST, SECCHI_CONST, LN100 / SECCHI_CONST)) P("") P(" Secchi Kd(490) z_1% z_10% loses this many m") P(" depth (1/m) (m) (m) per 10% more murk") P(" " + rule("-", 62)) secchi_ladder = [0.5, 1, 2, 3, 5, 7, 10, 15, 20, 25, 30, 35, 40, 50] for zsd in secchi_ladder: kd = SECCHI_CONST / zsd z0 = z_euphotic(kd) P(" %6.1f m %8.4f %7.1f %7.1f %10.2f" % (zsd, kd, z0, z_euphotic(kd, 0.10), z0 - z_euphotic(kd * 1.10))) P("") P(" Rule of thumb worth memorising: euphotic depth is about %.2f times" % (LN100 / SECCHI_CONST)) P(" whatever your Secchi reading was. Round it to 2.7.") P("") P(" Sensitivity to the conversion constant (Z_secchi = 20 m):") for name, c in sorted(SECCHI_ALTS.items(), key=lambda kv: kv[1]): kd = c / 20.0 P(" constant %.2f (%-22s) -> Kd = %.4f 1/m, z_1%% = %6.1f m" % (c, name, kd, z_euphotic(kd))) spread_lo = z_euphotic(2.0 / 20.0) spread_hi = z_euphotic(1.4 / 20.0) P(" Spread across published constants: %.1f to %.1f m, a factor of %.2f." % (spread_lo, spread_hi, spread_hi / spread_lo)) P(" This is why a Secchi disk is a screening tool and not a sensor.") P("") # ----------------------------------------------------------------------- P("PART D -- LOST LIT VOLUME FROM THE PAPER'S REPORTED AREAS") P(rule()) P(" Un-nesting the reported bands into disjoint bands:") P("") band_a_area = AREA_SHOAL_10 - AREA_SHOAL_50 band_b_area = AREA_SHOAL_50 - AREA_SHOAL_100 band_c_area = AREA_SHOAL_100 bands = [ ("shoaled 10-50 m", band_a_area, 10.0, 50.0, 30.0), ("shoaled 50-100 m", band_b_area, 50.0, 100.0, 75.0), ("shoaled >100 m", band_c_area, 100.0, 150.0, 100.0), ] P(" band area (km2) rep. shoaling lost volume (km3)") P(" " + rule("-", 70)) total_vol = 0.0 for name, area, _lo, _hi, rep in bands: vol = area * (rep / 1000.0) # km2 * km total_vol += vol P(" %-18s %14s %8.1f m %16s" % (name, f"{area:,}", rep, f"{vol:,.0f}")) P(" " + rule("-", 70)) P(" %-18s %14s %16s" % ("TOTAL", f"{AREA_SHOAL_10:,}", f"{total_vol:,.0f}")) P("") P(" Central estimate of lit habitat lost, 2003 to 2022:") P(" %.0f km3 ( %.3f million km3 )" % (total_vol, total_vol / 1e6)) P("") P(" As a fraction of things:") P(" of total ocean volume (%.3e km3) %.4f %%" % (OCEAN_VOL_KM3, 100 * total_vol / OCEAN_VOL_KM3)) epi_vol = OCEAN_AREA_KM2 * (EPIPELAGIC_M / 1000.0) P(" of the epipelagic 0-200 m (%.3e km3) %.3f %%" % (epi_vol, 100 * total_vol / epi_vol)) P(" smeared evenly over the whole ocean: a layer %.2f m thick" % (1000.0 * total_vol / OCEAN_AREA_KM2)) P(" smeared over just the darkened 21%%: a layer %.2f m thick" % (1000.0 * total_vol / AREA_KD_UP)) P("") P(" Rate, treating the change as linear over %d elapsed years:" % ELAPSED_YEARS) P(" %.0f km3 per year" % (total_vol / ELAPSED_YEARS)) P(" %.0f km3 per decade" % (total_vol / ELAPSED_YEARS * 10)) P(" %.4f %% of ocean volume per decade" % (100 * total_vol / ELAPSED_YEARS * 10 / OCEAN_VOL_KM3)) P(" %.0f km3 per day" % (total_vol / (ELAPSED_YEARS * 365.25))) P("") lake_superior = 12_100.0 # km3 P(" Comparisons, because km3 means nothing to anybody:") P(" Lake Superior holds %s km3. We lost %.0f Lake Superiors." % (f"{lake_superior:,.0f}", total_vol / lake_superior)) P(" That is %.1f Lake Superiors of lit water per year." % (total_vol / ELAPSED_YEARS / lake_superior)) P(" Mediterranean Sea volume is about 3,750,000 km3;") P(" our estimate is %.2f Mediterraneans." % (total_vol / 3_750_000.0)) P("") # ----------------------------------------------------------------------- P("PART D2 -- MONTE CARLO ON THE REPRESENTATIVE DEPTHS") P(rule()) P(" The only free choice above is where inside each band to put the") P(" representative shoaling. Sample it uniformly, 200,000 seeded draws.") P(" The top band is open-ended; we cap the sampler at 150 m, which is a") P(" guess, and it is the largest single source of spread in the answer.") P("") N = 200_000 if HAVE_NUMPY: da = rng.uniform(10.0, 50.0, N) db = rng.uniform(50.0, 100.0, N) dc = rng.uniform(100.0, 150.0, N) vols = (band_a_area * da + band_b_area * db + band_c_area * dc) / 1000.0 vols.sort() def q(p): return float(vols[int(p * (N - 1))]) med = float(_np.median(vols)) mean = float(vols.mean()) else: # pragma: no cover vols = [] for _ in range(N): v = (band_a_area * random.uniform(10.0, 50.0) + band_b_area * random.uniform(50.0, 100.0) + band_c_area * random.uniform(100.0, 150.0)) / 1000.0 vols.append(v) vols.sort() def q(p): return vols[int(p * (N - 1))] med = vols[N // 2] mean = sum(vols) / N P(" draws %s" % f"{N:,}") P(" mean %16s km3" % f"{mean:,.0f}") P(" median %16s km3" % f"{med:,.0f}") P(" 68%% interval %16s to %s km3" % (f"{q(0.16):,.0f}", f"{q(0.84):,.0f}")) P(" 95%% interval %16s to %s km3" % (f"{q(0.025):,.0f}", f"{q(0.975):,.0f}")) P(" our central estimate sits at the %.1fth percentile" % (100.0 * sum(1 for v in (vols if not HAVE_NUMPY else vols) if v <= total_vol) / N)) P("") P(" median as %% of ocean volume %.4f %%" % (100 * med / OCEAN_VOL_KM3)) P(" median per decade %s km3" % f"{med / ELAPSED_YEARS * 10:,.0f}") P("") # ----------------------------------------------------------------------- P("PART E -- WHAT Kd CHANGE WOULD PRODUCE A 50 m SHOALING?") P(rule()) P(" Inverting Beer-Lambert at the 1% level. If a patch of water starts") P(" at Kd0 and loses exactly 50 m of euphotic depth, the new Kd is") P(" Kd1 = ln(100) / ( ln(100)/Kd0 - 50 )") P(" and this is only possible at all if z_1%(Kd0) > 50 m,") P(" i.e. Kd0 < %.5f 1/m." % (LN100 / 50.0)) P("") P(" Kd0 z_1% (m) Kd1 needed change in Kd relative change") P(" " + rule("-", 68)) for kd0 in [0.020, 0.025, 0.030, 0.035, 0.040, 0.050, 0.060, 0.070, 0.080]: z0 = z_euphotic(kd0) if z0 <= 50.0 + 1e-9: P(" %6.3f %9.1f impossible -- already shallower than 50 m" % (kd0, z0)) continue kd1 = LN100 / (z0 - 50.0) P(" %6.3f %9.1f %9.4f %+10.4f %+8.1f %%" % (kd0, z0, kd1, kd1 - kd0, 100 * (kd1 - kd0) / kd0)) P("") P(" Read that column of percentages. In the clearest water a 50 m loss") P(" needs a Kd rise of about %.0f%%. In water at Kd = 0.08 it needs" % (100 * (LN100 / (z_euphotic(0.02) - 50.0) - 0.02) / 0.02)) P(" %.0f%%, which no satellite would call a subtle trend." % (100 * (LN100 / (z_euphotic(0.08) - 50.0) - 0.08) / 0.08)) P(" NOTE: the paper's >50 m shoaling figures are computed against the") P(" Calanus light threshold, which sits far below the 1% level, so the") P(" Kd changes THEY imply are smaller than this table. This table is a") P(" 1%-convention exercise and should be read as such.") P("") # ----------------------------------------------------------------------- P(rule("=")) P("HEADLINE NUMBERS USED IN THE ARTICLE") P(rule("=")) P(" euphotic depth = 2.709 x Secchi depth") P(" clearest gyre loses %.1f m per +0.005 1/m; river plume loses %.4f m" % (z_euphotic(0.02) - z_euphotic(0.025), z_euphotic(1.0) - z_euphotic(1.005))) P(" ratio of those two = %.0f : 1" % ratio) P(" a 10%% rise in Kd always costs %.2f%% of euphotic depth" % frac_loss) P(" lit volume lost (central) = %s km3" % f"{total_vol:,.0f}") P(" lit volume lost (MC median) = %s km3" % f"{med:,.0f}") P(" 68%% interval = %s to %s km3" % (f"{q(0.16):,.0f}", f"{q(0.84):,.0f}")) P(" = %.4f%% of ocean volume, %.2f%% of the epipelagic" % (100 * total_vol / OCEAN_VOL_KM3, 100 * total_vol / epi_vol)) P(" = a layer %.2f m thick over the entire ocean" % (1000.0 * total_vol / OCEAN_AREA_KM2)) P(" = %s km3 per decade" % f"{total_vol / ELAPSED_YEARS * 10:,.0f}") P(" = %.0f Lake Superiors, or %.1f per year" % (total_vol / lake_superior, total_vol / ELAPSED_YEARS / lake_superior)) P(rule("=")) text = "\n".join(out) print(text) return text if __name__ == "__main__": main()