Equations

Mathematical Reference

This page presents the mathematical formulation of the iCOW model from Ceres, Forest, and Keller (2019). All equations are implemented in ICOW.Core.

TipHow to Read This Page

This is a reference for the model’s equations, not a tutorial. If you are new to iCOW, start with the Quickstart for a hands-on introduction. Each section below includes a brief note on the physical reasoning behind the equation, followed by the precise formula.

Decision Levers

The model has five decision levers, stored in the FloodDefenses struct:

Lever	Symbol	Description	Units	Constraints
Withdrawal	\(W\)	Height below which city is relocated (absolute)	m	\(0 \leq W < H_{\text{city}}\)
Resistance Height	\(R\)	Height of flood-proofing above \(W\) (relative)	m	\(R \geq 0\)
Resistance Percentage	\(P\)	Fraction of buildings made resistant	–	\(0 \leq P < 1\)
Dike Height	\(D\)	Height of dike above its base (relative)	m	\(D \geq 0\)
Dike Base	\(B\)	Height of dike base above \(W\) (relative)	m	\(B \geq 0\)

\(W\) is the only absolute lever (measured from sea level). All other levers (\(R\), \(B\), \(D\)) are relative to \(W\). The absolute elevation of the dike base is \(W + B\); the dike top is at \(W + B + D\).

The constraint \(W + B + D \leq H_{\text{city}}\) ensures that defenses do not exceed the city’s maximum elevation.

Zone Definitions

The city is partitioned into five zones based on lever settings:

Zone	Elevation Range	Width	Description	Value Ratio
0	\(0\) to \(W\)	\(W\)	Withdrawn (no value)	0
1	\(W\) to \(W + \min(R, B)\)	\(\min(R, B)\)	Resistant	\(r_{\text{unprot}} = 0.95\)
2	\(W + \min(R, B)\) to \(W + B\)	\(B - R\) (if \(R < B\))	Unprotected gap	\(r_{\text{unprot}} = 0.95\)
3	\(W + B\) to \(W + B + D\)	\(D\)	Dike protected	\(r_{\text{prot}} = 1.1\)
4	\(W + B + D\) to \(H_{\text{city}}\)	\(H_{\text{city}} - W - B - D\)	Above dike	\(1.0\)

Zone 2 only exists when \(R < B\) (resistance does not reach the dike base).

Resistance-only strategies (\(B = 0\), \(D = 0\)): When there is no dike, Zone 1 extends to \(W + R\) (full resistance height), and Zones 2–3 have zero width. Zone 4 covers the remainder from \(W + R\) to \(H_{\text{city}}\).

NotePhysical Intuition: Why Five Zones?

The zones reflect how different parts of the city experience flooding differently. Zone 0 has been vacated entirely (no damage possible). Zones 1–2 are below the dike – they flood whenever surge exceeds the seawall, but resistant buildings in Zone 1 sustain less damage. Zone 3 is behind the dike – it floods only if the dike is overtopped or fails. Zone 4 is above the dike crest and floods only in extreme events. This zoning is what allows iCOW to capture graduated damage rather than a simple flood/no-flood binary.

Cross-Section Diagram

using ICOW
using CairoMakie
CairoMakie.activate!(; type="svg")

# Example defenses: W=2, R=1.5, P=0.5, D=4, B=3
fd = FloodDefenses(W=2.0, R=1.5, P=0.5, D=4.0, B=3.0)
H_city = 17.0

W, R, B, D = fd.W, fd.R, fd.B, fd.D
z1_top = W + min(R, B)
z2_top = W + B
z3_top = W + B + D

colors = (
    z0="#CC0099",  # dark magenta
    z1="#008B8B",  # dark cyan
    z2="#0000CD",  # medium blue
    z3="#FF8C00",  # dark orange
    z4="#228B22",  # forest green
)

