fixed plot and BC

docs/tutorials/standalone/Soil/layered_soil.jl

@@ -8,7 +8,7 @@
 # Note that all data used in this tutorial is available in their
 # publication.
 
-using Plots
+using CairoMakie
 import ClimaUtilities.SpaceVaryingInputs: SpaceVaryingInput
 import SciMLBase
 import ClimaTimeSteppers as CTS
@@ -30,7 +30,7 @@ dt = Float64(30);
 # Define the domain
 zmax = FT(0)
 zmin = FT(-1.1)
-nelems = 75
+nelems = 50
 Δ = FT((zmax - zmin) / nelems / 2)
 soil_domain = Column(; zlim = (zmin, zmax), nelements = nelems);
 
@@ -82,12 +82,9 @@ params = ClimaLand.Soil.RichardsParameters(;
 # From here on out, everything should look familiar if you've already gone
 # through the other soil tutorials.
 # Set Boundary conditions:
-# At the top, we use the observed value of Ksat at the top of the domain.
-# Setting the flux to be -Ksat is approximating the top as saturated.
-function top_flux_function(p, t)
-    return -0.0001033
-end
-top_bc = ClimaLand.Soil.WaterFluxBC(top_flux_function)
+# At the top, we set moisture to equivalent of 5cm hydraulic head
+# The bottom is free drainage
+top_bc = ClimaLand.Soil.MoistureStateBC((p, t) -> 0.46705)
 bottom_bc = ClimaLand.Soil.FreeDrainage()
 boundary_fluxes = (; top = top_bc, bottom = bottom_bc)
 soil = Soil.RichardsModel{FT}(;
@@ -100,8 +97,21 @@ soil = Soil.RichardsModel{FT}(;
 # Initial the state vectors, and set initial conditions
 Y, p, cds = initialize(soil);
 
-# Initial conditions
-Y.soil.ϑ_l .= 0.0353; # read from Figure 4 of Huang et al.
+# Initial conditions read from data
+data = readdlm("./docs/tutorials/standalone/Soil/sv_62_measurements.csv", '\t')
+data_z = Float64.(-1 .* data[2:end, 1]) ./ 100
+data_ic = Float64.(data[2:end, 2])
+data_8min = Float64.(data[2:end, 3])
+data_16min = Float64.(data[2:end, 4])
+data_24min = Float64.(data[2:end, 5])
+data_32min = Float64.(data[2:end, 6])
+data_40min = Float64.(data[2:end, 7])
+data_60min = Float64.(data[2:end, 8])
+Y.soil.ϑ_l .= SpaceVaryingInput(
+    reverse(data_z),
+    reverse(data_ic),
+    soil_domain.space.subsurface,
+)
 
 # We also set the initial conditions of the auxiliary state here:
 set_initial_cache! = make_set_initial_cache(soil)
@@ -139,28 +149,112 @@ prob = SciMLBase.ODEProblem(
 saveat = [0.0, 8.0, 16.0, 24.0, 32.0, 40.0, 60.0] .* 60 # chosen to compare with data in plots in paper
 sol = SciMLBase.solve(prob, ode_algo; dt = dt, saveat = saveat);
 
-z = parent(ClimaCore.Fields.coordinate_field(soil_domain.space.subsurface).z)
+z = parent(ClimaCore.Fields.coordinate_field(soil_domain.space.subsurface).z)[:]
 ϑ_l = [parent(sol.u[k].soil.ϑ_l) for k in 1:length(sol.t)]
-plot(ϑ_l[1], z, label = "initial", color = "grey", aspect_ratio = 0.8)
-plot!(ϑ_l[2], z, label = "8min", color = "orange")
-plot!(ϑ_l[3], z, label = "16min", color = "red")
-plot!(ϑ_l[4], z, label = "24min", color = "teal")
-plot!(ϑ_l[5], z, label = "32min", color = "blue")
-plot!(ϑ_l[6], z, label = "40min", color = "purple")
-plot!(ϑ_l[7], z, label = "60min", color = "green")
-scatter!(porosity, depth, label = "Porosity")
-plot!(legend = :bottomright)
-
-plot!(xlim = [0, 0.7])
-
-plot!(
-    ylim = [-1.1, 0],
+
+# Hydrus data
+
+hydrus = readdlm("./docs/tutorials/standalone/Soil/sv_62_hydrus.csv", '\t')
+hydrus_z = Float64.(-1 .* hydrus[:, 1]) ./ 100
+hydrus_ic = Float64.(hydrus[:, 2])
+hydrus_8min = Float64.(hydrus[:, 3])
+hydrus_16min = Float64.(hydrus[:, 4])
+hydrus_24min = Float64.(hydrus[:, 5])
+hydrus_32min = Float64.(hydrus[:, 6])
+hydrus_40min = Float64.(hydrus[:, 7])
+hydrus_60min = Float64.(hydrus[:, 8])
+
+
+fig = CairoMakie.Figure(size = (600, 800), fontsize = 30)
+ax = CairoMakie.Axis(
+    fig[1, 1],
+    xlabel = "Volumetric Water Content",
+    ylabel = "Depth (m)",
+    xgridvisible = false,
+    ygridvisible = false,
     yticks = [-1.1, -1, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1],
 )
+lines!(ax, ϑ_l[1][:], z, label = "ClimaLand", color = "black", linewidth=3)
+lines!(ax, ϑ_l[2][:], z, color = "blue", linewidth=3)
+lines!(ax, ϑ_l[3][:], z, color = "green", linewidth=3)
+lines!(ax, ϑ_l[4][:], z, color = "orange", linewidth=3)
+lines!(ax, ϑ_l[5][:], z, color = "purple", linewidth=3)
+lines!(ax, ϑ_l[6][:], z, color = "cyan", linewidth=3)
+lines!(ax, ϑ_l[7][:], z, color = "brown", linewidth=3)
 
-plot!(ylabel = "Depth (m)")
+lines!(
+    ax,
+    hydrus_ic,
+    hydrus_z,
+    label = "Hydrus",
+    color = "black",
+    linestyle = :dash,
+)
+lines!(ax, hydrus_8min, hydrus_z, color = "blue", linestyle = :dash, linewidth=3)
+lines!(ax, hydrus_16min, hydrus_z, color = "green", linestyle = :dash, linewidth=3)
+lines!(ax, hydrus_24min, hydrus_z, color = "orange", linestyle = :dash, linewidth=3)
+lines!(ax, hydrus_32min, hydrus_z, color = "purple", linestyle = :dash, linewidth=3)
+lines!(ax, hydrus_40min, hydrus_z, color = "cyan", linestyle = :dash, linewidth=3)
+lines!(ax, hydrus_60min, hydrus_z, color = "brown", linestyle = :dash, linewidth=3)
 
-plot!(xlabel = "Volumeteric Water Content")
+scatter!(
+    ax,
+    data_ic,
+    data_z,
+    label = "0 min",
+    color = "black",
+    marker = :xcross,
+)
+scatter!(
+    ax,
+    data_8min,
+    data_z,
+    color = "blue",
+    label = "8 min",
+    marker = :xcross,
+)
+scatter!(
+    ax,
+    data_16min,
+    data_z,
+    color = "green",
+    label = "16 min",
+    marker = :xcross,
+)
+scatter!(
+    ax,
+    data_24min,
+    data_z,
+    color = "orange",
+    label = "24 min",
+    marker = :xcross,
+)
+scatter!(
+    ax,
+    data_32min,
+    data_z,
+    color = "purple",
+    label = "32 min",
+    marker = :xcross,
+)
+scatter!(
+    ax,
+    data_40min,
+    data_z,
+    color = "cyan",
+    label = "40 min",
+    marker = :xcross,
+)
+scatter!(
+    ax,
+    data_60min,
+    data_z,
+    color = "brown",
+    label = "60 min",
+    marker = :xcross,
+)
+axislegend(ax, position = :rb, framevisible=false)
+limits!(ax, 0, 0.7, -1.1, 0)
 
-savefig("./sv62_alpha_2_inf_updated_data_climaland.png")
+CairoMakie.save("./sv62_alpha_2_inf_updated_data_climaland_moisturebc.png", fig)
 # ![](sv62_alpha_2_inf_updated_data_climaland.png)

src/standalone/Soil/boundary_conditions.jl

@@ -247,11 +247,18 @@ function boundary_flux!(
 
     # Lastly, we need an effective conductivity to compute the flux that results from
     # the gradient in pressure.
-    # currently we approximate this as equal to the center value at the top layer (K_c)
-    # More accurate would be to compute the mean between K_c and K evaluated at the boundary
-    # condition.
-    K_eff = ClimaLand.Domains.top_center_to_surface(p.soil.K)
-
+    # We use the arithmetic mean between K_c (at the first cell center) and K_bc evaluated at the boundary
+    # condition, below.
+
+    K_eff = p.soil.K_eff
+    K_sat_bc = ClimaLand.Domains.top_center_to_surface(model.parameters.K_sat)
+    @. K_eff = hydraulic_conductivity(
+        hcm_bc,
+        K_sat_bc,
+        max((θ_bc - θ_r_bc) / (ν_bc - θ_r_bc), 1),
+    ) # This does not take into account ice or temperature
+    K_center = ClimaLand.Domains.top_center_to_surface(p.soil.K)
+    @. K_eff = (K_eff + K_center) / 2
     # Pass in (ψ_bc .+ Δz) as to account for contribution of gravity (∂(ψ+z)/∂z
     @. bc_field = ClimaLand.diffusive_flux(
         K_eff,
@@ -378,8 +385,8 @@ function ClimaLand.set_dfluxBCdY!(
 )
     (; ν, hydrology_cm, S_s, θ_r) = model.parameters
 
-    # Copy center variables to top face space
-    KN = Domains.top_center_to_surface(p.soil.K)
+    # Copy center variables to top face space, except hydraulic conductivity
+    K_eff = p.soil.K_eff
     hydrology_cmN = Domains.top_center_to_surface(hydrology_cm)
     ϑ_lN = Domains.top_center_to_surface(Y.soil.ϑ_l)
     νN = Domains.top_center_to_surface(ν)
@@ -396,9 +403,10 @@ function ClimaLand.set_dfluxBCdY!(
 
     # Update dfluxBCdY at the top boundary in place
     # Calculate the value and convert it to a Covariant3Vector
+    # Ignore derivative of K_eff with respect to theta.
     @. p.soil.dfluxBCdY =
         covariant3_unit_vector(local_geometry_faceN) *
-        (KN * dψdϑ(hydrology_cmN, ϑ_lN, νN, θ_rN, S_sN) / Δz)
+        (K_eff * dψdϑ(hydrology_cmN, ϑ_lN, νN, θ_rN, S_sN) / Δz)
     return nothing
 end
 
@@ -852,7 +860,7 @@ top boundary.
 These variables are updated in place in `boundary_flux!`.
 """
 boundary_vars(bc::MoistureStateBC, ::ClimaLand.TopBoundary) =
-    (:top_bc, :top_bc_wvec, :dfluxBCdY)
+    (:top_bc, :K_eff, :top_bc_wvec, :dfluxBCdY)
 
 """
     boundary_var_domain_names(::MoistureStateBC, ::ClimaLand.TopBoundary)
@@ -861,7 +869,7 @@ An extension of the `boundary_var_domain_names` method for MoistureStateBC at th
 top boundary.
 """
 boundary_var_domain_names(bc::MoistureStateBC, ::ClimaLand.TopBoundary) =
-    (:surface, :surface, :surface)
+    (:surface, :surface, :surface, :surface)
 """
     boundary_var_types(::RichardsModel{FT},
                         ::MoistureStateBC,
@@ -872,10 +880,11 @@ An extension of the `boundary_var_types` method for MoistureStateBC at the
     top boundary.
 """
 boundary_var_types(
-    model::RichardsModel{FT},
+    model::AbstractSoilModel{FT},
     bc::MoistureStateBC,
     ::ClimaLand.TopBoundary,
 ) where {FT} = (
+    FT,
     FT,
     ClimaCore.Geometry.WVector{FT},
     ClimaCore.Geometry.Covariant3Vector{FT},