Case Study: Marshall Fire

The Marshall Fire ignited on December 30, 2021 near Boulder, Colorado during an extreme Chinook wind event. Sustained winds of 100+ mph drove the fire through grassland into suburban neighborhoods, burning over 6,000 acres and destroying more than 1,000 homes in just a few hours. It remains the most destructive wildfire in Colorado history by structure count.

This case study demonstrates how to integrate real geospatial data sources with Wildfires.jl’s component interface to simulate the Marshall Fire’s early spread.

Loading Geospatial Data

Pre-processed LANDFIRE terrain rasters and HRRR wind data are stored as GeoTIFF files. We load them with Rasters.jl and wrap them with GeoSurrogates.RasterWrap for fast bilinear interpolation.

# Load LANDFIRE terrain rasters
slope_raster     = Raster(joinpath(DATADIR, "slope.tif"))
aspect_raster    = Raster(joinpath(DATADIR, "aspect.tif"))
elevation_raster = Raster(joinpath(DATADIR, "elevation.tif"))
fuel_raster      = Raster(joinpath(DATADIR, "fuel.tif"))

# Wrap with bilinear interpolation for smooth spatial lookups
slope_wrap = RasterWrap(slope_raster)
aspect_wrap = RasterWrap(aspect_raster)
elevation_wrap = RasterWrap(elevation_raster)
fuel_wrap = RasterWrap(fuel_raster)
nothing

HRRR Wind Field

HRRR provides 10m u/v wind components and surface gusts at ~3 km resolution. We load nine snapshots at 15-minute intervals (17:00–19:00 UTC) from the NOAA HRRR archive on AWS — three hourly analyses plus sub-hourly forecasts — and build a WindField that interpolates in both space and time.

struct WindField
    times::Vector{DateTime}
    u_wraps::Vector{<:RasterWrap}
    v_wraps::Vector{<:RasterWrap}
    gust_wraps::Vector{<:RasterWrap}
end

function load_wind_field(datadir, suffixes, times)
    u_wraps    = [RasterWrap(Raster(joinpath(datadir, "wind_u_$s.tif"))) for s in suffixes]
    v_wraps    = [RasterWrap(Raster(joinpath(datadir, "wind_v_$s.tif"))) for s in suffixes]
    gust_wraps = [RasterWrap(Raster(joinpath(datadir, "wind_gust_$s.tif"))) for s in suffixes]
    WindField(times, u_wraps, v_wraps, gust_wraps)
end

function wind_at(wf::WindField, lon, lat, dt::DateTime)
    times = wf.times
    if dt <= first(times)
        return _wind_snapshot(wf, 1, lon, lat)
    elseif dt >= last(times)
        return _wind_snapshot(wf, length(times), lon, lat)
    end
    i = searchsortedlast(times, dt)
    t1, t2 = times[i], times[i + 1]
    alpha = Dates.value(dt - t1) / Dates.value(t2 - t1)
    w1 = _wind_snapshot(wf, i, lon, lat)
    w2 = _wind_snapshot(wf, i + 1, lon, lat)
    return (u    = (1 - alpha) * w1.u    + alpha * w2.u,
            v    = (1 - alpha) * w1.v    + alpha * w2.v,
            gust = (1 - alpha) * w1.gust + alpha * w2.gust)
end

function _wind_snapshot(wf::WindField, i, lon, lat)
    (u    = predict(wf.u_wraps[i], (lon, lat)),
     v    = predict(wf.v_wraps[i], (lon, lat)),
     gust = predict(wf.gust_wraps[i], (lon, lat)))
end

wind_field = load_wind_field(DATADIR, WIND_SUFFIXES, WIND_TIMES)
nothing

Custom Components

To plug real data into the simulation, we define custom AbstractFuel, AbstractWind, and AbstractTerrain subtypes. These convert from the simulation’s local meter-based grid to geographic (lon/lat) coordinates for spatial lookups.

using Wildfires.Rothermel: NFFL_MODELS, nffl_model

struct MarshallFuel <: AbstractFuel
    fuel_wrap::RasterWrap
    lon0::Float64
    lat0::Float64
end

const NON_FUEL_CODES = Set([92, 93, 98, 99])  # snow/ice, ag, water, barren

function (f::MarshallFuel)(t, x, y)
    lon = f.lon0 + x / M_PER_DEG_LON
    lat = f.lat0 + y / M_PER_DEG_LAT
    code = round(Int, predict(f.fuel_wrap, (lon, lat)))
    model = nffl_model(code)
    # Unrecognized codes (including 91/urban) fall back to SHORT_GRASS
    return model !== nothing ? model : SHORT_GRASS
end

function is_non_fuel(f::MarshallFuel, x, y)
    lon = f.lon0 + x / M_PER_DEG_LON
    lat = f.lat0 + y / M_PER_DEG_LAT
    code = round(Int, predict(f.fuel_wrap, (lon, lat)))
    return code ∈ NON_FUEL_CODES
