List of submissions
It follows the code and plot of the RainGauge of all RainMaker submissions to /submissions sorted in alphabetical order of the filename.
A Nicusan: Evo opt unbounded (disqualified)
path: /submissions/anicusan_evoopt_unbounded.jl
rank: 1. of 14 submissions
author = "A Nicusan"
description = "Evo opt unbounded (disqualified)"
using SpeedyWeather, RainMaker
const PARAMETER_KEYS = (
:orography_scale, # [1], default: 1, scale of global orography
:mountain_height, # [m], default: 0, height of an additional mountain
:mountain_size, # [˚], default: 1, horizontal size of mountain
:mountain_lon, # [˚E], default: -80, longitude of that mountain
:mountain_lat, # [˚N], default: 40.45, latitude of that mountain
:temperature_equator, # [K], default: 300, sea surface temperature at the equator
:temperature_pole, # [K], default: 273, sea surfaec temperature at the poles
:temperature_usa, # [K], default: 0, land surface temperature anomaly over the USA
:temperature_pa, # [K], default: 0, land surface temperature anomaly in Pennsylvania
:zonal_wind, # [m/s], default: 35, zonal wind speed
)
const PARAMETER_DEFAULTS = [1, 0, 1, -80, 40.45, 300, 273, 0, 0, 35]
function max_precipitation(parameters::AbstractVector)
parameter_tuple = NamedTuple{PARAMETER_KEYS}(parameters)
return max_precipitation(parameter_tuple)
end
function max_precipitation(parameters::NamedTuple)
# define resolution. Use trunc=42, 63, 85, 127, ... for higher resolution, cubically slower
spectral_grid = SpectralGrid(trunc=31, nlayers=8)
# Define AquaPlanet ocean, for idealised sea surface temperatures
# but don't change land-sea mask = retain real ocean basins
ocean = AquaPlanet(spectral_grid,
temp_equator=parameters.temperature_equator,
temp_poles=parameters.temperature_pole)
land_temperature = ConstantLandTemperature(spectral_grid)
land = LandModel(spectral_grid; temperature=land_temperature)
initial_conditions = InitialConditions(
vordiv = ZonalWind(u₀=parameters.zonal_wind),
temp = JablonowskiTemperature(u₀=parameters.zonal_wind),
pres = PressureOnOrography(),
humid = ConstantRelativeHumidity())
orography = EarthOrography(spectral_grid, scale=parameters.orography_scale)
# construct model
model = PrimitiveWetModel(spectral_grid; ocean, land, initial_conditions, orography)
# Add rain gauge, locate in Pittsburgh PA
rain_gauge = RainGauge(spectral_grid, lond=-80, latd=40.45)
add!(model, rain_gauge)
# Initialize
simulation = initialize!(model, time=DateTime(2025, 7, 22))
# Add additional mountain
H = parameters.mountain_height
λ₀, φ₀, σ = parameters.mountain_lon, parameters.mountain_lat, parameters.mountain_size
set!(model, orography=(λ,φ) -> H*exp(-spherical_distance((λ,φ), (λ₀,φ₀), radius=360/2π)^2/2σ^2), add=true)
# land sea surface temperature anomalies
# 1. USA
set!(simulation, soil_temperature=
(λ, φ, k) -> (30 < φ < 50) && (240 < λ < 285) ? parameters.temperature_usa : 0, add=true)
# 2. Pennsylvania
A = parameters.temperature_pa
λ_az, φ_az, σ_az = -80, 40.45, 4 # location [˚], size [˚] of Azores
set!(simulation, soil_temperature=
(λ, φ, k) -> A*exp(-spherical_distance((λ,φ), (λ_az,φ_az), radius=360/2π)^2/2σ_az^2), add=true)
# Run simulation for 20 days
run!(simulation, period=Day(20))
# skip first 5 days, as is done in the RainMaker challenge
RainMaker.skip!(rain_gauge, Day(5))
# evaluate rain gauge
lsc = rain_gauge.accumulated_rain_large_scale
conv = rain_gauge.accumulated_rain_convection
total_precip = maximum(lsc) + maximum(conv)
return rain_gauge, total_precip
end
best_params = [
1.680189,
240.442519,
10.640524,
166.748171,
60.635959,
256.528781,
268.849961,
18.577591,
19.591949,
30.854462,
]
rain_gauge, total_precip = max_precipitation(best_params)
A Nicusan: Evo opt bounded
path: /submissions/anicusan_evoopt_bounded.jl
rank: 2. of 14 submissions
author = "A Nicusan"
description = "Evo opt bounded"
using SpeedyWeather, RainMaker
const PARAMETER_KEYS = (
:orography_scale, # [1], default: 1, scale of global orography
:mountain_height, # [m], default: 0, height of an additional mountain
:mountain_size, # [˚], default: 1, horizontal size of mountain
:mountain_lon, # [˚E], default: -80, longitude of that mountain
:mountain_lat, # [˚N], default: 40.45, latitude of that mountain
:temperature_equator, # [K], default: 300, sea surface temperature at the equator
:temperature_pole, # [K], default: 273, sea surfaec temperature at the poles
:temperature_usa, # [K], default: 0, land surface temperature anomaly over the USA
:temperature_pa, # [K], default: 0, land surface temperature anomaly in Pennsylvania
:zonal_wind, # [m/s], default: 35, zonal wind speed
)
const PARAMETER_DEFAULTS = [1, 0, 1, -80, 40.45, 300, 273, 0, 0, 35]
function max_precipitation(parameters::AbstractVector)
parameter_tuple = NamedTuple{PARAMETER_KEYS}(parameters)
return max_precipitation(parameter_tuple)
end
function max_precipitation(parameters::NamedTuple)
# define resolution. Use trunc=42, 63, 85, 127, ... for higher resolution, cubically slower
spectral_grid = SpectralGrid(trunc=31, nlayers=8)
# Define AquaPlanet ocean, for idealised sea surface temperatures
# but don't change land-sea mask = retain real ocean basins
ocean = AquaPlanet(spectral_grid,
temp_equator=parameters.temperature_equator,
temp_poles=parameters.temperature_pole)
land_temperature = ConstantLandTemperature(spectral_grid)
land = LandModel(spectral_grid; temperature=land_temperature)
initial_conditions = InitialConditions(
vordiv = ZonalWind(u₀=parameters.zonal_wind),
temp = JablonowskiTemperature(u₀=parameters.zonal_wind),
pres = PressureOnOrography(),
humid = ConstantRelativeHumidity())
orography = EarthOrography(spectral_grid, scale=parameters.orography_scale)
# construct model
model = PrimitiveWetModel(spectral_grid; ocean, land, initial_conditions, orography)
# Add rain gauge, locate in Pittsburgh PA
rain_gauge = RainGauge(spectral_grid, lond=-80, latd=40.45)
add!(model, rain_gauge)
# Initialize
simulation = initialize!(model, time=DateTime(2025, 7, 22))
# Add additional mountain
H = parameters.mountain_height
λ₀, φ₀, σ = parameters.mountain_lon, parameters.mountain_lat, parameters.mountain_size
set!(model, orography=(λ,φ) -> H*exp(-spherical_distance((λ,φ), (λ₀,φ₀), radius=360/2π)^2/2σ^2), add=true)
# land sea surface temperature anomalies
# 1. USA
set!(simulation, soil_temperature=
(λ, φ, k) -> (30 < φ < 50) && (240 < λ < 285) ? parameters.temperature_usa : 0, add=true)
# 2. Pennsylvania
A = parameters.temperature_pa
λ_az, φ_az, σ_az = -80, 40.45, 4 # location [˚], size [˚] of Azores
set!(simulation, soil_temperature=
(λ, φ, k) -> A*exp(-spherical_distance((λ,φ), (λ_az,φ_az), radius=360/2π)^2/2σ_az^2), add=true)
# Run simulation for 20 days
run!(simulation, period=Day(20))
# skip first 5 days, as is done in the RainMaker challenge
RainMaker.skip!(rain_gauge, Day(5))
# evaluate rain gauge
lsc = rain_gauge.accumulated_rain_large_scale
conv = rain_gauge.accumulated_rain_convection
total_precip = maximum(lsc) + maximum(conv)
return rain_gauge, total_precip
end
best_params = [
1.999196,
2304.509728,
14.889662,
74.446324,
-70.338707,
280.224610,
299.998274,
4.999593,
4.996833,
33.312800,
]
rain_gauge, total_precip = max_precipitation(best_params)
Chris Rackauckas: Brute Force OptimizationBBO DE
path: /submissions/chris_brute_force.jl
rank: 3. of 14 submissions
author = "Chris Rackauckas"
description = "Brute Force OptimizationBBO DE"
#=
using RainMakerChallenge2025
using OptimizationBBO, Optimization
function objective(params, p)
return -max_precipitation(params)
end
lower_bounds = [0., -2000., 0., -180., -90., 270., 270., -5., -5., 5.]
upper_bounds = [2., 5000., 30., 180., 90., 300., 300., 5., 5., 50.]
optf = OptimizationFunction(objective)
prob = Optimization.OptimizationProblem(optf, [1, 0, 1, -80, 40.45, 300, 273, 0, 0, 35],
lb=lower_bounds, ub=upper_bounds)
println("Starting optimization...")
sol = solve(prob, BBO_adaptive_de_rand_1_bin_radiuslimited(); maxiters=10_000)
println("Optimized parameters: ", sol.u)
println("Predicted max precipitation from surrogate: ", max_precipitation(sol.u))
=#
using SpeedyWeather, RainMaker
const PARAMETER_KEYS = (
:orography_scale, # [1], default: 1, scale of global orography
:mountain_height, # [m], default: 0, height of an additional azores mountain
:mountain_size, # [˚], default: 1, horizontal size of mountain in Pittsburgh
:mountain_lon, # [˚E], default: -80, longitude of that mountain in Pittsburgh
:mountain_lat, # [˚N], default: 40.45, latitude of that mountain in Pittsburgh
:temperature_equator, # [K], default: 300, sea surface temperature at the equator
:temperature_pole, # [K], default: 273, sea surfaec temperature at the poles
:temperature_usa, # [K], default: 0, land surface temperature anomaly over the USA
:temperature_pa, # [K], default: 0, land surface temperature anomaly in Pennsylvania
:zonal_wind, # [m/s], default: 35, zonal wind speed
)
const PARAMETER_DEFAULTS = [1, 0, 1, -80, 40.45, 300, 273, 0, 0, 35]
function max_precipitation(parameters::AbstractVector)
parameter_tuple = NamedTuple{PARAMETER_KEYS}(parameters)
return max_precipitation(parameter_tuple)
end
function max_precipitation(parameters::NamedTuple)
# define resolution. Use trunc=42, 63, 85, 127, ... for higher resolution, cubically slower
spectral_grid = SpectralGrid(trunc=31, nlayers=8)
# Define AquaPlanet ocean, for idealised sea surface temperatures
# but don't change land-sea mask = retain real ocean basins
ocean = AquaPlanet(spectral_grid,
temp_equator=parameters.temperature_equator,
temp_poles=parameters.temperature_pole)
land_temperature = ConstantLandTemperature(spectral_grid)
land = LandModel(spectral_grid; temperature=land_temperature)
initial_conditions = InitialConditions(
vordiv = ZonalWind(u₀=parameters.zonal_wind),
temp = JablonowskiTemperature(u₀=parameters.zonal_wind),
pres = PressureOnOrography(),
humid = ConstantRelativeHumidity())
orography = EarthOrography(spectral_grid, scale=parameters.orography_scale)
# construct model
model = PrimitiveWetModel(spectral_grid; ocean, land, initial_conditions, orography)
# Add rain gauge, locate in Pittsburgh PA
rain_gauge = RainGauge(spectral_grid, lond=-80, latd=40.45)
add!(model, rain_gauge)
# Initialize
simulation = initialize!(model, time=DateTime(2025, 7, 22))
# Add additional mountain
H = parameters.mountain_height
λ₀, φ₀, σ = parameters.mountain_lon, parameters.mountain_lat, parameters.mountain_size
set!(model, orography=(λ,φ) -> H*exp(-spherical_distance((λ,φ), (λ₀,φ₀), radius=360/2π)^2/2σ^2), add=true)
# land sea surface temperature anomalies
# 1. USA
set!(simulation, soil_temperature=
(λ, φ, k) -> (30 < φ < 50) && (240 < λ < 285) ? parameters.temperature_usa : 0, add=true)
# 2. Pennsylvania
A = parameters.temperature_pa
λ_az, φ_az, σ_az = -80, 40.45, 4 # location [˚], size [˚] of Azores
set!(simulation, soil_temperature=
(λ, φ, k) -> A*exp(-spherical_distance((λ,φ), (λ_az,φ_az), radius=360/2π)^2/2σ_az^2), add=true)
# Run simulation for 20 days
run!(simulation, period=Day(20))
# skip first 5 days, as is done in the RainMaker challenge
RainMaker.skip!(rain_gauge, Day(5))
# evaluate rain gauge
lsc = rain_gauge.accumulated_rain_large_scale
conv = rain_gauge.accumulated_rain_convection
total_precip = maximum(lsc) + maximum(conv)
return rain_gauge, total_precip
end
final_params = [1.9851412873870289, 2485.737125917869, 14.118878320714353, 73.66989910502846, -72.40588406991382, 280.9029383973223, 299.2432527466829, 4.773285514316327, 4.586811470190235, 33.398388938656204]
rain_gauge, total_precip = max_precipitation(final_params)
total_precip
Anas: NN Surrogate
path: /submissions/anas-nn_surrogate_v2.jl
rank: 4. of 14 submissions
author = "Anas"
description = "NN Surrogate"
using SpeedyWeather, RainMaker
const PARAMETER_KEYS = (
:orography_scale, # [1], default: 1, scale of global orography
:mountain_height, # [m], default: 0, height of an additional azores mountain
:mountain_size, # [˚], default: 1, horizontal size of mountain in Pittsburgh
:mountain_lon, # [˚E], default: -80, longitude of that mountain in Pittsburgh
:mountain_lat, # [˚N], default: 40.45, latitude of that mountain in Pittsburgh
:temperature_equator, # [K], default: 300, sea surface temperature at the equator
:temperature_pole, # [K], default: 273, sea surfaec temperature at the poles
:temperature_usa, # [K], default: 0, land surface temperature anomaly over the USA
:temperature_pa, # [K], default: 0, land surface temperature anomaly in Pennsylvania
:zonal_wind, # [m/s], default: 35, zonal wind speed
)
const PARAMETER_DEFAULTS = [1, 0, 1, -80, 40.45, 300, 273, 0, 0, 35]
function max_precipitation(parameters::AbstractVector)
parameter_tuple = NamedTuple{PARAMETER_KEYS}(parameters)
return max_precipitation(parameter_tuple)
end
function max_precipitation(parameters::NamedTuple)
# define resolution. Use trunc=42, 63, 85, 127, ... for higher resolution, cubically slower
spectral_grid = SpectralGrid(trunc=31, nlayers=8)
# Define AquaPlanet ocean, for idealised sea surface temperatures
# but don't change land-sea mask = retain real ocean basins
ocean = AquaPlanet(spectral_grid,
temp_equator=parameters.temperature_equator,
temp_poles=parameters.temperature_pole)
land_temperature = ConstantLandTemperature(spectral_grid)
land = LandModel(spectral_grid; temperature=land_temperature)
initial_conditions = InitialConditions(
vordiv = ZonalWind(u₀=parameters.zonal_wind),
temp = JablonowskiTemperature(u₀=parameters.zonal_wind),
pres = PressureOnOrography(),
humid = ConstantRelativeHumidity())
orography = EarthOrography(spectral_grid, scale=parameters.orography_scale)
# construct model
model = PrimitiveWetModel(spectral_grid; ocean, land, initial_conditions, orography)
# Add rain gauge, locate in Pittsburgh PA
rain_gauge = RainGauge(spectral_grid, lond=-80, latd=40.45)
add!(model, rain_gauge)
# Initialize
simulation = initialize!(model, time=DateTime(2025, 7, 22))
# Add additional mountain
H = parameters.mountain_height
λ₀, φ₀, σ = parameters.mountain_lon, parameters.mountain_lat, parameters.mountain_size
set!(model, orography=(λ,φ) -> H*exp(-spherical_distance((λ,φ), (λ₀,φ₀), radius=360/2π)^2/2σ^2), add=true)
# land sea surface temperature anomalies
# 1. USA
set!(simulation, soil_temperature=
(λ, φ, k) -> (30 < φ < 50) && (240 < λ < 285) ? parameters.temperature_usa : 0, add=true)
# 2. Pennsylvania
A = parameters.temperature_pa
λ_az, φ_az, σ_az = -80, 40.45, 4 # location [˚], size [˚] of Azores
set!(simulation, soil_temperature=
(λ, φ, k) -> A*exp(-spherical_distance((λ,φ), (λ_az,φ_az), radius=360/2π)^2/2σ_az^2), add=true)
# Run simulation for 20 days
run!(simulation, period=Day(20))
# skip first 5 days, as is done in the RainMaker challenge
RainMaker.skip!(rain_gauge, Day(5))
# evaluate rain gauge
lsc = rain_gauge.accumulated_rain_large_scale
conv = rain_gauge.accumulated_rain_convection
total_precip = maximum(lsc) + maximum(conv)
return rain_gauge, total_precip
end
final_params =
[ 1.1229311227798462,
3772.9730010032654,
4.986201524734497,
-120.26883065700531,
49.61676836013794,
276.7463958263397,
282.4692580103874,
4.013188481330872,
3.889443278312683,
26.725226491689682 ]
rain_gauge, total_precip = max_precipitation(final_params)
total_precip
Anas: 200mm+ RainMaker
path: /submissions/anas_submission.jl
rank: 5. of 14 submissions
author = "Anas"
description = "200mm+ RainMaker"
using SpeedyWeather, RainMaker
const PARAMETER_KEYS = (
:orography_scale, # [1], default: 1, scale of global orography
:mountain_height, # [m], default: 0, height of an additional azores mountain
:mountain_size, # [˚], default: 1, horizontal size of mountain in Pittsburgh
:mountain_lon, # [˚E], default: -80, longitude of that mountain in Pittsburgh
:mountain_lat, # [˚N], default: 40.45, latitude of that mountain in Pittsburgh
:temperature_equator, # [K], default: 300, sea surface temperature at the equator
:temperature_pole, # [K], default: 273, sea surfaec temperature at the poles
:temperature_usa, # [K], default: 0, land surface temperature anomaly over the USA
:temperature_pa, # [K], default: 0, land surface temperature anomaly in Pennsylvania
:zonal_wind, # [m/s], default: 35, zonal wind speed
)
const PARAMETER_DEFAULTS = [1, 0, 1, -80, 40.45, 300, 273, 0, 0, 35]
function max_precipitation(parameters::AbstractVector)
parameter_tuple = NamedTuple{PARAMETER_KEYS}(parameters)
return max_precipitation(parameter_tuple)
end
function max_precipitation(parameters::NamedTuple)
# define resolution. Use trunc=42, 63, 85, 127, ... for higher resolution, cubically slower
spectral_grid = SpectralGrid(trunc=31, nlayers=8)
# Define AquaPlanet ocean, for idealised sea surface temperatures
# but don't change land-sea mask = retain real ocean basins
ocean = AquaPlanet(spectral_grid,
temp_equator=parameters.temperature_equator,
temp_poles=parameters.temperature_pole)
land_temperature = ConstantLandTemperature(spectral_grid)
land = LandModel(spectral_grid; temperature=land_temperature)
initial_conditions = InitialConditions(
vordiv = ZonalWind(u₀=parameters.zonal_wind),
temp = JablonowskiTemperature(u₀=parameters.zonal_wind),
pres = PressureOnOrography(),
humid = ConstantRelativeHumidity())
orography = EarthOrography(spectral_grid, scale=parameters.orography_scale)
# construct model
model = PrimitiveWetModel(spectral_grid; ocean, land, initial_conditions, orography)
# Add rain gauge, locate in Pittsburgh PA
rain_gauge = RainGauge(spectral_grid, lond=-80, latd=40.45)
add!(model, rain_gauge)
# Initialize
simulation = initialize!(model, time=DateTime(2025, 7, 22))
# Add additional mountain
H = parameters.mountain_height
λ₀, φ₀, σ = parameters.mountain_lon, parameters.mountain_lat, parameters.mountain_size
set!(model, orography=(λ,φ) -> H*exp(-spherical_distance((λ,φ), (λ₀,φ₀), radius=360/2π)^2/2σ^2), add=true)
# land sea surface temperature anomalies
# 1. USA
set!(simulation, soil_temperature=
(λ, φ, k) -> (30 < φ < 50) && (240 < λ < 285) ? parameters.temperature_usa : 0, add=true)
# 2. Pennsylvania
A = parameters.temperature_pa
λ_az, φ_az, σ_az = -80, 40.45, 4 # location [˚], size [˚] of Azores
set!(simulation, soil_temperature=
(λ, φ, k) -> A*exp(-spherical_distance((λ,φ), (λ_az,φ_az), radius=360/2π)^2/2σ_az^2), add=true)
# Run simulation for 20 days
run!(simulation, period=Day(20))
# skip first 5 days, as is done in the RainMaker challenge
RainMaker.skip!(rain_gauge, Day(5))
# evaluate rain gauge
lsc = rain_gauge.accumulated_rain_large_scale
conv = rain_gauge.accumulated_rain_convection
total_precip = maximum(lsc) + maximum(conv)
return rain_gauge, total_precip
end
final_params =
[2.0163631439208984,
5013.73028755188,
-0.017217400018125772,
-22.203487157821655,
-19.583306908607483,
282.971847653389,
300.2009081840515,
2.134748101234436,
4.705338478088379,
34.92057800292969 ]
rain_gauge, total_precip = max_precipitation(final_params)
total_precip
Gabriel Konar-Steenberg: Brute force built on Anas surrogate
path: /submissions/gabrielks_brute_force_1.jl
rank: 6. of 14 submissions
author = "Gabriel Konar-Steenberg"
description = "Brute force built on Anas surrogate"
using SpeedyWeather, RainMaker
"""
Quickly written code based on
https://github.com/AnasAbdelR/RainMakerChallenge2025.jl/blob/main/examples/surrogate_nn_ex.jl.
We load the pretrained surrogate model and do a brute force search over many optimization
starting points for the one that results in the best performance on the real model.
"""
#=
using JLSO
using LinearAlgebra
using Distributions
using NNlib
using Random
using Optimization
using OptimizationOptimJL
using Flux
using SpeedyWeather
using RainMaker
using RainMakerChallenge2025
using Logging
minmaxnorm(data, lb, ub, norm_min=0., norm_max=1.) = @. norm_min + (data - lb) * (norm_max - norm_min) / (ub - lb)
minmaxdenorm(data, lb, ub, norm_min=0., norm_max=1.) = @. lb + (data - norm_min) * (ub - lb) / (norm_max - norm_min)
function outer_iter(seed)
path = joinpath(@__DIR__, "data", "10kdata.jlso")
data = JLSO.load(path)
output_data = data[:d][:outputs]
inputs_lb = [0., -2000., 0., -180., -90., 270., 270., -5., -5., 5.]
inputs_ub = [2., 5000., 30., 180., 90., 300., 300., 5., 5., 50.]
output_data = reshape(output_data, (1,:))
outputs_lb, outputs_ub = extrema(output_data)
# Load the surrogate model
path = joinpath(@__DIR__, "models", "surrogate_nn_ex_128_500-v2.jlso")
surrogate = JLSO.load(path)[:surrogate]
# Simple SciML Optimization equivalent
function objective_func_simple(x, p)
X = reshape(x, 10, :) |> gpu # Same as your Flux code
y_pred = surrogate(X)
obj = -sum(y_pred) # Same objective
λ = 100.0f0
lower_violation = max.(0.0f0, -x)
upper_violation = max.(0.0f0, x .- 1.0f0)
penalty = λ * (sum(lower_violation.^2) + sum(upper_violation.^2))
return obj + penalty
end
# Same initial setup as your Flux code
Random.seed!(seed)
x0 = vec(rand(Float32, 10, 10)) # Flatten to vector for Optimization.jl
optf = Optimization.OptimizationFunction(objective_func_simple, Optimization.AutoZygote())
prob = Optimization.OptimizationProblem(optf, x0)
# println("Starting SciML optimization (1000 samples)...")
sol = solve(prob, LBFGS(), maxiters=20000)
# Extract results exactly like your Flux code
X_candidate_sciml = reshape(sol.u, 10, :) |> gpu
idx = argmax(surrogate(X_candidate_sciml))
best_param_sample_sciml = minmaxdenorm(X_candidate_sciml[:,idx[2]]|>cpu, inputs_lb, inputs_ub)
# @show best_param_sample_sciml
pred_vals_sciml = minmaxdenorm(surrogate(X_candidate_sciml)[:,idx[2]]|>cpu, outputs_lb, outputs_ub)
# @show pred_vals_sciml
actual_precip_sciml = max_precipitation(best_param_sample_sciml)
# @show actual_precip_sciml
return best_param_sample_sciml, pred_vals_sciml, actual_precip_sciml
end
function main()
global_logger(SimpleLogger(stderr, Logging.Warn))
Threads.@threads for i in 1:1000
best_param_sample_sciml, pred_vals_sciml, actual_precip_sciml = outer_iter(i)
println("Iteration $i")
if actual_precip_sciml > 210.0
println("Best Parameters: $best_param_sample_sciml\nPredicted Values: $pred_vals_sciml\nActual Precipitation: $actual_precip_sciml")
end
end
end
main()
=#
const PARAMETER_KEYS = (
:orography_scale, # [1], default: 1, scale of global orography
:mountain_height, # [m], default: 0, height of an additional azores mountain
:mountain_size, # [˚], default: 1, horizontal size of mountain in Pittsburgh
:mountain_lon, # [˚E], default: -80, longitude of that mountain in Pittsburgh
:mountain_lat, # [˚N], default: 40.45, latitude of that mountain in Pittsburgh
:temperature_equator, # [K], default: 300, sea surface temperature at the equator
:temperature_pole, # [K], default: 273, sea surfaec temperature at the poles
:temperature_usa, # [K], default: 0, land surface temperature anomaly over the USA
:temperature_pa, # [K], default: 0, land surface temperature anomaly in Pennsylvania
:zonal_wind, # [m/s], default: 35, zonal wind speed
)
const PARAMETER_DEFAULTS = [1, 0, 1, -80, 40.45, 300, 273, 0, 0, 35]
function max_precipitation(parameters::AbstractVector)
parameter_tuple = NamedTuple{PARAMETER_KEYS}(parameters)
return max_precipitation(parameter_tuple)
end
function max_precipitation(parameters::NamedTuple)
# define resolution. Use trunc=42, 63, 85, 127, ... for higher resolution, cubically slower
spectral_grid = SpectralGrid(trunc=31, nlayers=8)
# Define AquaPlanet ocean, for idealised sea surface temperatures
# but don't change land-sea mask = retain real ocean basins
ocean = AquaPlanet(spectral_grid,
temp_equator=parameters.temperature_equator,
temp_poles=parameters.temperature_pole)
land_temperature = ConstantLandTemperature(spectral_grid)
land = LandModel(spectral_grid; temperature=land_temperature)
initial_conditions = InitialConditions(
vordiv = ZonalWind(u₀=parameters.zonal_wind),
temp = JablonowskiTemperature(u₀=parameters.zonal_wind),
pres = PressureOnOrography(),
humid = ConstantRelativeHumidity())
orography = EarthOrography(spectral_grid, scale=parameters.orography_scale)
# construct model
model = PrimitiveWetModel(spectral_grid; ocean, land, initial_conditions, orography)
# Add rain gauge, locate in Pittsburgh PA
rain_gauge = RainGauge(spectral_grid, lond=-80, latd=40.45)
add!(model, rain_gauge)
# Initialize
simulation = initialize!(model, time=DateTime(2025, 7, 22))
# Add additional mountain
H = parameters.mountain_height
λ₀, φ₀, σ = parameters.mountain_lon, parameters.mountain_lat, parameters.mountain_size
set!(model, orography=(λ,φ) -> H*exp(-spherical_distance((λ,φ), (λ₀,φ₀), radius=360/2π)^2/2σ^2), add=true)
# land sea surface temperature anomalies
# 1. USA
set!(simulation, soil_temperature=
(λ, φ, k) -> (30 < φ < 50) && (240 < λ < 285) ? parameters.temperature_usa : 0, add=true)
# 2. Pennsylvania
A = parameters.temperature_pa
λ_az, φ_az, σ_az = -80, 40.45, 4 # location [˚], size [˚] of Azores
set!(simulation, soil_temperature=
(λ, φ, k) -> A*exp(-spherical_distance((λ,φ), (λ_az,φ_az), radius=360/2π)^2/2σ_az^2), add=true)
# Run simulation for 20 days
run!(simulation, period=Day(20))
# skip first 5 days, as is done in the RainMaker challenge
RainMaker.skip!(rain_gauge, Day(5))
# evaluate rain gauge
lsc = rain_gauge.accumulated_rain_large_scale
conv = rain_gauge.accumulated_rain_convection
total_precip = maximum(lsc) + maximum(conv)
return rain_gauge, total_precip
end
final_params =
[2.0163488388061523,
5013.733625411987,
-0.017316662706434727,
-22.202972173690796,
-19.587678909301758,
282.97233670949936,
300.200936794281,
2.133525013923645,
4.704113602638245,
34.923005402088165]
rain_gauge, total_precip = max_precipitation(final_params)
total_precip # produces 212.7933 on RainMakerChallenge2025.jl's evaluator, 211.040 on this evaluator -- TODO investigate the discrepancy
Theo D: Kriging Surrogate with Latin Hypercube Sampling
path: /submissions/theo_juliacon.jl
rank: 7. of 14 submissions
author = "Theo D"
description = "Kriging Surrogate with Latin Hypercube Sampling"
#=
using JLSO
using EvoTrees
using OptimizationBBO, Optimization
path = joinpath(dirname(@__DIR__), "data", "10kdata.jlso")
data = JLSO.load(path)
x = data[:d].inputs
y = data[:d].outputs
X_train = x'
y_train = vec(y)
config = EvoTreeRegressor(
nrounds=200,
max_depth=4,
eta=0.1,
rng=123
)
println("Training EvoTrees model...")
model = EvoTrees.fit_evotree(config; x_train=X_train, y_train=y_train)
println("Model trained successfully")
function surrogate_model(params)
pred = EvoTrees.predict(model, reshape(params, 1, :))
return Float64(pred[1])
end
function objective(params, p)
return -surrogate_model(params)
end
lower_bounds = [0., -2000., 0., -180., -90., 270., 270., -5., -5., 5.]
upper_bounds = [2., 5000., 30., 180., 90., 300., 300., 5., 5., 50.]
optf = OptimizationFunction(objective)
prob = Optimization.OptimizationProblem(optf, [1, 0, 1, -80, 40.45, 300, 273, 0, 0, 35],
lb=lower_bounds, ub=upper_bounds)
println("Starting optimization...")
sol = solve(prob, BBO_adaptive_de_rand_1_bin_radiuslimited(); maxtime=60)
println("Optimized parameters: ", sol.u)
println("Predicted max precipitation from surrogate: ", surrogate_model(sol.u))
=#
using SpeedyWeather, RainMaker
const PARAMETER_KEYS = (
:orography_scale, # [1], default: 1, scale of global orography
:mountain_height, # [m], default: 0, height of an additional azores mountain
:mountain_size, # [˚], default: 1, horizontal size of mountain in Pittsburgh
:mountain_lon, # [˚E], default: -80, longitude of that mountain in Pittsburgh
:mountain_lat, # [˚N], default: 40.45, latitude of that mountain in Pittsburgh
:temperature_equator, # [K], default: 300, sea surface temperature at the equator
:temperature_pole, # [K], default: 273, sea surfaec temperature at the poles
:temperature_usa, # [K], default: 0, land surface temperature anomaly over the USA
:temperature_pa, # [K], default: 0, land surface temperature anomaly in Pennsylvania
:zonal_wind, # [m/s], default: 35, zonal wind speed
)
const PARAMETER_DEFAULTS = [1, 0, 1, -80, 40.45, 300, 273, 0, 0, 35]
function max_precipitation(parameters::AbstractVector)
parameter_tuple = NamedTuple{PARAMETER_KEYS}(parameters)
return max_precipitation(parameter_tuple)
end
function max_precipitation(parameters::NamedTuple)
# define resolution. Use trunc=42, 63, 85, 127, ... for higher resolution, cubically slower
spectral_grid = SpectralGrid(trunc=31, nlayers=8)
# Define AquaPlanet ocean, for idealised sea surface temperatures
# but don't change land-sea mask = retain real ocean basins
ocean = AquaPlanet(spectral_grid,
temp_equator=parameters.temperature_equator,
temp_poles=parameters.temperature_pole)
land_temperature = ConstantLandTemperature(spectral_grid)
land = LandModel(spectral_grid; temperature=land_temperature)
initial_conditions = InitialConditions(
vordiv = ZonalWind(u₀=parameters.zonal_wind),
temp = JablonowskiTemperature(u₀=parameters.zonal_wind),
pres = PressureOnOrography(),
humid = ConstantRelativeHumidity())
orography = EarthOrography(spectral_grid, scale=parameters.orography_scale)
# construct model
model = PrimitiveWetModel(spectral_grid; ocean, land, initial_conditions, orography)
# Add rain gauge, locate in Pittsburgh PA
rain_gauge = RainGauge(spectral_grid, lond=-80, latd=40.45)
add!(model, rain_gauge)
# Initialize
simulation = initialize!(model, time=DateTime(2025, 7, 22))
# Add additional mountain
H = parameters.mountain_height
λ₀, φ₀, σ = parameters.mountain_lon, parameters.mountain_lat, parameters.mountain_size
set!(model, orography=(λ,φ) -> H*exp(-spherical_distance((λ,φ), (λ₀,φ₀), radius=360/2π)^2/2σ^2), add=true)
# land sea surface temperature anomalies
# 1. USA
set!(simulation, soil_temperature=
(λ, φ, k) -> (30 < φ < 50) && (240 < λ < 285) ? parameters.temperature_usa : 0, add=true)
# 2. Pennsylvania
A = parameters.temperature_pa
λ_az, φ_az, σ_az = -80, 40.45, 4 # location [˚], size [˚] of Azores
set!(simulation, soil_temperature=
(λ, φ, k) -> A*exp(-spherical_distance((λ,φ), (λ_az,φ_az), radius=360/2π)^2/2σ_az^2), add=true)
# Run simulation for 20 days
run!(simulation, period=Day(20))
# skip first 5 days, as is done in the RainMaker challenge
RainMaker.skip!(rain_gauge, Day(5))
# evaluate rain gauge
lsc = rain_gauge.accumulated_rain_large_scale
conv = rain_gauge.accumulated_rain_convection
total_precip = maximum(lsc) + maximum(conv)
return rain_gauge, total_precip
end
final_params = [1.9728003644919363, 2506.981756969998, 10.799751688314835, -119.61854020127022, 55.932955161047225, 287.9571828050346, 291.00468964730595, 0.5460025649020176, 4.157291328686813, 25.87763330165432]
rain_gauge, total_precip = max_precipitation(final_params)
total_precip
Cameron Bieganek: Latin hypercube search, no surrogate.
path: /submissions/bieganek_latin_hypercube.jl
rank: 8. of 14 submissions
author = "Cameron Bieganek"
description = "Latin hypercube search, no surrogate."
#=
using RainMakerChallenge2025
using Surrogates
lb = [0., -2000., 0., -180., -90., 270., 270., -5., -5., 5.]
ub = [2., 5000., 30., 180., 90., 300., 300., 5., 5., 50.]
n = 800
params = collect.(sample(n, lb, ub, LatinHypercubeSample()))
rainfall = max_precipitation.(params)
maximum(rainfall)
i = argmax(rainfall)
best_params = params[i]
=#
using SpeedyWeather, RainMaker
const PARAMETER_KEYS = (
:orography_scale, # [1], default: 1, scale of global orography
:mountain_height, # [m], default: 0, height of an additional azores mountain
:mountain_size, # [˚], default: 1, horizontal size of mountain in Pittsburgh
:mountain_lon, # [˚E], default: -80, longitude of that mountain in Pittsburgh
:mountain_lat, # [˚N], default: 40.45, latitude of that mountain in Pittsburgh
:temperature_equator, # [K], default: 300, sea surface temperature at the equator
:temperature_pole, # [K], default: 273, sea surfaec temperature at the poles
:temperature_usa, # [K], default: 0, land surface temperature anomaly over the USA
:temperature_pa, # [K], default: 0, land surface temperature anomaly in Pennsylvania
:zonal_wind, # [m/s], default: 35, zonal wind speed
)
const PARAMETER_DEFAULTS = [1, 0, 1, -80, 40.45, 300, 273, 0, 0, 35]
function max_precipitation(parameters::AbstractVector)
parameter_tuple = NamedTuple{PARAMETER_KEYS}(parameters)
return max_precipitation(parameter_tuple)
end
function max_precipitation(parameters::NamedTuple)
# define resolution. Use trunc=42, 63, 85, 127, ... for higher resolution, cubically slower
spectral_grid = SpectralGrid(trunc=31, nlayers=8)
# Define AquaPlanet ocean, for idealised sea surface temperatures
# but don't change land-sea mask = retain real ocean basins
ocean = AquaPlanet(spectral_grid,
temp_equator=parameters.temperature_equator,
temp_poles=parameters.temperature_pole)
land_temperature = ConstantLandTemperature(spectral_grid)
land = LandModel(spectral_grid; temperature=land_temperature)
initial_conditions = InitialConditions(
vordiv = ZonalWind(u₀=parameters.zonal_wind),
temp = JablonowskiTemperature(u₀=parameters.zonal_wind),
pres = PressureOnOrography(),
humid = ConstantRelativeHumidity())
orography = EarthOrography(spectral_grid, scale=parameters.orography_scale)
# construct model
model = PrimitiveWetModel(spectral_grid; ocean, land, initial_conditions, orography)
# Add rain gauge, locate in Pittsburgh PA
rain_gauge = RainGauge(spectral_grid, lond=-80, latd=40.45)
add!(model, rain_gauge)
# Initialize
simulation = initialize!(model, time=DateTime(2025, 7, 22))
# Add additional mountain
H = parameters.mountain_height
λ₀, φ₀, σ = parameters.mountain_lon, parameters.mountain_lat, parameters.mountain_size
set!(model, orography=(λ,φ) -> H*exp(-spherical_distance((λ,φ), (λ₀,φ₀), radius=360/2π)^2/2σ^2), add=true)
# land sea surface temperature anomalies
# 1. USA
set!(simulation, soil_temperature=
(λ, φ, k) -> (30 < φ < 50) && (240 < λ < 285) ? parameters.temperature_usa : 0, add=true)
# 2. Pennsylvania
A = parameters.temperature_pa
λ_az, φ_az, σ_az = -80, 40.45, 4 # location [˚], size [˚] of Azores
set!(simulation, soil_temperature=
(λ, φ, k) -> A*exp(-spherical_distance((λ,φ), (λ_az,φ_az), radius=360/2π)^2/2σ_az^2), add=true)
# Run simulation for 20 days
run!(simulation, period=Day(20))
# skip first 5 days, as is done in the RainMaker challenge
RainMaker.skip!(rain_gauge, Day(5))
# evaluate rain gauge
lsc = rain_gauge.accumulated_rain_large_scale
conv = rain_gauge.accumulated_rain_convection
total_precip = maximum(lsc) + maximum(conv)
return rain_gauge, total_precip
end
final_params = [
1.85625,
2230.625,
3.0937500000000004,
-76.72500000000001,
6.187500000000014,
290.34375,
298.33125,
-3.46875,
2.6312500000000005,
32.534375,
]
rain_gauge, total_precip = max_precipitation(final_params)
total_precip
Chris Rackauckas: Evotrees.jl Surrogate, Claude generated
path: /submissions/chris_evotrees.jl
rank: 9. of 14 submissions
author = "Chris Rackauckas"
description = "Evotrees.jl Surrogate, Claude generated"
#=
using JLSO
using EvoTrees
using OptimizationBBO, Optimization
path = joinpath(dirname(@__DIR__), "data", "10kdata.jlso")
data = JLSO.load(path)
x = data[:d].inputs
y = data[:d].outputs
X_train = x'
y_train = vec(y)
config = EvoTreeRegressor(
nrounds=200,
max_depth=4,
eta=0.1,
rng=123
)
println("Training EvoTrees model...")
model = EvoTrees.fit_evotree(config; x_train=X_train, y_train=y_train)
println("Model trained successfully")
function surrogate_model(params)
pred = EvoTrees.predict(model, reshape(params, 1, :))
return Float64(pred[1])
end
function objective(params, p)
return -surrogate_model(params)
end
lower_bounds = [0., -2000., 0., -180., -90., 270., 270., -5., -5., 5.]
upper_bounds = [2., 5000., 30., 180., 90., 300., 300., 5., 5., 50.]
optf = OptimizationFunction(objective)
prob = Optimization.OptimizationProblem(optf, [1, 0, 1, -80, 40.45, 300, 273, 0, 0, 35],
lb=lower_bounds, ub=upper_bounds)
println("Starting optimization...")
sol = solve(prob, BBO_adaptive_de_rand_1_bin_radiuslimited(); maxtime=60)
println("Optimized parameters: ", sol.u)
println("Predicted max precipitation from surrogate: ", surrogate_model(sol.u))
=#
using SpeedyWeather, RainMaker
const PARAMETER_KEYS = (
:orography_scale, # [1], default: 1, scale of global orography
:mountain_height, # [m], default: 0, height of an additional azores mountain
:mountain_size, # [˚], default: 1, horizontal size of mountain in Pittsburgh
:mountain_lon, # [˚E], default: -80, longitude of that mountain in Pittsburgh
:mountain_lat, # [˚N], default: 40.45, latitude of that mountain in Pittsburgh
:temperature_equator, # [K], default: 300, sea surface temperature at the equator
:temperature_pole, # [K], default: 273, sea surfaec temperature at the poles
:temperature_usa, # [K], default: 0, land surface temperature anomaly over the USA
:temperature_pa, # [K], default: 0, land surface temperature anomaly in Pennsylvania
:zonal_wind, # [m/s], default: 35, zonal wind speed
)
const PARAMETER_DEFAULTS = [1, 0, 1, -80, 40.45, 300, 273, 0, 0, 35]
function max_precipitation(parameters::AbstractVector)
parameter_tuple = NamedTuple{PARAMETER_KEYS}(parameters)
return max_precipitation(parameter_tuple)
end
function max_precipitation(parameters::NamedTuple)
# define resolution. Use trunc=42, 63, 85, 127, ... for higher resolution, cubically slower
spectral_grid = SpectralGrid(trunc=31, nlayers=8)
# Define AquaPlanet ocean, for idealised sea surface temperatures
# but don't change land-sea mask = retain real ocean basins
ocean = AquaPlanet(spectral_grid,
temp_equator=parameters.temperature_equator,
temp_poles=parameters.temperature_pole)
land_temperature = ConstantLandTemperature(spectral_grid)
land = LandModel(spectral_grid; temperature=land_temperature)
initial_conditions = InitialConditions(
vordiv = ZonalWind(u₀=parameters.zonal_wind),
temp = JablonowskiTemperature(u₀=parameters.zonal_wind),
pres = PressureOnOrography(),
humid = ConstantRelativeHumidity())
orography = EarthOrography(spectral_grid, scale=parameters.orography_scale)
# construct model
model = PrimitiveWetModel(spectral_grid; ocean, land, initial_conditions, orography)
# Add rain gauge, locate in Pittsburgh PA
rain_gauge = RainGauge(spectral_grid, lond=-80, latd=40.45)
add!(model, rain_gauge)
# Initialize
simulation = initialize!(model, time=DateTime(2025, 7, 22))
# Add additional mountain
H = parameters.mountain_height
λ₀, φ₀, σ = parameters.mountain_lon, parameters.mountain_lat, parameters.mountain_size
set!(model, orography=(λ,φ) -> H*exp(-spherical_distance((λ,φ), (λ₀,φ₀), radius=360/2π)^2/2σ^2), add=true)
# land sea surface temperature anomalies
# 1. USA
set!(simulation, soil_temperature=
(λ, φ, k) -> (30 < φ < 50) && (240 < λ < 285) ? parameters.temperature_usa : 0, add=true)
# 2. Pennsylvania
A = parameters.temperature_pa
λ_az, φ_az, σ_az = -80, 40.45, 4 # location [˚], size [˚] of Azores
set!(simulation, soil_temperature=
(λ, φ, k) -> A*exp(-spherical_distance((λ,φ), (λ_az,φ_az), radius=360/2π)^2/2σ_az^2), add=true)
# Run simulation for 20 days
run!(simulation, period=Day(20))
# skip first 5 days, as is done in the RainMaker challenge
RainMaker.skip!(rain_gauge, Day(5))
# evaluate rain gauge
lsc = rain_gauge.accumulated_rain_large_scale
conv = rain_gauge.accumulated_rain_convection
total_precip = maximum(lsc) + maximum(conv)
return rain_gauge, total_precip
end
final_params = [1.99, -1494.66, 0.44, -175.36, -84.63, 274.32, 277.01, 4.70,
4.96, 27.44]
rain_gauge, total_precip = max_precipitation(final_params)
total_precip
Shirin Ermis: Aqua-planet simulation with a mountain
path: /submissions/aquaplanet_mountain.jl
rank: 10. of 14 submissions
author = "Shirin Ermis"
description = "Aqua-planet simulation with a mountain"
using SpeedyWeather, RainMaker
spectral_grid = SpectralGrid(trunc=31, nlayers=10)
model = PrimitiveWetModel(spectral_grid)
# Set up aqauaplanet but add large mountain in "North Sea" after initialization!
ocean = AquaPlanet(spectral_grid, temp_equator=302, temp_poles=300)
land_sea_mask = AquaPlanetMask(spectral_grid)
model = PrimitiveWetModel(spectral_grid; ocean, land_sea_mask)
# Add rain gauge
rain_gauge = RainGauge(spectral_grid, lond=-80, latd=40.45)
add!(model, rain_gauge)
# Initialize and run simulation
simulation = initialize!(model, time=DateTime(2000, 9, 1))
# Add mountain now! details for mountain
H, λ₀, φ₀, σ = 4000, 2, 51, 5 # height, lon, lat position, and width
set!(model, orography=(λ,φ) -> H*exp(-spherical_distance((λ,φ), (λ₀,φ₀), radius=360/2π)^2/2σ^2))
# Run simulation for 20 days
run!(simulation, period=Day(20))
# Plot the results (not needed for submission but doesn't hurt!)
using CairoMakie
# heatmap(model.orography.orography, title="My orogoraphy [m]") # check orography
RainMaker.plot(rain_gauge, rate_Δt=Hour(1))
Milan: Aquaplanet
path: /submissions/aquaplanet.jl
rank: 11. of 14 submissions
author = "Milan"
description = "Aquaplanet"
using SpeedyWeather, RainMaker
spectral_grid = SpectralGrid(trunc=31, nlayers=8)
# define aquaplanet
ocean = AquaPlanet(spectral_grid, temp_equator=302, temp_poles=273)
land_sea_mask = AquaPlanetMask(spectral_grid)
orography = NoOrography(spectral_grid)
model = PrimitiveWetModel(spectral_grid; ocean, land_sea_mask, orography)
rain_gauge = RainGauge(spectral_grid, lond=-80, latd=40.45)
add!(model, rain_gauge)
simulation = initialize!(model)
run!(simulation, period=Day(20))
Tim Reichelt: North Sea mountain
path: /submissions/north_sea_mountain.jl
rank: 12. of 14 submissions
author = "Tim Reichelt"
description = "North Sea mountain"
using SpeedyWeather, RainMaker
spectral_grid = SpectralGrid(trunc=31, nlayers=8)
model = PrimitiveWetModel(spectral_grid)
rain_gauge = RainGauge(spectral_grid, lond=-80, latd=40.45)
add!(model, rain_gauge)
simulation = initialize!(model)
# add a massive mountain at 51.75°N, 0°W, *after* model initialization
# using spherical_distance for geodesic distances, use radius=360/2π for distance in degrees
H, λ₀, φ₀, σ = 4000, 2, 51, 5 # height, lon, lat position, and width
set!(model, orography=(λ,φ) -> H*exp(-spherical_distance((λ,φ), (λ₀,φ₀), radius=360/2π)^2/2σ^2))
run!(simulation, period=Day(20))
Milan: default
path: /submissions/default.jl
rank: 13. of 14 submissions
author = "Milan"
description = "default"
using SpeedyWeather, RainMaker
spectral_grid = SpectralGrid(trunc=31, nlayers=8)
model = PrimitiveWetModel(spectral_grid)
rain_gauge = RainGauge(spectral_grid, lond=-80, latd=40.45)
add!(model, rain_gauge)
simulation = initialize!(model)
run!(simulation, period=Day(20))
Milan: Atlantic mountain
path: /submissions/atlantic_mountain.jl
rank: 14. of 14 submissions
author = "Milan"
description = "Atlantic mountain"
using SpeedyWeather, RainMaker
spectral_grid = SpectralGrid(trunc=31, nlayers=8)
model = PrimitiveWetModel(spectral_grid)
rain_gauge = RainGauge(spectral_grid, lond=-80, latd=40.45)
add!(model, rain_gauge)
simulation = initialize!(model)
# add a massive mountain at 50°N, 35°W, *after* model initialization
H, λ₀, φ₀, σ = 4000, 325, 50, 5 # height, lon, lat position, and width
set!(model, orography=(λ,φ) -> H*exp((-(λ-λ₀)^2 - (φ-φ₀)^2)/2σ^2), add=true)
run!(simulation, period=Day(20))