Lab 5: Sea-Level Rise


CEVE 421/521


Fri., Feb. 16

There are two objectives of this lab:

  1. To familiarize ourselves with an increasingly complex model of our house-elevation problem
  2. To conduct exploratory modeling to understand the implications of different parameter values and how they affect our decision-making


The usual

As always:

  1. Clone the lab repository to your computer
  2. Open the lab repository in VS Code
  3. Open the Julia REPL and activate, then instantiate, the lab environment
  4. Make sure you can render: quarto render template.qmd in the terminal.
    • If you run into issues, try running ] build IJulia in the Julia REPL (] enters the package manager).
    • If you still have issues, try opening up That should trigger VS Code to give you the option to install the Python extension, which you should do. Then you should be able to open a menu in the bottom right of your screen to select which Python installation you want VS Code to use.

Load packages

using CSV
using DataFrames
using DataFramesMeta
using Distributions
using Plots
using StatsPlots
using Unitful

Precompiling CSV
  ✓ WeakRefStrings
  ✓ FilePathsBase
  ✓ CSV
  3 dependencies successfully precompiled in 8 seconds. 20 already precompiled.
Precompiling DataFramesMeta
  ✓ DataFramesMeta
  1 dependency successfully precompiled in 1 seconds. 29 already precompiled.
Precompiling Distributions
  ✓ StatsFuns
  ✓ StatsFuns → StatsFunsChainRulesCoreExt
  ✓ Distributions
  3 dependencies successfully precompiled in 4 seconds. 39 already precompiled.
Precompiling DistributionsTestExt
  ✓ Distributions → DistributionsTestExt
  ✓ Distributions → DistributionsChainRulesCoreExt
  2 dependencies successfully precompiled in 1 seconds. 43 already precompiled.
Precompiling StatsPlots
  ✓ AbstractFFTs → AbstractFFTsChainRulesCoreExt
  ✓ Distances → DistancesChainRulesCoreExt
  ✓ Adapt → AdaptStaticArraysExt
  ✓ StaticArrays → StaticArraysChainRulesCoreExt
  ✓ Interpolations
  ✓ Interpolations → InterpolationsUnitfulExt
  ✓ KernelDensity
  ✓ StatsPlots
  8 dependencies successfully precompiled in 7 seconds. 222 already precompiled.

Local package

We’re starting to accumulate a lot of code describing our model. A good way to store this model is by creating a local package. I have created a package called HouseElevation that contains the model code. You don’t need to do anything special to install it, and you don’t need to edit the code, though I’d encourage you to have a look around!

When we work with local packages, it’s common to use another package called Revise. This is a cool package that will automatically propagate any changes you make to the package to any code that uses the package. You don’t need to worry about this for now – just load them.

using Revise
using HouseElevation
Precompiling Revise
  ✓ JuliaInterpreter
  ✓ LoweredCodeUtils
  ✓ Revise
  3 dependencies successfully precompiled in 11 seconds. 4 already precompiled.
Precompiling HouseElevation
  ✓ HouseElevation
  1 dependency successfully precompiled in 2 seconds. 79 already precompiled.

Building the model

We’ve added a bit of complexity to our model. In this section, we walk through each of the sections of the model.


We will consider a single house, and will ignore uncertainty in the depth-damage function or other house parameters

  • Neglect uncertainty in depth-damage function
  • Consider a single building
  • We’re going to put all relevant information into a House object:
    • Depth-damage function
    • Area
    • Cost (USD)
    • Elevation relative to gauge
    • Metadata

We can create a House as follows – note that we’re using a let...end block to create the House object. This means that any variables defined inside the block are not available outside the block, which is a good way to avoid “polluting the global namespace.”

house = let
    haz_fl_dept ="data/haz_fl_dept.csv", DataFrame) # read in the file
    desc = "one story, Contents, fresh water, short duration"
    row = @rsubset(haz_fl_dept, :Description == desc)[1, :] # select the row I want
    area = 500u"ft^2"
    height_above_gauge = 10u"ft"

We can then use the House object to calculate the damage to the house for a given flood depth. Let’s convert the damage to dollars by multiplying the fraction (given by our depth-damage function) by the value of the house. For example:

    depths = uconvert.(u"ft", (-7.0u"ft"):(1.0u"inch"):(30.0u"ft"))
    damages = house.ddf.(depths) ./ 100
    damages_1000_usd = damages .* house.value_usd ./ 1000
        xlabel="Flood Depth",
        ylabel="Damage (Thousand USD)",
        size=(800, 400),
        yformatter=:plain, # prevents scientific notation

We can also use the House object to calculate the cost of raising the house to a given elevation. We use the elevation_cost function like this:

elevation_cost(house, 10u"ft")

and again we can plot this.

    elevations = 0u"ft":0.25u"ft":14u"ft"
    costs = [elevation_cost(house, eᵢ) for eᵢ in elevations]
        costs ./ 1_000;
        ylabel="Cost (Thousand USD)",
        size=(800, 400),
        yformatter=:plain, # prevents scientific notation

Sea-level rise

We will sample many different scenarios of sea-level rise

We’re modeling sea-level rise following the approach of Oddo et al. (2017). Essentially, we use five parameters: \(a\), \(b\), \(c\), \(t^*\), and \(c^*\). The local sea-level in year \(t\) is given by equation 6 of Oddo et al. (2017):

\[ \mathrm{SLR}= a + b(t - 2000) + c (t - 2000)^2 + c^* \, \mathbb{I} (t > t^*) (t - t^*) \]

The authors note:

In this model, the parameters \(a\), \(b\), and \(c\) represent the reasonably well-characterized process of thermosteric expansion as a second-order polynomial. It also accounts for more poorly understood processes, including potential abrupt sealevel rise consistent with sudden changes in ice flow dynamics. Here, \(c^*\) represents an increase in the rate of sea-level rise that takes place at some uncertain time, \(t^*\), in the future.

This is, of course, a highly simplified model. However, the parameters can be calibrated to match historical sea-level rise (i.e., throwing out any parameter values that don’t match the historical record) and use a statistical inversion method to estimate the parameters. One could also calibrate the parameters to match other, more complex, physics-based models. We’ll use Monte Carlo simulations from Oddo et al. (2017), available on GitHub. These were actually calibrated for the Netherlands, but we’ll pretend that sea-level rise in your location matches (which – as we know – it doesn’t).

slr_scenarios = let
    df ="data/slr_oddo.csv", DataFrame)
    [Oddo17SLR(a, b, c, tstar, cstar) for (a, b, c, tstar, cstar) in eachrow(df)]
println("There are $(length(slr_scenarios)) parameter sets")

We can plot these scenarios to get a sense of the range of sea-level rise we might expect.

    years = 1900:2150
    p = plot(;
        ylabel="Mean sea-level (ft)\nwith respect to the year 2000",
        label="Oddo et al. (2017)",
    for s in rand(slr_scenarios, 250)
        plot!(p, years, s.(years); color=:lightgrey, alpha=0.5, linewidth=0.5)