fig = Figure(size=(700, 250))
ax = Axis(fig[1, 1],
    xlabel="Elevation (m)",
    limits=((-0.5, H_city + 0.5), (-0.4, 1)),
    xticks=0:2:H_city)
hideydecorations!(ax)
hidespines!(ax, :t, :l, :r)

bar_y = 0.5
bar_h = 0.35

function zone_bar!(ax, x_lo, x_hi, color, label)
    poly!(ax, Point2f[(x_lo, bar_y - bar_h / 2), (x_hi, bar_y - bar_h / 2),
            (x_hi, bar_y + bar_h / 2), (x_lo, bar_y + bar_h / 2)],
        color=color, strokecolor=:white, strokewidth=1)
    text!(ax, (x_lo + x_hi) / 2, bar_y, text=label, color=:white, fontsize=9,
        align=(:center, :center), font=:bold)
end

zone_bar!(ax, 0, W, colors.z0, "Zone 0\nWithdrawn")
zone_bar!(ax, W, z1_top, colors.z1, "Zone 1\nResistant")
if R < B
    zone_bar!(ax, z1_top, z2_top, colors.z2, "Zone 2\nGap")
end
zone_bar!(ax, z2_top, z3_top, colors.z3, "Zone 3\nDike-protected")
zone_bar!(ax, z3_top, H_city, colors.z4, "Zone 4\nAbove-dike")

function bracket!(ax, x_lo, x_hi, label, y_pos)
    lines!(ax, [x_lo, x_lo], [y_pos, y_pos - 0.08], color=:gray30, linewidth=1)
    lines!(ax, [x_hi, x_hi], [y_pos, y_pos - 0.08], color=:gray30, linewidth=1)
    lines!(ax, [x_lo, x_hi], [y_pos - 0.08, y_pos - 0.08], color=:gray30, linewidth=1)
    text!(ax, (x_lo + x_hi) / 2, y_pos - 0.12, text=label, fontsize=9,
        align=(:center, :top), color=:gray30)
end

y_base = bar_y - bar_h / 2 - 0.02
bracket!(ax, 0, W, "W", y_base)
bracket!(ax, W, W + R, "R", y_base)
bracket!(ax, W, z2_top, "B", y_base - 0.18)
bracket!(ax, z2_top, z3_top, "D", y_base)

fig

Zone Value Ratios

The value ratios \(r_{\text{prot}} = 1.1\) and \(r_{\text{unprot}} = 0.95\) mean that zone values do not sum exactly to \(V_w\). This is intentional:

• Protected areas appreciate: Dike protection increases land/building values
• Unprotected areas depreciate: Residual flood risk reduces property values

Cost Equations

Withdrawal Cost (Equation 1)

Relocating people and buildings from the lowest-lying areas. The cost increases nonlinearly with \(W\) because the denominator \(H_{\text{city}} - W\) shrinks: as more of the city is vacated, the remaining land becomes more congested and valuable, making each additional meter of withdrawal more expensive.

\[ C_W = \frac{V_{\text{city}} \cdot W \cdot f_w}{H_{\text{city}} - W} \]

City Value After Withdrawal (Equation 2)

Withdrawal permanently removes some city value. The loss is small (\(f_l = 0.01\) per meter) because most value migrates to higher ground rather than being destroyed.

\[ V_w = V_{\text{city}} \cdot \left(1 - \frac{f_l \cdot W}{H_{\text{city}}}\right) \]

Resistance Cost Fraction (Equation 3)

The cost of flood-proofing a fraction \(P\) of buildings. The linear term (\(f_{\text{lin}} \cdot P\)) captures the base cost. The exponential term activates above a threshold \(t_{\text{exp}} = 0.4\) and grows sharply as \(P \to 1\), reflecting the practical difficulty of retrofitting the last holdout buildings (historic structures, buildings with complex foundations, etc.).

\[ f_{\text{cR}} = f_{\text{adj}} \cdot \left( f_{\text{lin}} \cdot P + \frac{f_{\text{exp}} \cdot \max(0, P - t_{\text{exp}})}{1 - P} \right) \]