end
nothing

struct MarshallWind <: AbstractWind
    field::WindField
    start_time::DateTime
    lon0::Float64
    lat0::Float64
    adjustment::Float64  # 10m wind → midflame adjustment
end

function (w::MarshallWind)(t, x, y)
    lon = w.lon0 + x / M_PER_DEG_LON
    lat = w.lat0 + y / M_PER_DEG_LAT
    dt = w.start_time + Second(round(Int, t * 60))
    wnd = wind_at(w.field, lon, lat, dt)
    # 10m mean wind speed (Rothermel is calibrated for sustained wind)
    speed_kmh = hypot(wnd.u, wnd.v) * 3.6 * w.adjustment
    # FROM direction in math convention (CCW from +x/east):
    # HRRR u = eastward, v = northward; negate to get FROM direction
    direction = atan(-wnd.v, -wnd.u)
    return (speed_kmh, direction)
end

struct MarshallTerrain <: AbstractTerrain
    slp_wrap::RasterWrap
    asp_wrap::RasterWrap
    lon0::Float64
    lat0::Float64
end

function (t::MarshallTerrain)(_, x, y)
    lon = t.lon0 + x / M_PER_DEG_LON
    lat = t.lat0 + y / M_PER_DEG_LAT
    slp_deg = predict(t.slp_wrap, (lon, lat))
    slope = tan(deg2rad(clamp(slp_deg, 0.0, 80.0)))
    # LANDFIRE aspect is azimuth (°CW from north) — convert to math convention (rad CCW from east)
    aspect = π/2 - deg2rad(predict(t.asp_wrap, (lon, lat)))
    return (slope, aspect)
end
nothing

Simulation

We set up a level set grid covering the Marshall Fire extent (~8.7 km x 6.3 km at 30 m resolution) and run a 2-hour simulation with spatially varying LANDFIRE fuel data and the Chinook wind field.

# Grid setup: 290 x 210 cells at 30m
nx, ny = 290, 210
dx = 30.0

# Offsets so ignition point (0, 0) maps to IGNITION_POINT in lon/lat
x0 = -(IGNITION_POINT.lon - EXTENT.X[1]) * M_PER_DEG_LON
y0 = -(IGNITION_POINT.lat - EXTENT.Y[1]) * M_PER_DEG_LAT

g = LevelSetGrid(nx, ny, dx=dx, x0=x0, y0=y0)
ignite!(g, 0.0, 0.0, 60.0)

# Components
fuel = MarshallFuel(fuel_wrap, IGNITION_POINT.lon, IGNITION_POINT.lat)
wind = MarshallWind(wind_field, DateTime(2021, 12, 30, 17), IGNITION_POINT.lon, IGNITION_POINT.lat, 0.31)
terrain = MarshallTerrain(slope_wrap, aspect_wrap, IGNITION_POINT.lon, IGNITION_POINT.lat)
moisture = UniformMoisture(FuelClasses(d1=0.03, d10=0.04, d100=0.05, herb=0.30, wood=0.50))

# Mark non-fuel cells (water, agriculture, snow/ice, barren) as unburnable
xs = collect(xcoords(g))
ys = collect(ycoords(g))
for j in eachindex(xs), i in eachindex(ys)
    if is_non_fuel(fuel, xs[j], ys[i])
        g.t_ignite[i, j] = NaN
    end
end

model = RothermelModel(fuel, wind, moisture, terrain, EllipticalBlending())
trace = Trace(g, 20)
nothing

println("Simulation complete: t = $(round(g.t, digits=1)) min, burned = $(count(<(0), g.φ)) cells")

Simulation complete: t = 44.8 min, burned = 4840 cells

Component Overview

After simulation, fireplot(grid, model) provides a multi-panel view of the fire state alongside all environmental inputs. The fire front (red contour) is overlaid on each panel for spatial reference.

fireplot(g, model)

Toggle individual panels or enable the spread rate field:

fireplot(g, model; terrain=false, spread_rate=true)

Animation

Fire spread animation with the observed perimeter (solid red) overlaid. The simulated perimeter is shown with a black outline and semi-transparent red fill over terrain elevation.

Results

Note

This simulation uses spatially varying LANDFIRE fuel data (NFFL 1–13 fuel models) with constant moisture and covers only the first ~2 hours of a fire that burned for over 24 hours. Wind data comes from HRRR 10m mean wind at 15-minute resolution, scaled to midflame height (adjustment factor 0.31). Non-fuel LANDFIRE cells (water, agriculture, snow/ice, barren) are marked as unburnable. LANDFIRE urban cells (code 91) are not in the NFFL fuel model table, so they fall back to SHORT_GRASS — a rough approximation that allows fire to spread through developed areas, consistent with the Marshall Fire’s rapid structure-to-structure ignition. The real fire’s spread was further intensified by extremely dry conditions (record low humidity), decades of grass accumulation in the wildland-urban interface, and ember transport — factors not captured in this basic Rothermel surface fire model.