# Soil heat with SFCC and constant BCs

In this example, we use the preset `SoilHeatTile`

to construct a `Tile`

consisting of a soil column with heat conduction and zero-flux boundary conditions. This is a useful test case for checking energy conservation since we can guarantee that no energy is being added or removed at the boundaries.

`using CryoGrid`

Select default grid and initial temperature profile.

```
grid = CryoGrid.Presets.DefaultGrid_10cm
tempprofile = TemperatureProfile(
0.0u"m" => -20.0u"°C",
1000.0u"m" => 10.0u"°C",
);
```

Here we use a simple single layer model with default soil parameters (50% porosity, no organic).

```
soilprofile = SoilProfile(
0.0u"m" => MineralOrganic()
);
```

Here we specify the soil freezing characteristic curve (SFCC) formulation of Painter and Karra (2014). The van Genuchten parameters `α=0.1`

and `n=1.8`

correspond to a silty soil.

```
sfcc = PainterKarra(swrc=VanGenuchten(α=1.0, n=2.0));
import Plots
Plots.plot(-2.0u"°C":0.01u"K":0.0u"°C", sfcc)
```

Enthalpy form of the heat transfer operator (i.e. prognostic :H). In this case, this is equivalent to the shorthand `SoilHeatTile(:H, ...)`

. However, it's worth demonstrating how the operator can be explicitly specified.

```
heatop = Heat.Diffusion1D(:H)
initT = initializer(:T, tempprofile)
tile = CryoGrid.Presets.SoilHeatTile(
heatop,
# 10 W/m^2 in and out, i.e. net zero flux
ConstantBC(HeatBalance, CryoGrid.Neumann, 10.0u"W/m^2"),
ConstantBC(HeatBalance, CryoGrid.Neumann, 10.0u"W/m^2"),
soilprofile,
initT;
grid,
freezecurve=sfcc,
);
```

Define the simulation time span.

```
tspan = (DateTime(2010,1,1),DateTime(2010,12,31))
u0, du0 = initialcondition!(tile, tspan);
```

Construct and solve the `CryoGridProblem`

using the built-in forward Euler integrator.

```
prob = CryoGridProblem(tile, u0, tspan, saveat=3600.0, savevars=(:T,:H))
sol = @time solve(prob, dtmax=900.0);
out = CryoGridOutput(sol)
```

```
CryoGridOutput with 2913 time steps (2010-01-01T00:00:00 to 2010-12-31T00:00:00) and 2 variables:
H => DimArray of Quantity{Float64, 𝐌 𝐋^-1 𝐓^-2, Unitful.FreeUnits{(J, m^-3), 𝐌 𝐋^-1 𝐓^-2, nothing}} with dimensions (178, 2913)
T => DimArray of Quantity{Float64, 𝚯, Unitful.FreeUnits{(K,), 𝚯, Unitful.Affine{-5463//20}}} with dimensions (178, 2913)
```

Compute total energy balance error.

```
Htot = Diagnostics.integrate(out.H, grid)
mass_balance_error = Htot[end] - Htot[1]
```

`0.0006561279296875 J`

Plot it!

```
import Plots
zs = [1,45,95,195,495,795,995]u"cm"
Diagnostics.plot_at_depths(:T, out, zs, ylabel="Temperature", leg=false, size=(800,500), dpi=150)
```

Plot the energy balance error over time.

`Plots.plot(uconvert.(u"MJ", Htot .- Htot[1]), title="", ylabel="Energy balance error")`

*This page was generated using Literate.jl.*