The factor \(f_{\text{adj}} = 1.25\) is present in the original implementation but not prominently shown in the paper.

Resistance Cost (Equations 4–5)

The total resistance cost depends on how far up the buildings are protected (\(R\)) and what fraction of buildings are treated (\(P\), via \(f_{\text{cR}}\) above). The cost also includes basement depth \(b\) because flood-proofing must extend below ground level.

When \(R < B\) (unconstrained):

\[ C_R = \frac{V_w \cdot f_{\text{cR}} \cdot R \cdot (R/2 + b)}{H_{\text{bldg}} \cdot (H_{\text{city}} - W)} \]

When \(R \geq B\) (resistance capped at dike base):

\[ C_R = \frac{V_w \cdot f_{\text{cR}} \cdot B \cdot (R - B/2 + b)}{H_{\text{bldg}} \cdot (H_{\text{city}} - W)} \]

Note

\(R > B\) is a dominated strategy: physical protection is capped at \(B\) but cost increases with \(R\). The optimizer will naturally discover this.

Dike Volume (Equation 6)

The dike is a three-dimensional structure: a trapezoidal seawall along the coastline with two tapered wings running inland up the slope. The wings taper because the terrain rises – less dike height is needed further inland.

\[ V_{\text{dike}} = \underbrace{W_{\text{city}} \cdot h_d \left( w_d + \frac{h_d}{s^2} \right)}_{\text{Main Seawall}} + \underbrace{\frac{h_d^2}{S} \left( w_d + \frac{2 h_d}{3 s^2} \right)}_{\text{Side Wings}} \]

Where \(h_d = D + D_{\text{startup}}\) is the effective dike height and \(S = H_{\text{city}} / D_{\text{city}}\) is the terrain slope.

Dike Cost (Equation 7)

\[ C_D = V_{\text{dike}} \cdot c_d \]

Total Investment

\[ C_{\text{total}} = C_W + C_R + C_D \]

Damage Equations

Effective Surge Height

Raw surge heights are adjusted for wave runup (\(f_{\text{runup}} = 1.1\)) and reduced by the existing seawall. Surges below the seawall cause no flooding.

\[ h_{\text{eff}} = \begin{cases} 0 & \text{if } h_{\text{raw}} \leq H_{\text{seawall}} \\ h_{\text{raw}} \cdot f_{\text{runup}} - H_{\text{seawall}} & \text{if } h_{\text{raw}} > H_{\text{seawall}} \end{cases} \]

Dike Failure Probability (Equation 8)

The dike does not simply “hold or fail” at a fixed height. Instead, failure probability ramps linearly from a baseline \(p_{\text{min}} = 0.05\) (accounting for seepage, structural defects, etc.) to certainty as the surge approaches the dike crest. The threshold \(t_{\text{fail}} = 0.95\) means the ramp begins when the surge reaches 95% of the dike height.

\[ p_{\text{fail}} = \begin{cases} p_{\text{min}} & \text{if } h_{\text{surge}} < t_{\text{fail}} \cdot D \\ \frac{h_{\text{surge}} - t_{\text{fail}} \cdot D}{D(1 - t_{\text{fail}})} & \text{if } t_{\text{fail}} \cdot D \leq h_{\text{surge}} < D \\ 1.0 & \text{if } h_{\text{surge}} \geq D \end{cases} \]

Here \(h_{\text{surge}}\) is the surge height above the dike base (not absolute elevation).

Zone Damage (Equation 9)

Base damage for zone \(Z\):

\[ d_Z = Val_Z \cdot \frac{Vol_F}{Vol_Z} \cdot f_{\text{damage}} \]

The fraction \(Vol_F / Vol_Z\) captures how much of the zone is flooded (by volume, accounting for the wedge geometry). The factor \(f_{\text{damage}} = 0.39\) is the maximum fraction of a zone’s value destroyed when fully inundated.

