snowy land

.buildkite/longruns_gpu/pipeline.yml

@@ -31,6 +31,16 @@ steps:
   - group: "Global Land Models"
     steps:
 
+      - label: ":snow_capped_mountain: Snowy Land"
+        command:
+          - julia --color=yes --project=.buildkite experiments/long_runs/snowy_land.jl
+        artifact_paths: "snowy_land_longrun_gpu/*png"
+        agents:
+          slurm_gpus: 1
+          slurm_time: 03:30:00
+        env:
+          CLIMACOMMS_DEVICE: "CUDA"
+
       - label: ":seedling: Soil-Canopy"
         command:
           - julia --color=yes --project=.buildkite experiments/long_runs/land.jl
@@ -41,7 +51,7 @@ steps:
         env:
           CLIMACOMMS_DEVICE: "CUDA"
 
-      - label: ":sun: California regional simulation"
+      - label: ":sunglasses: California regional simulation"
         command:
           - julia --color=yes --project=.buildkite experiments/long_runs/land_region.jl
         artifact_paths: "california_longrun_gpu/*png"

NEWS.md

@@ -4,6 +4,9 @@ ClimaLand.jl Release Notes
 main
 --------
 
+- Integrated land model with snow, soil, canopy
+  PR[#863](https://github.com/CliMA/ClimaLand.jl/pull/834)
+
 v0.15.4
 -------
 - The foliage clumping index of the radiative transfer model now varies spatially, using MODIS data

experiments/integrated/fluxnet/US-Ha1/US-Ha1_simulation.jl

@@ -17,7 +17,7 @@ h_stem = FT(14) # m
 # TIME STEPPING:
 
 # Starting time:
-t0 = Float64(120 * 3600 * 24)# start mid year to avoid snow
+t0 = Float64(0)
 
 # Time step size:
 dt = Float64(450)

experiments/integrated/fluxnet/US-MOz/US-MOz_simulation.jl

@@ -17,7 +17,7 @@ h_leaf = FT(9.5) # m
 # TIME STEPPING:
 
 # Starting time:
-t0 = Float64(120 * 3600 * 24)# start mid year to avoid snow
+t0 = Float64(0)
 
 # Time step size:
-dt = Float64(900)
+dt = Float64(450)

experiments/integrated/fluxnet/US-NR1/US-NR1_parameters.jl

@@ -20,11 +20,11 @@ atmos_h = FT(21.5)
 
 # Soil parameters
 soil_ν = FT(0.45) # m3/m3
-soil_K_sat = FT(4e-7) # m/s, matches Natan
+soil_K_sat = FT(4e-7) # m/s,
 soil_S_s = FT(1e-3) # 1/m, guess
 soil_vg_n = FT(2.05) # unitless
 soil_vg_α = FT(0.04) # inverse meters
-θ_r = FT(0.067) # m3/m3, from Wang et al. 2021 https://doi.org/10.5194/gmd-14-6741-2021
+θ_r = FT(0.0) # m3/m3, 
 
 # Soil makeup
 ν_ss_quartz = FT(0.1)

experiments/integrated/fluxnet/US-NR1/US-NR1_simulation.jl

@@ -17,7 +17,7 @@ h_stem = FT(7.5) # m
 # TIME STEPPING:
 
 # Starting time:
-t0 = Float64(120 * 3600 * 24)# start mid year to avoid snow
+t0 = Float64(0)
 
 # Time step size:
-dt = Float64(450)
+dt = Float64(200)

experiments/integrated/fluxnet/US-Var/US-Var_simulation.jl

@@ -17,7 +17,7 @@ h_stem = FT(0) # m
 # TIME STEPPING:
 
 # Starting time:
-t0 = Float64(21 * 3600 * 24)# start day 21 of the year
+t0 = Float64(0)
 
 # Time step size:
 dt = Float64(900)

experiments/integrated/fluxnet/fluxnet_simulation.jl

@@ -2,8 +2,9 @@
 fluxtower sites. These work in tandem with the site-specific timing parameters
 found in the {site-ID}_simulation.jl files in each site directory."""
 
-N_spinup_days = 30
-N_days = N_spinup_days + 30
+N_spinup_days = 15
+N_days = N_spinup_days + 340
+
 tf = Float64(t0 + 3600 * 24 * N_days)
 t_spinup = Float64(t0 + N_spinup_days * 3600 * 24)
 

experiments/integrated/fluxnet/run_fluxnet.jl

@@ -12,6 +12,7 @@ using StatsBase
 
 using ClimaLand
 using ClimaLand.Domains: Column
+using ClimaLand.Snow
 using ClimaLand.Soil
 using ClimaLand.Soil.Biogeochemistry
 using ClimaLand.Canopy
@@ -186,14 +187,19 @@ shared_params = SharedCanopyParameters{FT, typeof(earth_param_set)}(
 )
 
 canopy_model_args = (; parameters = shared_params, domain = canopy_domain)
+
+# Snow model
+snow_parameters = SnowParameters{FT}(dt; earth_param_set = earth_param_set);
+snow_args = (; parameters = snow_parameters, domain = canopy_domain);
+snow_model_type = Snow.SnowModel
 # Integrated plant hydraulics and soil model
 land_input = (
     atmos = atmos,
     radiation = radiation,
     soil_organic_carbon = Csom,
     runoff = ClimaLand.Soil.Runoff.SurfaceRunoff(),
 )
-land = SoilCanopyModel{FT}(;
+land = LandModel{FT}(;
     soilco2_type = soilco2_type,
     soilco2_args = soilco2_args,
     land_args = land_input,
@@ -202,6 +208,8 @@ land = SoilCanopyModel{FT}(;
     canopy_component_types = canopy_component_types,
     canopy_component_args = canopy_component_args,
     canopy_model_args = canopy_model_args,
+    snow_args = snow_args,
+    snow_model_type = snow_model_type,
 )
 
 Y, p, cds = initialize(land)
@@ -240,6 +248,10 @@ for i in 1:(n_stem + n_leaf)
 end
 
 Y.canopy.energy.T = drivers.TA.values[1 + Int(round(t0 / DATA_DT))] # Get atmos temperature at t0
+
+Y.snow.S .= 0.0
+Y.snow.U .= 0.0
+
 set_initial_cache! = make_set_initial_cache(land)
 set_initial_cache!(p, Y, t0);
 
@@ -309,6 +321,7 @@ short_names_1D = [
     "lhf", # LHF
     "ghf", # G
     "rn", # Rn
+    "swe",
 ]
 short_names_2D = [
     "swc", # swc_sfc or swc_5 or swc_10
@@ -323,11 +336,11 @@ hourly_diag_name_2D = short_names_2D .* "_1h_average"
 # diagnostic_as_vectors()[2] is a vector of a variable,
 # whereas diagnostic_as_vectors()[1] is a vector or time associated with that variable.
 # We index to only extract the period post-spinup.
-SIF, AR, g_stomata, GPP, canopy_T, SW_u, LW_u, ER, ET, β, SHF, LHF, G, Rn = [
-    ClimaLand.Diagnostics.diagnostic_as_vectors(d_writer, diag_name)[2][(N_spinup_days * 24):end]
-    for diag_name in hourly_diag_name
-]
-
+SIF, AR, g_stomata, GPP, canopy_T, SW_u, LW_u, ER, ET, β, SHF, LHF, G, Rn, SWE =
+    [
+        ClimaLand.Diagnostics.diagnostic_as_vectors(d_writer, diag_name)[2][(N_spinup_days * 24):end]
+        for diag_name in hourly_diag_name
+    ]
 
 swc, soil_T, si = [
     ClimaLand.Diagnostics.diagnostic_as_vectors(
@@ -510,7 +523,7 @@ if isfile(
     )
 end
 
-# Water content in soil
+# Water content in soil and snow
 # Soil water content
 # Current resolution has the first layer at 0.1 cm, the second at 5cm.
 plt1 = Plots.plot(size = (1500, 800))
@@ -546,7 +559,23 @@ if drivers.SWC.status != absent
         label = "Data",
     )
 end
-plt2 = Plots.plot(
+plt2 = Plots.plot(size = (1500, 800))
+Plots.plot!(
+    plt2,
+    model_times ./ 3600 ./ 24,
+    SWE,
+    label = "Model",
+    xlim = [
+        minimum(model_times ./ 3600 ./ 24),
+        maximum(model_times ./ 3600 ./ 24),
+    ],
+    ylim = [0.0, 0.15],
+    xlabel = "Days",
+    ylabel = "SWE [m]",
+    color = "blue",
+    margin = 10Plots.mm,
+)
+plt3 = Plots.plot(
     data_times ./ 3600 ./ 24,
     (drivers.P.values .* (-1e3 * 24 * 3600) .* (1 .- snow_frac))[data_id_post_spinup],
     label = "Rain (data)",
@@ -560,13 +589,13 @@ plt2 = Plots.plot(
     size = (1500, 400),
 )
 Plots.plot!(
-    plt2,
+    plt3,
     data_times ./ 3600 ./ 24,
     (drivers.P.values .* (-1e3 * 24 * 3600) .* snow_frac)[data_id_post_spinup],
     label = "Snow (data)",
     ylabel = "Precipitation [mm/day]",
 )
-Plots.plot(plt2, plt1, layout = grid(2, 1, heights = [0.3, 0.7]))
+Plots.plot(plt3, plt2, plt1, layout = grid(3, 1, heights = [0.2, 0.4, 0.4]))
 Plots.savefig(joinpath(savedir, "ground_water_content.png"))
 
 plt1 = Plots.plot(size = (1500, 800))