Quickstart

Your First iCOW Analysis

This tutorial walks through a complete flood risk analysis using ICOW.jl. By the end, you will have configured a coastal city, tried several protection strategies, and compared their investment-vs-damage tradeoffs.

What You’ll Need

Julia 1.10 or later with ICOW.jl installed (see Home for installation)
Basic familiarity with Julia syntax (variables, functions, arrays)
No prior knowledge of flood risk modeling is assumed

Setup

Code

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

The Setting

Imagine you are advising a coastal city on flood protection investments. The city sits on a wedge of land that rises from sea level to $H_{\text{city}} = 17$ meters. It has a short existing seawall, but major storm surges can overtop it and cause serious damage.

You need to answer: what combination of withdrawal, resistance, and dikes should the city invest in?

Step 1: Configure the City

An EADConfig holds the 28 physical and economic parameters that define the city. The defaults describe a Manhattan-like setting.

config = EADConfig()

EADConfig{Float64}(1.5e12, 30.0, 17.0, 2000.0, 43000.0, 1.75, 2.0, 3.0, 0.5, 1000.0, 1.1, 0.95, 1.0, 0.01, 1.25, 0.35, 0.115, 0.4, 3.0, 0.39, 0.03, 1.5, 0.95, 0.05, 1.1, 4.0e9, 1.0, 1.01, 100, QuadratureIntegrator{Float64}(1.0e-6))

The most important parameters to understand are:

Parameter	Value	What It Means
`V_city`	$1.5T	Total value of all buildings and infrastructure
`H_city`	17.0 m	Elevation from waterfront to highest point
`H_seawall`	1.75 m	Existing seawall height (surges below this cause no flooding)
`H_bldg`	30.0 m	Average building height (normalizes resistance costs)

You can customize any parameter:

# A smaller, lower-lying city
config_small = EADConfig(; V_city=5e11, H_city=10.0, H_seawall=1.0)

EADConfig{Float64}(5.0e11, 30.0, 10.0, 2000.0, 43000.0, 1.0, 2.0, 3.0, 0.5, 1000.0, 1.1, 0.95, 1.0, 0.01, 1.25, 0.35, 0.115, 0.4, 3.0, 0.39, 0.03, 1.5, 0.95, 0.05, 1.1, 4.0e9, 1.0, 1.01, 100, QuadratureIntegrator{Float64}(1.0e-6))

For this tutorial, we will use the defaults.

Step 2: Define the Hazard

Storm surge heights follow a Generalized Extreme Value (GEV) distribution, which is the natural statistical model for annual maxima. An EADScenario specifies the GEV parameters, a discount rate, and a mean sea level time series (which can represent sea level rise):

n_years = config.n_years
scenario = EADScenario(;
    surge_loc=0.935,      # GEV location (m): calibrated to The Battery, NYC
    surge_scale=0.260,    # GEV scale (m)
    surge_shape=0.030,    # GEV shape: slightly heavy-tailed
    discount_rate=0.02,   # 2% annual discount rate
    mean_sea_level=collect(range(0.0; step=0.01, length=n_years)),  # 1 cm/yr SLR
)

EADScenario{Float64}(SimOptDecisions.ContinuousParameter{Float64}(0.935, (-5.0, 30.0)), SimOptDecisions.ContinuousParameter{Float64}(0.26, (0.01, 20.0)), SimOptDecisions.ContinuousParameter{Float64}(0.03, (-1.0, 1.0)), SimOptDecisions.ContinuousParameter{Float64}(0.02, (0.0, 1.0)), SimOptDecisions.TimeSeriesParameter{Float64, Int64}([1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  91, 92, 93, 94, 95, 96, 97, 98, 99, 100], [0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09  …  0.9, 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99]))

GEV Shape Parameter

The shape parameter controls the tail behavior:

shape = 0: Gumbel distribution (light-tailed, exponential decay)
shape > 0: Frechet distribution (heavy-tailed, more extreme surges possible)
shape < 0: Weibull distribution (bounded upper tail)