Zone-specific modifiers:

• Zone 1 (resistant): \(d_1 = \text{base} \cdot (1 - P)\) – only non-resistant buildings are damaged
• Zone 3 (dike-protected): \(d_3 = \text{base} \cdot f_{\text{dike}}\) where \(f_{\text{dike}} = f_{\text{intact}}\) if dike holds, \(f_{\text{failed}}\) if dike fails
• Zones 2, 4: standard damage (no modifier)

NoteWhy \(f_{\text{failed}} = 1.5 > 1\)?

When a dike fails, damage behind it exceeds what would occur without a dike. This is because the sudden rush of water through a breach causes more destruction than gradual inundation – and because residents behind the dike may not have evacuated or prepared, having trusted the dike to hold.

Zone Values

Standard case (with dike):

\[ Val_{Z1} = V_w \cdot r_{\text{unprot}} \cdot \frac{\min(R, B)}{H_{\text{city}} - W} \]

\[ Val_{Z2} = V_w \cdot r_{\text{unprot}} \cdot \frac{\max(0, B - R)}{H_{\text{city}} - W} \]

\[ Val_{Z3} = V_w \cdot r_{\text{prot}} \cdot \frac{D}{H_{\text{city}} - W} \]

\[ Val_{Z4} = V_w \cdot \frac{H_{\text{city}} - W - B - D}{H_{\text{city}} - W} \]

Resistance-only case (\(B = 0\), \(D = 0\)): When there is no dike, Zone 1 uses the full resistance height \(R\) without the \(r_{\text{unprot}}\) multiplier:

\[ Val_{Z1} = V_w \cdot \frac{R}{H_{\text{city}} - W}, \quad Val_{Z2} = Val_{Z3} = 0, \quad Val_{Z4} = V_w \cdot \frac{H_{\text{city}} - W - R}{H_{\text{city}} - W} \]

Threshold Damage

When total damage exceeds a threshold \(d_{\text{thresh}}\), cascading societal effects amplify losses. This captures phenomena like insurance market collapse, mass displacement, and compound vulnerabilities that make extreme events disproportionately destructive.

\[ d_{\text{total}} = \begin{cases} \sum d_Z & \text{if } \sum d_Z \leq d_{\text{thresh}} \\ \sum d_Z + \left( f_{\text{thresh}} \cdot (\sum d_Z - d_{\text{thresh}}) \right)^{\gamma_{\text{thresh}}} & \text{if } \sum d_Z > d_{\text{thresh}} \end{cases} \]

Expected Annual Damage Integration

EAD mode uses two-level integration:

\[ \text{EAD} = \int p_h(h) \cdot \mathbb{E}[\text{damage} \mid h] \, dh \]

The inner expectation (analytical) integrates over dike failure:

\[ \mathbb{E}[\text{damage} \mid h] = p_{\text{fail}}(h) \cdot d_{\text{failed}}(h) + (1 - p_{\text{fail}}(h)) \cdot d_{\text{intact}}(h) \]

The outer expectation is computed numerically via adaptive quadrature or Monte Carlo sampling.

Parameters

All 28 parameters are fields in StochasticConfig / EADConfig. Default values are read directly from the struct definition:

using ICOW.EAD

config = EADConfig()

# Metadata: (symbol, category, units, description)
meta = Dict(
    :V_city       => (raw"$V_{\text{city}}$",       "Geometry",   "USD",       "Total replacement value of all city infrastructure and buildings"),
    :H_bldg       => (raw"$H_{\text{bldg}}$",       "Geometry",   "m",         "Average building height, used to normalize resistance costs per floor"),
    :H_city       => (raw"$H_{\text{city}}$",        "Geometry",   "m",         "Total elevation change from the seawall to the highest point of the city"),
    :D_city       => (raw"$D_{\text{city}}$",        "Geometry",   "m",         "Horizontal distance from the seawall to the city's peak elevation"),
    :W_city       => (raw"$W_{\text{city}}$",        "Geometry",   "m",         "Length of the coastline along which the seawall and dike are built"),
    :H_seawall    => (raw"$H_{\text{seawall}}$",     "Geometry",   "m",         "Height of the existing seawall; surges below this cause no flooding"),
    :D_startup    => (raw"$D_{\text{startup}}$",     "Dike",       "m",         "Minimum effective dike height representing fixed mobilization costs (added to D)"),
    :w_d          => (raw"$w_d$",                     "Dike",       "m",         "Width of the flat top of the dike cross-section"),
    :s_dike       => (raw"$s$",                       "Dike",       "m/m",       raw"Dike side slope parameter (horizontal run per unit rise = $1/s^2$)"),
    :c_d          => (raw"$c_d$",                     "Dike",       "USD/m³",    "Unit cost of dike construction per cubic meter"),
    :r_prot       => (raw"$r_{\text{prot}}$",        "Zones",      "--",        "Value multiplier for dike-protected land (Zone 3)"),
    :r_unprot     => (raw"$r_{\text{unprot}}$",      "Zones",      "--",        "Value multiplier for unprotected land (Zones 1--2)"),
    :f_w          => (raw"$f_w$",                     "Withdrawal", "--",        "Scales the cost of relocating buildings and infrastructure from the withdrawn zone"),
    :f_l          => (raw"$f_l$",                     "Withdrawal", "--",        "Fraction of city value permanently lost per meter of withdrawal height"),
    :f_adj        => (raw"$f_{\text{adj}}$",         "Resistance", "--",        "Overall multiplier on resistance costs (present in C++ but not prominent in paper)"),
    :f_lin        => (raw"$f_{\text{lin}}$",         "Resistance", "--",        raw"Weight of the linear component of resistance cost as a function of $P$"),
    :f_exp        => (raw"$f_{\text{exp}}$",         "Resistance", "--",        raw"Weight of the exponential component that makes high $P$ values increasingly expensive"),
    :t_exp        => (raw"$t_{\text{exp}}$",         "Resistance", "--",        raw"Threshold of $P$ above which the exponential cost component activates"),
    :b_basement   => (raw"$b$",                       "Resistance", "m",         "Depth of basements below ground; flood-proofing must extend this far below the surface"),
    :f_damage     => (raw"$f_{\text{damage}}$",      "Damage",     "--",        "Maximum fraction of a zone's value destroyed when fully inundated"),
    :f_intact     => (raw"$f_{\text{intact}}$",      "Damage",     "--",        "Damage multiplier for Zone 3 when the dike holds (seepage and wave overtopping)"),
    :f_failed     => (raw"$f_{\text{failed}}$",      "Damage",     "--",        "Damage multiplier for Zone 3 when the dike fails (amplified by sudden breach flooding)"),
    :t_fail       => (raw"$t_{\text{fail}}$",        "Damage",     "--",        raw"Fraction of dike height $D$ at which failure probability begins to rise above $p_{\text{min}}$"),
    :p_min        => (raw"$p_{\text{min}}$",         "Damage",     "--",        "Baseline probability of dike failure even when surge is well below the dike crest"),
    :f_runup      => (raw"$f_{\text{runup}}$",       "Damage",     "--",        "Multiplier on raw surge height to account for wave runup effects at the seawall"),
    :d_thresh     => (raw"$d_{\text{thresh}}$",      "Threshold",  "USD",       "Damage level above which cascading societal effects amplify losses"),
    :f_thresh     => (raw"$f_{\text{thresh}}$",      "Threshold",  "--",        raw"Scales the excess damage above $d_{\text{thresh}}$ before applying the power-law penalty"),
    :gamma_thresh => (raw"$\gamma_{\text{thresh}}$", "Threshold",  "--",        "Exponent of the power-law penalty applied to excess damage above the threshold"),
)