For most coastal locations, the shape parameter is near zero or slightly positive.

Step 3: Try a Protection Strategy

A StaticPolicy specifies how to allocate protection using five fractions, all between 0 and 1. The “stick-breaking” parameterization guarantees that all physical constraints are satisfied automatically (see Architecture for details).

Let’s try a strategy that invests in a moderate dike with some resistance:

policy_mixed = StaticPolicy(;
    a_frac=0.3,   # use 30% of city elevation as protection budget
    w_frac=0.1,   # 10% of budget goes to withdrawal
    b_frac=0.2,   # 20% of remaining budget goes to dike base
    r_frac=0.1,   # resistance height = 10% of H_city
    P=0.5,        # flood-proof 50% of buildings
)

StaticPolicy{Float64}(SimOptDecisions.ContinuousParameter{Float64}(0.3, (0.0, 1.0)), SimOptDecisions.ContinuousParameter{Float64}(0.1, (0.0, 1.0)), SimOptDecisions.ContinuousParameter{Float64}(0.2, (0.0, 1.0)), SimOptDecisions.ContinuousParameter{Float64}(0.1, (0.0, 1.0)), SimOptDecisions.ContinuousParameter{Float64}(0.5, (0.0, 0.99)))

What does this translate to in meters? The policy converts to a FloodDefenses struct:

fd = FloodDefenses(policy_mixed, config)
fd

FloodDefenses(W=0.51m, R=1.7000000000000002m, P=50.0%, D=3.6719999999999997m, B=0.918m)

Step 4: Run the Simulation

The simulate function computes the total investment cost and expected flood damage over the planning horizon (100 years by default), discounted to present value:

rng = MersenneTwister(42)
outcome = simulate(config, scenario, policy_mixed, rng)
outcome

EADOutcome{Float64}(SimOptDecisions.ContinuousParameter{Float64}(5.455527272132463e10, (-Inf, Inf)), SimOptDecisions.ContinuousParameter{Float64}(2.1096777425625317e9, (-Inf, Inf)))

The outcome tells you the present-value investment cost and expected damage. The total cost is their sum:

total_cost(outcome)

5.666495046388716e10

Step 5: Compare Strategies

One strategy is not enough – you need to compare alternatives. Here are five archetypes ranging from “do nothing” to “heavy investment”:

strategies = [
    ("No additional protection", StaticPolicy(; a_frac=0.0, w_frac=0.0, b_frac=0.0, r_frac=0.0, P=0.0)),
    ("Dike only",                StaticPolicy(; a_frac=0.5, w_frac=0.0, b_frac=0.0, r_frac=0.0, P=0.0)),
    ("Withdraw + resist",        StaticPolicy(; a_frac=0.2, w_frac=0.5, b_frac=1.0, r_frac=0.3, P=0.7)),
    ("Mixed",                    StaticPolicy(; a_frac=0.3, w_frac=0.1, b_frac=0.2, r_frac=0.1, P=0.5)),
    ("High dike + resist",       StaticPolicy(; a_frac=0.5, w_frac=0.0, b_frac=0.1, r_frac=0.1, P=0.3)),
]
policy_list = last.(strategies)
strategy_names = first.(strategies)

5-element Vector{String}:
 "No additional protection"
 "Dike only"
 "Withdraw + resist"
 "Mixed"
 "High dike + resist"

The explore() function evaluates all policies against a scenario in one call and returns a structured dataset:

result = explore(config, [scenario], policy_list; progress=false)