# Emit a Markdown table so Quarto renders the math
header = "| Category | Symbol | Field | Default | Units | Description |"
sep    = "|----------|--------|-------|---------|-------|-------------|"
println(header)
println(sep)
for fn in fieldnames(EADConfig)
    haskey(meta, fn) || continue
    sym, cat, units, desc = meta[fn]
    val = getfield(config, fn)
    println("| $cat | $sym | `$fn` | $val | $units | $desc |")
end
println()
println(": {tbl-colwidths=\"[10,12,10,8,6,54]\"}")

Category	Symbol	Field	Default	Units	Description
Geometry	\(V_{\text{city}}\)	V_city	1.5e12	USD	Total replacement value of all city infrastructure and buildings
Geometry	\(H_{\text{bldg}}\)	H_bldg	30.0	m	Average building height, used to normalize resistance costs per floor
Geometry	\(H_{\text{city}}\)	H_city	17.0	m	Total elevation change from the seawall to the highest point of the city
Geometry	\(D_{\text{city}}\)	D_city	2000.0	m	Horizontal distance from the seawall to the city’s peak elevation
Geometry	\(W_{\text{city}}\)	W_city	43000.0	m	Length of the coastline along which the seawall and dike are built
Geometry	\(H_{\text{seawall}}\)	H_seawall	1.75	m	Height of the existing seawall; surges below this cause no flooding
Dike	\(D_{\text{startup}}\)	D_startup	2.0	m	Minimum effective dike height representing fixed mobilization costs (added to D)
Dike	\(w_d\)	w_d	3.0	m	Width of the flat top of the dike cross-section
Dike	\(s\)	s_dike	0.5	m/m	Dike side slope parameter (horizontal run per unit rise = \(1/s^2\))
Dike	\(c_d\)	c_d	1000.0	USD/m³	Unit cost of dike construction per cubic meter
Zones	\(r_{\text{prot}}\)	r_prot	1.1	–	Value multiplier for dike-protected land (Zone 3)
Zones	\(r_{\text{unprot}}\)	r_unprot	0.95	–	Value multiplier for unprotected land (Zones 1–2)
Withdrawal	\(f_w\)	f_w	1.0	–	Scales the cost of relocating buildings and infrastructure from the withdrawn zone
Withdrawal	\(f_l\)	f_l	0.01	–	Fraction of city value permanently lost per meter of withdrawal height
Resistance	\(f_{\text{adj}}\)	f_adj	1.25	–	Overall multiplier on resistance costs (present in C++ but not prominent in paper)
Resistance	\(f_{\text{lin}}\)	f_lin	0.35	–	Weight of the linear component of resistance cost as a function of \(P\)
Resistance	\(f_{\text{exp}}\)	f_exp	0.115	–	Weight of the exponential component that makes high \(P\) values increasingly expensive
Resistance	\(t_{\text{exp}}\)	t_exp	0.4	–	Threshold of \(P\) above which the exponential cost component activates
Resistance	\(b\)	b_basement	3.0	m	Depth of basements below ground; flood-proofing must extend this far below the surface
Damage	\(f_{\text{damage}}\)	f_damage	0.39	–	Maximum fraction of a zone’s value destroyed when fully inundated
Damage	\(f_{\text{intact}}\)	f_intact	0.03	–	Damage multiplier for Zone 3 when the dike holds (seepage and wave overtopping)
Damage	\(f_{\text{failed}}\)	f_failed	1.5	–	Damage multiplier for Zone 3 when the dike fails (amplified by sudden breach flooding)
Damage	\(t_{\text{fail}}\)	t_fail	0.95	–	Fraction of dike height \(D\) at which failure probability begins to rise above \(p_{\text{min}}\)
Damage	\(p_{\text{min}}\)	p_min	0.05	–	Baseline probability of dike failure even when surge is well below the dike crest
Damage	\(f_{\text{runup}}\)	f_runup	1.1	–	Multiplier on raw surge height to account for wave runup effects at the seawall
Threshold	\(d_{\text{thresh}}\)	d_thresh	4.0e9	USD	Damage level above which cascading societal effects amplify losses
Threshold	\(f_{\text{thresh}}\)	f_thresh	1.0	–	Scales the excess damage above \(d_{\text{thresh}}\) before applying the power-law penalty
Threshold	\(\gamma_{\text{thresh}}\)	gamma_thresh	1.01	–	Exponent of the power-law penalty applied to excess damage above the threshold