YAXArray Dataset
Shared Axes: 
None
Variables with additional axes:
  Additional Axes: 
  (↓ policy Sampled{Int64} 1:5 ForwardOrdered Regular Points,
  → scenario Sampled{Int64} 1:1 ForwardOrdered Regular Points)
  Variables: 
  expected_damage, investment

  Additional Axes: 
  (↓ policy Sampled{Int64} 1:5 ForwardOrdered Regular Points)
  Variables: 
  policy_P, policy_a_frac, policy_b_frac, policy_r_frac, policy_w_frac

  Additional Axes: 
  (↓ scenario Sampled{Int64} 1:1 ForwardOrdered Regular Points)
  Variables: 
  scenario_discount_rate, scenario_mean_sea_level[100], scenario_mean_sea_level[10], scenario_mean_sea_level[11], scenario_mean_sea_level[12], scenario_mean_sea_level[13], scenario_mean_sea_level[14], scenario_mean_sea_level[15], scenario_mean_sea_level[16], scenario_mean_sea_level[17], scenario_mean_sea_level[18], scenario_mean_sea_level[19], scenario_mean_sea_level[1], scenario_mean_sea_level[20], scenario_mean_sea_level[21], scenario_mean_sea_level[22], scenario_mean_sea_level[23], scenario_mean_sea_level[24], scenario_mean_sea_level[25], scenario_mean_sea_level[26], scenario_mean_sea_level[27], scenario_mean_sea_level[28], scenario_mean_sea_level[29], scenario_mean_sea_level[2], scenario_mean_sea_level[30], scenario_mean_sea_level[31], scenario_mean_sea_level[32], scenario_mean_sea_level[33], scenario_mean_sea_level[34], scenario_mean_sea_level[35], scenario_mean_sea_level[36], scenario_mean_sea_level[37], scenario_mean_sea_level[38], scenario_mean_sea_level[39], scenario_mean_sea_level[3], scenario_mean_sea_level[40], scenario_mean_sea_level[41], scenario_mean_sea_level[42], scenario_mean_sea_level[43], scenario_mean_sea_level[44], scenario_mean_sea_level[45], scenario_mean_sea_level[46], scenario_mean_sea_level[47], scenario_mean_sea_level[48], scenario_mean_sea_level[49], scenario_mean_sea_level[4], scenario_mean_sea_level[50], scenario_mean_sea_level[51], scenario_mean_sea_level[52], scenario_mean_sea_level[53], scenario_mean_sea_level[54], scenario_mean_sea_level[55], scenario_mean_sea_level[56], scenario_mean_sea_level[57], scenario_mean_sea_level[58], scenario_mean_sea_level[59], scenario_mean_sea_level[5], scenario_mean_sea_level[60], scenario_mean_sea_level[61], scenario_mean_sea_level[62], scenario_mean_sea_level[63], scenario_mean_sea_level[64], scenario_mean_sea_level[65], scenario_mean_sea_level[66], scenario_mean_sea_level[67], scenario_mean_sea_level[68], scenario_mean_sea_level[69], scenario_mean_sea_level[6], scenario_mean_sea_level[70], scenario_mean_sea_level[71], scenario_mean_sea_level[72], scenario_mean_sea_level[73], scenario_mean_sea_level[74], scenario_mean_sea_level[75], scenario_mean_sea_level[76], scenario_mean_sea_level[77], scenario_mean_sea_level[78], scenario_mean_sea_level[79], scenario_mean_sea_level[7], scenario_mean_sea_level[80], scenario_mean_sea_level[81], scenario_mean_sea_level[82], scenario_mean_sea_level[83], scenario_mean_sea_level[84], scenario_mean_sea_level[85], scenario_mean_sea_level[86], scenario_mean_sea_level[87], scenario_mean_sea_level[88], scenario_mean_sea_level[89], scenario_mean_sea_level[8], scenario_mean_sea_level[90], scenario_mean_sea_level[91], scenario_mean_sea_level[92], scenario_mean_sea_level[93], scenario_mean_sea_level[94], scenario_mean_sea_level[95], scenario_mean_sea_level[96], scenario_mean_sea_level[97], scenario_mean_sea_level[98], scenario_mean_sea_level[99], scenario_mean_sea_level[9], scenario_surge_loc, scenario_surge_scale, scenario_surge_shape

Now plot the investment vs. expected damage tradeoff:

inv_vals = vec(result[:investment].data) ./ 1e9
dmg_vals = vec(result[:expected_damage].data) ./ 1e9

fig = Figure(; size=(650, 420))
ax = Axis(fig[1, 1];
    xlabel="Investment (\$B)",
    ylabel="Expected Damage (\$B)",
    title="Investment vs Expected Damage")
scatter!(ax, inv_vals, dmg_vals; markersize=14, color=:steelblue)
for (i, name) in enumerate(strategy_names)
    text!(ax, inv_vals[i], dmg_vals[i]; text=name, fontsize=10, offset=(8, 4))
end
fig

Reading the Plot

Points in the lower-left are cheapest overall (low investment + low damage). The ideal strategy minimizes the sum of investment and expected damage – but this depends on what you’re optimizing for. A risk-averse decision-maker might prefer a strategy with slightly higher expected cost but much lower tail risk.

Step 6: What if the Hazard Is Different?

So far we assumed a specific surge distribution. But in practice, the GEV parameters are uncertain – this is a key aspect of decision making under deep uncertainty.

Let’s see how our strategies perform under a heavier-tailed surge distribution (positive shape parameter):

scenario_heavy = EADScenario(;
    surge_loc=0.935,
    surge_scale=0.260,
    surge_shape=0.15,   # heavier tail: more extreme surges possible
    discount_rate=0.02,
    mean_sea_level=collect(range(0.0; step=0.01, length=n_years)),
)

scenarios_both = [scenario, scenario_heavy]
result_both = explore(config, scenarios_both, policy_list; progress=false)

YAXArray Dataset
Shared Axes: 
None
Variables with additional axes:
  Additional Axes: 
  (↓ policy Sampled{Int64} 1:5 ForwardOrdered Regular Points,
  → scenario Sampled{Int64} 1:2 ForwardOrdered Regular Points)
  Variables: 
  expected_damage, investment

  Additional Axes: 
  (↓ policy Sampled{Int64} 1:5 ForwardOrdered Regular Points)
  Variables: 
  policy_P, policy_a_frac, policy_b_frac, policy_r_frac, policy_w_frac

  Additional Axes: 
  (↓ scenario Sampled{Int64} 1:2 ForwardOrdered Regular Points)
  Variables: 
  scenario_discount_rate, scenario_mean_sea_level[100], scenario_mean_sea_level[10], scenario_mean_sea_level[11], scenario_mean_sea_level[12], scenario_mean_sea_level[13], scenario_mean_sea_level[14], scenario_mean_sea_level[15], scenario_mean_sea_level[16], scenario_mean_sea_level[17], scenario_mean_sea_level[18], scenario_mean_sea_level[19], scenario_mean_sea_level[1], scenario_mean_sea_level[20], scenario_mean_sea_level[21], scenario_mean_sea_level[22], scenario_mean_sea_level[23], scenario_mean_sea_level[24], scenario_mean_sea_level[25], scenario_mean_sea_level[26], scenario_mean_sea_level[27], scenario_mean_sea_level[28], scenario_mean_sea_level[29], scenario_mean_sea_level[2], scenario_mean_sea_level[30], scenario_mean_sea_level[31], scenario_mean_sea_level[32], scenario_mean_sea_level[33], scenario_mean_sea_level[34], scenario_mean_sea_level[35], scenario_mean_sea_level[36], scenario_mean_sea_level[37], scenario_mean_sea_level[38], scenario_mean_sea_level[39], scenario_mean_sea_level[3], scenario_mean_sea_level[40], scenario_mean_sea_level[41], scenario_mean_sea_level[42], scenario_mean_sea_level[43], scenario_mean_sea_level[44], scenario_mean_sea_level[45], scenario_mean_sea_level[46], scenario_mean_sea_level[47], scenario_mean_sea_level[48], scenario_mean_sea_level[49], scenario_mean_sea_level[4], scenario_mean_sea_level[50], scenario_mean_sea_level[51], scenario_mean_sea_level[52], scenario_mean_sea_level[53], scenario_mean_sea_level[54], scenario_mean_sea_level[55], scenario_mean_sea_level[56], scenario_mean_sea_level[57], scenario_mean_sea_level[58], scenario_mean_sea_level[59], scenario_mean_sea_level[5], scenario_mean_sea_level[60], scenario_mean_sea_level[61], scenario_mean_sea_level[62], scenario_mean_sea_level[63], scenario_mean_sea_level[64], scenario_mean_sea_level[65], scenario_mean_sea_level[66], scenario_mean_sea_level[67], scenario_mean_sea_level[68], scenario_mean_sea_level[69], scenario_mean_sea_level[6], scenario_mean_sea_level[70], scenario_mean_sea_level[71], scenario_mean_sea_level[72], scenario_mean_sea_level[73], scenario_mean_sea_level[74], scenario_mean_sea_level[75], scenario_mean_sea_level[76], scenario_mean_sea_level[77], scenario_mean_sea_level[78], scenario_mean_sea_level[79], scenario_mean_sea_level[7], scenario_mean_sea_level[80], scenario_mean_sea_level[81], scenario_mean_sea_level[82], scenario_mean_sea_level[83], scenario_mean_sea_level[84], scenario_mean_sea_level[85], scenario_mean_sea_level[86], scenario_mean_sea_level[87], scenario_mean_sea_level[88], scenario_mean_sea_level[89], scenario_mean_sea_level[8], scenario_mean_sea_level[90], scenario_mean_sea_level[91], scenario_mean_sea_level[92], scenario_mean_sea_level[93], scenario_mean_sea_level[94], scenario_mean_sea_level[95], scenario_mean_sea_level[96], scenario_mean_sea_level[97], scenario_mean_sea_level[98], scenario_mean_sea_level[99], scenario_mean_sea_level[9], scenario_surge_loc, scenario_surge_scale, scenario_surge_shape