Deviations from Ceres, Forest, and Keller (2019)

This implementation fixes several bugs in the original C++ reference code and makes one parameter recalibration. Full details, including line-by-line C++ bug analysis, are in _background/equations.md in the source repository (not published as part of these docs).

Dike volume formula (Equation 6)

The paper’s Equation 6 contains a polynomial \(T\) under a square root that is numerically unstable for realistic terrain slopes. With \(S = H_{\text{city}} / D_{\text{city}} = 0.0085\), the polynomial denominators involve \(S^4 \approx 5 \times 10^{-9}\), causing \(T < 0\) and the square root to be undefined.

The original C++ code never actually evaluated this formula correctly due to Bug #1: pow(T, 1/2) performs integer division in C++, yielding pow(T, 0) = 1. The complex geometric term was silently ignored in all published results.

This implementation replaces Equation 6 with a simplified geometric derivation that is numerically stable and physically equivalent: a trapezoidal seawall prism along the coastline plus two tapered wing prisms running inland up the city slope.

Dike unit cost (\(c_d\))

The paper and C++ both specify \(c_d = 10\) dollars per cubic meter. However, the C++ CalculateDikeCost function contains Bug #4: it uses CityWidth (43,000 m, the coastline length) where the formula calls for WidthDikeTop (3 m, the dike crest width). This inflates the wing volume term by a factor of $$14,000, which combined with the wrong slope definition (Bug #5) produces dike costs roughly \(100 \times\) larger than the correct geometric formula.

Since the paper’s published results were generated with these bugs, its tradeoff analysis effectively assumes dike costs of 10–20 billion dollars for a Manhattan-scale seawall. To restore that cost scale with the corrected volume formula, we set \(c_d = 1000\) dollars per cubic meter. This is also more physically realistic: 10 dollars per cubic meter covers only raw earth fill, while 1000 dollars per cubic meter better reflects the full cost of engineered flood defense construction.

Other C++ bug fixes

The original C++ reference code contains eight bugs total. In addition to the dike-related bugs above:

• Bug #2: Uses constant array index dh=5 instead of local variable ch in the fourth term of the \(T\) polynomial.
• Bug #3: Algebraic error in the fourth term: -4*ch2+ch2/sd^2 should be -3*ch2/sd^2.
• Bug #5: City slope computed as CityLength/CityWidth (= 0.047) instead of H_city/D_city (= 0.0085).
• Bug #6: Resistance cost (Equations 4–5) uses zone value \(Val_{Z1}\) instead of \(V_w\) (value after withdrawal).
• Bug #7: \(V_w\) computed as \(V_{\text{city}} - C_W\) instead of Equation 2.
• Bug #8: When \(B < 0.1\) m, the C++ code sets \(R = 0\), preventing “resistance-only” strategies (flood-proofing without a dike). The correct behavior is to use Equation 4 (unconstrained) when there is no dike (\(B = 0\) and \(D = 0\)).

This implementation fixes all eight bugs to match the paper’s intended formulas.

References

Ceres, Robert L., Chris E. Forest, and Klaus Keller. 2019. “Optimization of Multiple Storm Surge Risk Mitigation Strategies for an Island City On a Wedge.” Environmental Modelling & Software 119 (September): 341–53. https://doi.org/10.1016/j.envsoft.2019.06.011.