fig = Figure(; size=(700, 420))
ax = Axis(fig[1, 1];
    xlabel="Investment (\$B)",
    ylabel="Expected Damage (\$B)",
    title="Strategy Performance Under Two Hazard Scenarios")

# Light-tailed scenario
inv_1 = result_both[:investment].data[:, 1] ./ 1e9
dmg_1 = result_both[:expected_damage].data[:, 1] ./ 1e9
scatter!(ax, inv_1, dmg_1; markersize=12, color=:steelblue, label="Baseline (shape=0.03)")

# Heavy-tailed scenario
inv_2 = result_both[:investment].data[:, 2] ./ 1e9
dmg_2 = result_both[:expected_damage].data[:, 2] ./ 1e9
scatter!(ax, inv_2, dmg_2; markersize=12, color=:firebrick, marker=:diamond, label="Frechet (shape=0.15)")

# Connect same policies across scenarios
for i in 1:length(policy_list)
    lines!(ax, [inv_1[i], inv_2[i]], [dmg_1[i], dmg_2[i]]; color=:gray70, linewidth=1)
    text!(ax, inv_2[i], dmg_2[i]; text=strategy_names[i], fontsize=9, offset=(8, 4))
end
axislegend(ax; position=:rt)
fig

Key Insight

Investment costs stay the same across scenarios (they depend only on the policy), but expected damage shifts substantially. Strategies that looked adequate under a light-tailed hazard may leave the city badly exposed under heavier tails. This is why analyzing sensitivity to uncertain hazard parameters is essential before committing to a strategy.

Next Steps

You’ve completed a basic iCOW analysis. Here’s where to go from here:

Understand the math: Equations describes every cost and damage formula
Explore tail risk: Stochastic Walkthrough shows how stochastic dike failure creates hidden variance
Understand the design: Architecture explains the module structure and policy reparameterization

Setup

The Setting

Step 1: Configure the City

Step 2: Define the Hazard

Step 3: Try a Protection Strategy

Step 4: Run the Simulation

Step 5: Compare Strategies

Step 6: What if the Hazard Is Different?

Next Steps

References