diff --git a/Manifest.toml b/Manifest.toml
index ff43819..5752fb0 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -2,7 +2,7 @@
 
 julia_version = "1.9.2"
 manifest_format = "2.0"
-project_hash = "919dab67366ad3399e882328af9dc0cc9ff8b04f"
+project_hash = "c88ef8733706dadcee8b23c8aa1dbaecea03d700"
 
 [[deps.AbstractFFTs]]
 deps = ["LinearAlgebra"]
@@ -127,6 +127,12 @@ git-tree-sha1 = "0eee5eb66b1cf62cd6ad1b460238e60e4b09400c"
 uuid = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
 version = "0.2.4"
 
+[[deps.CommonSubexpressions]]
+deps = ["MacroTools", "Test"]
+git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7"
+uuid = "bbf7d656-a473-5ed7-a52c-81e309532950"
+version = "0.3.0"
+
 [[deps.Compat]]
 deps = ["UUIDs"]
 git-tree-sha1 = "e460f044ca8b99be31d35fe54fc33a5c33dd8ed7"
@@ -187,6 +193,28 @@ version = "1.0.0"
 deps = ["Printf"]
 uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
 
+[[deps.DiffResults]]
+deps = ["StaticArraysCore"]
+git-tree-sha1 = "782dd5f4561f5d267313f23853baaaa4c52ea621"
+uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
+version = "1.1.0"
+
+[[deps.DiffRules]]
+deps = ["IrrationalConstants", "LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"]
+git-tree-sha1 = "23163d55f885173722d1e4cf0f6110cdbaf7e272"
+uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
+version = "1.15.1"
+
+[[deps.Distances]]
+deps = ["LinearAlgebra", "Statistics", "StatsAPI"]
+git-tree-sha1 = "b6def76ffad15143924a2199f72a5cd883a2e8a9"
+uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
+version = "0.10.9"
+weakdeps = ["SparseArrays"]
+
+    [deps.Distances.extensions]
+    DistancesSparseArraysExt = "SparseArrays"
+
 [[deps.Distributed]]
 deps = ["Random", "Serialization", "Sockets"]
 uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
@@ -233,6 +261,16 @@ version = "1.16.1"
 [[deps.FileWatching]]
 uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"
 
+[[deps.ForwardDiff]]
+deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions"]
+git-tree-sha1 = "cf0fe81336da9fb90944683b8c41984b08793dad"
+uuid = "f6369f11-7733-5829-9624-2563aa707210"
+version = "0.10.36"
+weakdeps = ["StaticArrays"]
+
+    [deps.ForwardDiff.extensions]
+    ForwardDiffStaticArraysExt = "StaticArrays"
+
 [[deps.Future]]
 deps = ["Random"]
 uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"
@@ -311,9 +349,9 @@ version = "0.4.33"
 
 [[deps.JLLWrappers]]
 deps = ["Artifacts", "Preferences"]
-git-tree-sha1 = "a7e91ef94114d5bc8952bcaa8d6ff952cf709808"
+git-tree-sha1 = "7e5d6779a1e09a36db2a7b6cff50942a0a7d0fca"
 uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
-version = "1.4.2"
+version = "1.5.0"
 
 [[deps.JSON3]]
 deps = ["Dates", "Mmap", "Parsers", "PrecompileTools", "StructTypes", "UUIDs"]
@@ -389,9 +427,9 @@ uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
 [[deps.LogExpFunctions]]
 deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"]
-git-tree-sha1 = "c3ce8e7420b3a6e071e0fe4745f5d4300e37b13f"
+git-tree-sha1 = "7d6dd4e9212aebaeed356de34ccf262a3cd415aa"
 uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
-version = "0.3.24"
+version = "0.3.26"
 
     [deps.LogExpFunctions.extensions]
     LogExpFunctionsChainRulesCoreExt = "ChainRulesCore"
@@ -446,9 +484,9 @@ version = "5.3.1+0"
 
 [[deps.MacroTools]]
 deps = ["Markdown", "Random"]
-git-tree-sha1 = "42324d08725e200c23d4dfb549e0d5d89dede2d2"
+git-tree-sha1 = "9ee1618cbf5240e6d4e0371d6f24065083f60c48"
 uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
-version = "0.5.10"
+version = "0.5.11"
 
 [[deps.Markdown]]
 deps = ["Base64"]
@@ -478,6 +516,12 @@ git-tree-sha1 = "4263c4220f22e20729329838bf7e94a49d1ac32f"
 uuid = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
 version = "0.12.17"
 
+[[deps.NaNMath]]
+deps = ["OpenLibm_jll"]
+git-tree-sha1 = "0877504529a3e5c3343c6f8b4c0381e57e4387e4"
+uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
+version = "1.0.2"
+
 [[deps.NetCDF_jll]]
 deps = ["Artifacts", "Bzip2_jll", "HDF5_jll", "JLLWrappers", "LibCURL_jll", "Libdl", "XML2_jll", "Zlib_jll", "Zstd_jll"]
 git-tree-sha1 = "10c612c81eaffdd6b7c28a45a554cdd9d2f40ff1"
@@ -489,10 +533,10 @@ uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
 version = "1.2.0"
 
 [[deps.Oceananigans]]
-deps = ["Adapt", "CUDA", "Crayons", "CubedSphere", "Dates", "DocStringExtensions", "FFTW", "Glob", "IncompleteLU", "InteractiveUtils", "IterativeSolvers", "JLD2", "KernelAbstractions", "LinearAlgebra", "Logging", "MPI", "NCDatasets", "OffsetArrays", "OrderedCollections", "PencilArrays", "PencilFFTs", "Pkg", "Printf", "Random", "Rotations", "SeawaterPolynomials", "SparseArrays", "Statistics", "StructArrays"]
-git-tree-sha1 = "f4bebf2d3e79f81fe516afb3fd71ec829cd8e36f"
+deps = ["Adapt", "CUDA", "Crayons", "CubedSphere", "Dates", "Distances", "DocStringExtensions", "FFTW", "Glob", "IncompleteLU", "InteractiveUtils", "IterativeSolvers", "JLD2", "KernelAbstractions", "LinearAlgebra", "Logging", "MPI", "NCDatasets", "OffsetArrays", "OrderedCollections", "PencilArrays", "PencilFFTs", "Pkg", "Printf", "Random", "Rotations", "SeawaterPolynomials", "SparseArrays", "Statistics", "StructArrays"]
+git-tree-sha1 = "996eb159d0a1380d09beb7190ed0807d51ba0924"
 uuid = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"
-version = "0.86.0"
+version = "0.87.1"
 
 [[deps.OffsetArrays]]
 deps = ["Adapt"]
@@ -572,9 +616,9 @@ version = "0.3.2"
 
 [[deps.PrecompileTools]]
 deps = ["Preferences"]
-git-tree-sha1 = "9673d39decc5feece56ef3940e5dafba15ba0f81"
+git-tree-sha1 = "03b4c25b43cb84cee5c90aa9b5ea0a78fd848d2f"
 uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
-version = "1.1.2"
+version = "1.2.0"
 
 [[deps.Preferences]]
 deps = ["TOML"]
@@ -588,9 +632,9 @@ uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 
 [[deps.ProgressBars]]
 deps = ["Printf"]
-git-tree-sha1 = "9d84c8646109eb8bc7a006d59b157c64d5155c81"
+git-tree-sha1 = "b437cdb0385ed38312d91d9c00c20f3798b30256"
 uuid = "49802e3a-d2f1-5c88-81d8-b72133a6f568"
-version = "1.5.0"
+version = "1.5.1"
 
 [[deps.Quaternions]]
 deps = ["LinearAlgebra", "Random", "RealDot"]
@@ -641,6 +685,12 @@ git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
 uuid = "ae029012-a4dd-5104-9daa-d747884805df"
 version = "1.3.0"
 
+[[deps.RootSolvers]]
+deps = ["DocStringExtensions", "ForwardDiff"]
+git-tree-sha1 = "9fb3462240d2898be5d4acf8925e47f70ec64d07"
+uuid = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74"
+version = "0.3.5"
+
 [[deps.Roots]]
 deps = ["ChainRulesCore", "CommonSolve", "Printf", "Setfield"]
 git-tree-sha1 = "ff42754a57bb0d6dcfe302fd0d4272853190421f"
@@ -752,6 +802,12 @@ deps = ["LinearAlgebra", "SparseArrays"]
 uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 version = "1.9.0"
 
+[[deps.StatsAPI]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "45a7769a04a3cf80da1c1c7c60caf932e6f4c9f7"
+uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
+version = "1.6.0"
+
 [[deps.Strided]]
 deps = ["LinearAlgebra", "StridedViews", "TupleTools"]
 git-tree-sha1 = "137303f5e0a39f966b462c53ae2c5c6e34c4828b"
diff --git a/Project.toml b/Project.toml
index 633b53c..c6f556f 100644
--- a/Project.toml
+++ b/Project.toml
@@ -6,10 +6,10 @@ version = "0.1.0"
 [deps]
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 Oceananigans = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"
+RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 
 [compat]
-Oceananigans = "0.86"
 KernelAbstractions = "0.9"
 Roots = "2"
 
diff --git a/docs/Manifest.toml b/docs/Manifest.toml
index 03cb182..56447ea 100644
--- a/docs/Manifest.toml
+++ b/docs/Manifest.toml
@@ -2,16 +2,97 @@
 
 julia_version = "1.9.2"
 manifest_format = "2.0"
-project_hash = "c8b9437d4b095fc9a3fcf7a39d8c1fad8f05210e"
+project_hash = "c17176d00ba2b8c3e38a54d717165bc946d88851"
 
 [[deps.ANSIColoredPrinters]]
 git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c"
 uuid = "a4c015fc-c6ff-483c-b24f-f7ea428134e9"
 version = "0.0.1"
 
+[[deps.AbstractFFTs]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "d92ad398961a3ed262d8bf04a1a2b8340f915fef"
+uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
+version = "1.5.0"
+weakdeps = ["ChainRulesCore", "Test"]
+
+    [deps.AbstractFFTs.extensions]
+    AbstractFFTsChainRulesCoreExt = "ChainRulesCore"
+    AbstractFFTsTestExt = "Test"
+
+[[deps.AbstractLattices]]
+git-tree-sha1 = "f35684b7349da49fcc8a9e520e30e45dbb077166"
+uuid = "398f06c4-4d28-53ec-89ca-5b2656b7603d"
+version = "0.2.1"
+
+[[deps.AbstractTrees]]
+git-tree-sha1 = "faa260e4cb5aba097a73fab382dd4b5819d8ec8c"
+uuid = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
+version = "0.4.4"
+
+[[deps.Adapt]]
+deps = ["LinearAlgebra", "Requires"]
+git-tree-sha1 = "76289dc51920fdc6e0013c872ba9551d54961c24"
+uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
+version = "3.6.2"
+weakdeps = ["StaticArrays"]
+
+    [deps.Adapt.extensions]
+    AdaptStaticArraysExt = "StaticArrays"
+
+[[deps.Animations]]
+deps = ["Colors"]
+git-tree-sha1 = "e81c509d2c8e49592413bfb0bb3b08150056c79d"
+uuid = "27a7e980-b3e6-11e9-2bcd-0b925532e340"
+version = "0.4.1"
+
+[[deps.ArgTools]]
+uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
+version = "1.1.1"
+
+[[deps.ArrayInterface]]
+deps = ["Adapt", "LinearAlgebra", "Requires", "SparseArrays", "SuiteSparse"]
+git-tree-sha1 = "f83ec24f76d4c8f525099b2ac475fc098138ec31"
+uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
+version = "7.4.11"
+
+    [deps.ArrayInterface.extensions]
+    ArrayInterfaceBandedMatricesExt = "BandedMatrices"
+    ArrayInterfaceBlockBandedMatricesExt = "BlockBandedMatrices"
+    ArrayInterfaceCUDAExt = "CUDA"
+    ArrayInterfaceGPUArraysCoreExt = "GPUArraysCore"
+    ArrayInterfaceStaticArraysCoreExt = "StaticArraysCore"
+    ArrayInterfaceTrackerExt = "Tracker"
+
+    [deps.ArrayInterface.weakdeps]
+    BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
+    BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
+    CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+    GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
+    StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
+    Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c"
+
 [[deps.Artifacts]]
 uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
 
+[[deps.Automa]]
+deps = ["TranscodingStreams"]
+git-tree-sha1 = "ef9997b3d5547c48b41c7bd8899e812a917b409d"
+uuid = "67c07d97-cdcb-5c2c-af73-a7f9c32a568b"
+version = "0.8.4"
+
+[[deps.AxisAlgorithms]]
+deps = ["LinearAlgebra", "Random", "SparseArrays", "WoodburyMatrices"]
+git-tree-sha1 = "66771c8d21c8ff5e3a93379480a2307ac36863f7"
+uuid = "13072b0f-2c55-5437-9ae7-d433b7a33950"
+version = "1.0.1"
+
+[[deps.AxisArrays]]
+deps = ["Dates", "IntervalSets", "IterTools", "RangeArrays"]
+git-tree-sha1 = "16351be62963a67ac4083f748fdb3cca58bfd52f"
+uuid = "39de3d68-74b9-583c-8d2d-e117c070f3a9"
+version = "0.4.7"
+
 [[deps.Base64]]
 uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 
@@ -20,28 +101,212 @@ git-tree-sha1 = "43b1a4a8f797c1cddadf60499a8a077d4af2cd2d"
 uuid = "d1d4a3ce-64b1-5f1a-9ba4-7e7e69966f35"
 version = "0.1.7"
 
+[[deps.Bzip2_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "19a35467a82e236ff51bc17a3a44b69ef35185a2"
+uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
+version = "1.0.8+0"
+
+[[deps.CEnum]]
+git-tree-sha1 = "eb4cb44a499229b3b8426dcfb5dd85333951ff90"
+uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82"
+version = "0.4.2"
+
+[[deps.CRC32c]]
+uuid = "8bf52ea8-c179-5cab-976a-9e18b702a9bc"
+
+[[deps.CRlibm]]
+deps = ["CRlibm_jll"]
+git-tree-sha1 = "32abd86e3c2025db5172aa182b982debed519834"
+uuid = "96374032-68de-5a5b-8d9e-752f78720389"
+version = "1.0.1"
+
+[[deps.CRlibm_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "e329286945d0cfc04456972ea732551869af1cfc"
+uuid = "4e9b3aee-d8a1-5a3d-ad8b-7d824db253f0"
+version = "1.0.1+0"
+
+[[deps.Cairo]]
+deps = ["Cairo_jll", "Colors", "Glib_jll", "Graphics", "Libdl", "Pango_jll"]
+git-tree-sha1 = "d0b3f8b4ad16cb0a2988c6788646a5e6a17b6b1b"
+uuid = "159f3aea-2a34-519c-b102-8c37f9878175"
+version = "1.0.5"
+
+[[deps.CairoMakie]]
+deps = ["Base64", "Cairo", "Colors", "FFTW", "FileIO", "FreeType", "GeometryBasics", "LinearAlgebra", "Makie", "PrecompileTools", "SHA"]
+git-tree-sha1 = "30562a68ded3dabe80109caf6b4de73a48ac27bc"
+uuid = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+version = "0.10.8"
+
+[[deps.Cairo_jll]]
+deps = ["Artifacts", "Bzip2_jll", "CompilerSupportLibraries_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"]
+git-tree-sha1 = "4b859a208b2397a7a623a03449e4636bdb17bcf2"
+uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a"
+version = "1.16.1+1"
+
+[[deps.Calculus]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "f641eb0a4f00c343bbc32346e1217b86f3ce9dad"
+uuid = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
+version = "0.5.1"
+
+[[deps.ChainRulesCore]]
+deps = ["Compat", "LinearAlgebra", "SparseArrays"]
+git-tree-sha1 = "e30f2f4e20f7f186dc36529910beaedc60cfa644"
+uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+version = "1.16.0"
+
 [[deps.CodecZlib]]
 deps = ["TranscodingStreams", "Zlib_jll"]
 git-tree-sha1 = "02aa26a4cf76381be7f66e020a3eddeb27b0a092"
 uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
 version = "0.7.2"
 
+[[deps.ColorBrewer]]
+deps = ["Colors", "JSON", "Test"]
+git-tree-sha1 = "61c5334f33d91e570e1d0c3eb5465835242582c4"
+uuid = "a2cac450-b92f-5266-8821-25eda20663c8"
+version = "0.4.0"
+
+[[deps.ColorSchemes]]
+deps = ["ColorTypes", "ColorVectorSpace", "Colors", "FixedPointNumbers", "PrecompileTools", "Random"]
+git-tree-sha1 = "d9a8f86737b665e15a9641ecbac64deef9ce6724"
+uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
+version = "3.23.0"
+
+[[deps.ColorTypes]]
+deps = ["FixedPointNumbers", "Random"]
+git-tree-sha1 = "eb7f0f8307f71fac7c606984ea5fb2817275d6e4"
+uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
+version = "0.11.4"
+
+[[deps.ColorVectorSpace]]
+deps = ["ColorTypes", "FixedPointNumbers", "LinearAlgebra", "SpecialFunctions", "Statistics", "TensorCore"]
+git-tree-sha1 = "600cc5508d66b78aae350f7accdb58763ac18589"
+uuid = "c3611d14-8923-5661-9e6a-0046d554d3a4"
+version = "0.9.10"
+
+[[deps.Colors]]
+deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
+git-tree-sha1 = "fc08e5930ee9a4e03f84bfb5211cb54e7769758a"
+uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
+version = "0.12.10"
+
+[[deps.Combinatorics]]
+git-tree-sha1 = "08c8b6831dc00bfea825826be0bc8336fc369860"
+uuid = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
+version = "1.0.2"
+
+[[deps.CommonSubexpressions]]
+deps = ["MacroTools", "Test"]
+git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7"
+uuid = "bbf7d656-a473-5ed7-a52c-81e309532950"
+version = "0.3.0"
+
+[[deps.Compat]]
+deps = ["UUIDs"]
+git-tree-sha1 = "e460f044ca8b99be31d35fe54fc33a5c33dd8ed7"
+uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
+version = "4.9.0"
+weakdeps = ["Dates", "LinearAlgebra"]
+
+    [deps.Compat.extensions]
+    CompatLinearAlgebraExt = "LinearAlgebra"
+
+[[deps.CompilerSupportLibraries_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
+version = "1.0.5+0"
+
 [[deps.ConcurrentUtilities]]
 deps = ["Serialization", "Sockets"]
 git-tree-sha1 = "5372dbbf8f0bdb8c700db5367132925c0771ef7e"
 uuid = "f0e56b4a-5159-44fe-b623-3e5288b988bb"
 version = "2.2.1"
 
+[[deps.ConstructionBase]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "fe2838a593b5f776e1597e086dcd47560d94e816"
+uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
+version = "1.5.3"
+weakdeps = ["IntervalSets", "StaticArrays"]
+
+    [deps.ConstructionBase.extensions]
+    ConstructionBaseIntervalSetsExt = "IntervalSets"
+    ConstructionBaseStaticArraysExt = "StaticArrays"
+
+[[deps.Contour]]
+git-tree-sha1 = "d05d9e7b7aedff4e5b51a029dced05cfb6125781"
+uuid = "d38c429a-6771-53c6-b99e-75d170b6e991"
+version = "0.6.2"
+
+[[deps.DataAPI]]
+git-tree-sha1 = "8da84edb865b0b5b0100c0666a9bc9a0b71c553c"
+uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
+version = "1.15.0"
+
 [[deps.DataDeps]]
 deps = ["HTTP", "Libdl", "Reexport", "SHA", "p7zip_jll"]
 git-tree-sha1 = "6e8d74545d34528c30ccd3fa0f3c00f8ed49584c"
 uuid = "124859b0-ceae-595e-8997-d05f6a7a8dfe"
 version = "0.7.11"
 
+[[deps.DataStructures]]
+deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
+git-tree-sha1 = "3dbd312d370723b6bb43ba9d02fc36abade4518d"
+uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+version = "0.18.15"
+
+[[deps.DataValueInterfaces]]
+git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
+uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
+version = "1.0.0"
+
 [[deps.Dates]]
 deps = ["Printf"]
 uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
 
+[[deps.DelaunayTriangulation]]
+deps = ["DataStructures", "EnumX", "ExactPredicates", "Random", "SimpleGraphs"]
+git-tree-sha1 = "a1d8532de83f8ce964235eff1edeff9581144d02"
+uuid = "927a84f5-c5f4-47a5-9785-b46e178433df"
+version = "0.7.2"
+weakdeps = ["MakieCore"]
+
+    [deps.DelaunayTriangulation.extensions]
+    DelaunayTriangulationMakieCoreExt = "MakieCore"
+
+[[deps.DiffResults]]
+deps = ["StaticArraysCore"]
+git-tree-sha1 = "782dd5f4561f5d267313f23853baaaa4c52ea621"
+uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
+version = "1.1.0"
+
+[[deps.DiffRules]]
+deps = ["IrrationalConstants", "LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"]
+git-tree-sha1 = "23163d55f885173722d1e4cf0f6110cdbaf7e272"
+uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
+version = "1.15.1"
+
+[[deps.Distributed]]
+deps = ["Random", "Serialization", "Sockets"]
+uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
+
+[[deps.Distributions]]
+deps = ["FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SpecialFunctions", "Statistics", "StatsAPI", "StatsBase", "StatsFuns", "Test"]
+git-tree-sha1 = "938fe2981db009f531b6332e31c58e9584a2f9bd"
+uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
+version = "0.25.100"
+
+    [deps.Distributions.extensions]
+    DistributionsChainRulesCoreExt = "ChainRulesCore"
+    DistributionsDensityInterfaceExt = "DensityInterface"
+
+    [deps.Distributions.weakdeps]
+    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+    DensityInterface = "b429d917-457f-4dbc-8f4c-0cc954292b1d"
+
 [[deps.DocStringExtensions]]
 deps = ["LibGit2"]
 git-tree-sha1 = "2fb1e02f2b635d0845df5d7c167fec4dd739b00d"
@@ -54,38 +319,369 @@ git-tree-sha1 = "39fd748a73dce4c05a9655475e437170d8fb1b67"
 uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 version = "0.27.25"
 
+[[deps.Downloads]]
+deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"]
+uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
+version = "1.6.0"
+
+[[deps.DualNumbers]]
+deps = ["Calculus", "NaNMath", "SpecialFunctions"]
+git-tree-sha1 = "5837a837389fccf076445fce071c8ddaea35a566"
+uuid = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74"
+version = "0.6.8"
+
+[[deps.EarCut_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "e3290f2d49e661fbd94046d7e3726ffcb2d41053"
+uuid = "5ae413db-bbd1-5e63-b57d-d24a61df00f5"
+version = "2.2.4+0"
+
+[[deps.EnumX]]
+git-tree-sha1 = "bdb1942cd4c45e3c678fd11569d5cccd80976237"
+uuid = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
+version = "1.0.4"
+
+[[deps.ErrorfreeArithmetic]]
+git-tree-sha1 = "d6863c556f1142a061532e79f611aa46be201686"
+uuid = "90fa49ef-747e-5e6f-a989-263ba693cf1a"
+version = "0.5.2"
+
+[[deps.ExactPredicates]]
+deps = ["IntervalArithmetic", "Random", "StaticArraysCore", "Test"]
+git-tree-sha1 = "276e83bc8b21589b79303b9985c321024ffdf59c"
+uuid = "429591f6-91af-11e9-00e2-59fbe8cec110"
+version = "2.2.5"
+
 [[deps.ExceptionUnwrapping]]
 deps = ["Test"]
 git-tree-sha1 = "e90caa41f5a86296e014e148ee061bd6c3edec96"
 uuid = "460bff9d-24e4-43bc-9d9f-a8973cb893f4"
 version = "0.1.9"
 
+[[deps.Expat_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "4558ab818dcceaab612d1bb8c19cee87eda2b83c"
+uuid = "2e619515-83b5-522b-bb60-26c02a35a201"
+version = "2.5.0+0"
+
+[[deps.Extents]]
+git-tree-sha1 = "5e1e4c53fa39afe63a7d356e30452249365fba99"
+uuid = "411431e0-e8b7-467b-b5e0-f676ba4f2910"
+version = "0.1.1"
+
+[[deps.FFMPEG]]
+deps = ["FFMPEG_jll"]
+git-tree-sha1 = "b57e3acbe22f8484b4b5ff66a7499717fe1a9cc8"
+uuid = "c87230d0-a227-11e9-1b43-d7ebe4e7570a"
+version = "0.4.1"
+
+[[deps.FFMPEG_jll]]
+deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "PCRE2_jll", "Zlib_jll", "libaom_jll", "libass_jll", "libfdk_aac_jll", "libvorbis_jll", "x264_jll", "x265_jll"]
+git-tree-sha1 = "466d45dc38e15794ec7d5d63ec03d776a9aff36e"
+uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5"
+version = "4.4.4+1"
+
+[[deps.FFTW]]
+deps = ["AbstractFFTs", "FFTW_jll", "LinearAlgebra", "MKL_jll", "Preferences", "Reexport"]
+git-tree-sha1 = "b4fbdd20c889804969571cc589900803edda16b7"
+uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+version = "1.7.1"
+
+[[deps.FFTW_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "c6033cc3892d0ef5bb9cd29b7f2f0331ea5184ea"
+uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a"
+version = "3.3.10+0"
+
+[[deps.FastRounding]]
+deps = ["ErrorfreeArithmetic", "LinearAlgebra"]
+git-tree-sha1 = "6344aa18f654196be82e62816935225b3b9abe44"
+uuid = "fa42c844-2597-5d31-933b-ebd51ab2693f"
+version = "0.3.1"
+
+[[deps.FileIO]]
+deps = ["Pkg", "Requires", "UUIDs"]
+git-tree-sha1 = "299dc33549f68299137e51e6d49a13b5b1da9673"
+uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
+version = "1.16.1"
+
+[[deps.FileWatching]]
+uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"
+
+[[deps.FillArrays]]
+deps = ["LinearAlgebra", "Random"]
+git-tree-sha1 = "a20eaa3ad64254c61eeb5f230d9306e937405434"
+uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
+version = "1.6.1"
+weakdeps = ["SparseArrays", "Statistics"]
+
+    [deps.FillArrays.extensions]
+    FillArraysSparseArraysExt = "SparseArrays"
+    FillArraysStatisticsExt = "Statistics"
+
+[[deps.FiniteDiff]]
+deps = ["ArrayInterface", "LinearAlgebra", "Requires", "Setfield", "SparseArrays"]
+git-tree-sha1 = "c6e4a1fbe73b31a3dea94b1da449503b8830c306"
+uuid = "6a86dc24-6348-571c-b903-95158fe2bd41"
+version = "2.21.1"
+
+    [deps.FiniteDiff.extensions]
+    FiniteDiffBandedMatricesExt = "BandedMatrices"
+    FiniteDiffBlockBandedMatricesExt = "BlockBandedMatrices"
+    FiniteDiffStaticArraysExt = "StaticArrays"
+
+    [deps.FiniteDiff.weakdeps]
+    BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
+    BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
+    StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+
+[[deps.FixedPointNumbers]]
+deps = ["Statistics"]
+git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
+uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
+version = "0.8.4"
+
+[[deps.Fontconfig_jll]]
+deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Pkg", "Zlib_jll"]
+git-tree-sha1 = "21efd19106a55620a188615da6d3d06cd7f6ee03"
+uuid = "a3f928ae-7b40-5064-980b-68af3947d34b"
+version = "2.13.93+0"
+
+[[deps.Formatting]]
+deps = ["Printf"]
+git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
+uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
+version = "0.4.2"
+
+[[deps.ForwardDiff]]
+deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions"]
+git-tree-sha1 = "cf0fe81336da9fb90944683b8c41984b08793dad"
+uuid = "f6369f11-7733-5829-9624-2563aa707210"
+version = "0.10.36"
+weakdeps = ["StaticArrays"]
+
+    [deps.ForwardDiff.extensions]
+    ForwardDiffStaticArraysExt = "StaticArrays"
+
+[[deps.FreeType]]
+deps = ["CEnum", "FreeType2_jll"]
+git-tree-sha1 = "50351f83f95282cf903e968d7c6e8d44a5f83d0b"
+uuid = "b38be410-82b0-50bf-ab77-7b57e271db43"
+version = "4.1.0"
+
+[[deps.FreeType2_jll]]
+deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "Zlib_jll"]
+git-tree-sha1 = "d8db6a5a2fe1381c1ea4ef2cab7c69c2de7f9ea0"
+uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7"
+version = "2.13.1+0"
+
+[[deps.FreeTypeAbstraction]]
+deps = ["ColorVectorSpace", "Colors", "FreeType", "GeometryBasics"]
+git-tree-sha1 = "38a92e40157100e796690421e34a11c107205c86"
+uuid = "663a7486-cb36-511b-a19d-713bb74d65c9"
+version = "0.10.0"
+
+[[deps.FriBidi_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "aa31987c2ba8704e23c6c8ba8a4f769d5d7e4f91"
+uuid = "559328eb-81f9-559d-9380-de523a88c83c"
+version = "1.0.10+0"
+
+[[deps.Future]]
+deps = ["Random"]
+uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"
+
+[[deps.GPUArraysCore]]
+deps = ["Adapt"]
+git-tree-sha1 = "2d6ca471a6c7b536127afccfa7564b5b39227fe0"
+uuid = "46192b85-c4d5-4398-a991-12ede77f4527"
+version = "0.1.5"
+
+[[deps.GeoInterface]]
+deps = ["Extents"]
+git-tree-sha1 = "bb198ff907228523f3dee1070ceee63b9359b6ab"
+uuid = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
+version = "1.3.1"
+
+[[deps.GeometryBasics]]
+deps = ["EarCut_jll", "Extents", "GeoInterface", "IterTools", "LinearAlgebra", "StaticArrays", "StructArrays", "Tables"]
+git-tree-sha1 = "424a5a6ce7c5d97cca7bcc4eac551b97294c54af"
+uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
+version = "0.4.9"
+
+[[deps.Gettext_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "XML2_jll"]
+git-tree-sha1 = "9b02998aba7bf074d14de89f9d37ca24a1a0b046"
+uuid = "78b55507-aeef-58d4-861c-77aaff3498b1"
+version = "0.21.0+0"
+
+[[deps.Glib_jll]]
+deps = ["Artifacts", "Gettext_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE2_jll", "Pkg", "Zlib_jll"]
+git-tree-sha1 = "d3b3624125c1474292d0d8ed0f65554ac37ddb23"
+uuid = "7746bdde-850d-59dc-9ae8-88ece973131d"
+version = "2.74.0+2"
+
 [[deps.Glob]]
 git-tree-sha1 = "97285bbd5230dd766e9ef6749b80fc617126d496"
 uuid = "c27321d9-0574-5035-807b-f59d2c89b15c"
 version = "1.3.1"
 
+[[deps.Graphics]]
+deps = ["Colors", "LinearAlgebra", "NaNMath"]
+git-tree-sha1 = "d61890399bc535850c4bf08e4e0d3a7ad0f21cbd"
+uuid = "a2bd30eb-e257-5431-a919-1863eab51364"
+version = "1.1.2"
+
+[[deps.Graphite2_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "344bf40dcab1073aca04aa0df4fb092f920e4011"
+uuid = "3b182d85-2403-5c21-9c21-1e1f0cc25472"
+version = "1.3.14+0"
+
+[[deps.GridLayoutBase]]
+deps = ["GeometryBasics", "InteractiveUtils", "Observables"]
+git-tree-sha1 = "f57a64794b336d4990d90f80b147474b869b1bc4"
+uuid = "3955a311-db13-416c-9275-1d80ed98e5e9"
+version = "0.9.2"
+
+[[deps.Grisu]]
+git-tree-sha1 = "53bb909d1151e57e2484c3d1b53e19552b887fb2"
+uuid = "42e2da0e-8278-4e71-bc24-59509adca0fe"
+version = "1.0.2"
+
 [[deps.HTTP]]
 deps = ["Base64", "CodecZlib", "ConcurrentUtilities", "Dates", "ExceptionUnwrapping", "Logging", "LoggingExtras", "MbedTLS", "NetworkOptions", "OpenSSL", "Random", "SimpleBufferStream", "Sockets", "URIs", "UUIDs"]
 git-tree-sha1 = "cb56ccdd481c0dd7f975ad2b3b62d9eda088f7e2"
 uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
 version = "1.9.14"
 
+[[deps.HarfBuzz_jll]]
+deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "Graphite2_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg"]
+git-tree-sha1 = "129acf094d168394e80ee1dc4bc06ec835e510a3"
+uuid = "2e76f6c2-a576-52d4-95c1-20adfe4de566"
+version = "2.8.1+1"
+
+[[deps.HypergeometricFunctions]]
+deps = ["DualNumbers", "LinearAlgebra", "OpenLibm_jll", "SpecialFunctions"]
+git-tree-sha1 = "f218fe3736ddf977e0e772bc9a586b2383da2685"
+uuid = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
+version = "0.3.23"
+
 [[deps.IOCapture]]
 deps = ["Logging", "Random"]
 git-tree-sha1 = "d75853a0bdbfb1ac815478bacd89cd27b550ace6"
 uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
 version = "0.2.3"
 
+[[deps.ImageAxes]]
+deps = ["AxisArrays", "ImageBase", "ImageCore", "Reexport", "SimpleTraits"]
+git-tree-sha1 = "2e4520d67b0cef90865b3ef727594d2a58e0e1f8"
+uuid = "2803e5a7-5153-5ecf-9a86-9b4c37f5f5ac"
+version = "0.6.11"
+
+[[deps.ImageBase]]
+deps = ["ImageCore", "Reexport"]
+git-tree-sha1 = "b51bb8cae22c66d0f6357e3bcb6363145ef20835"
+uuid = "c817782e-172a-44cc-b673-b171935fbb9e"
+version = "0.1.5"
+
+[[deps.ImageCore]]
+deps = ["AbstractFFTs", "ColorVectorSpace", "Colors", "FixedPointNumbers", "Graphics", "MappedArrays", "MosaicViews", "OffsetArrays", "PaddedViews", "Reexport"]
+git-tree-sha1 = "acf614720ef026d38400b3817614c45882d75500"
+uuid = "a09fc81d-aa75-5fe9-8630-4744c3626534"
+version = "0.9.4"
+
+[[deps.ImageIO]]
+deps = ["FileIO", "IndirectArrays", "JpegTurbo", "LazyModules", "Netpbm", "OpenEXR", "PNGFiles", "QOI", "Sixel", "TiffImages", "UUIDs"]
+git-tree-sha1 = "bca20b2f5d00c4fbc192c3212da8fa79f4688009"
+uuid = "82e4d734-157c-48bb-816b-45c225c6df19"
+version = "0.6.7"
+
+[[deps.ImageMetadata]]
+deps = ["AxisArrays", "ImageAxes", "ImageBase", "ImageCore"]
+git-tree-sha1 = "355e2b974f2e3212a75dfb60519de21361ad3cb7"
+uuid = "bc367c6b-8a6b-528e-b4bd-a4b897500b49"
+version = "0.9.9"
+
+[[deps.Imath_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "3d09a9f60edf77f8a4d99f9e015e8fbf9989605d"
+uuid = "905a6f67-0a94-5f89-b386-d35d92009cd1"
+version = "3.1.7+0"
+
+[[deps.IndirectArrays]]
+git-tree-sha1 = "012e604e1c7458645cb8b436f8fba789a51b257f"
+uuid = "9b13fd28-a010-5f03-acff-a1bbcff69959"
+version = "1.0.0"
+
+[[deps.Inflate]]
+git-tree-sha1 = "5cd07aab533df5170988219191dfad0519391428"
+uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
+version = "0.1.3"
+
+[[deps.IntegerMathUtils]]
+git-tree-sha1 = "b8ffb903da9f7b8cf695a8bead8e01814aa24b30"
+uuid = "18e54dd8-cb9d-406c-a71d-865a43cbb235"
+version = "0.1.2"
+
+[[deps.IntelOpenMP_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "ad37c091f7d7daf900963171600d7c1c5c3ede32"
+uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0"
+version = "2023.2.0+0"
+
 [[deps.InteractiveUtils]]
 deps = ["Markdown"]
 uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 
+[[deps.Interpolations]]
+deps = ["Adapt", "AxisAlgorithms", "ChainRulesCore", "LinearAlgebra", "OffsetArrays", "Random", "Ratios", "Requires", "SharedArrays", "SparseArrays", "StaticArrays", "WoodburyMatrices"]
+git-tree-sha1 = "721ec2cf720536ad005cb38f50dbba7b02419a15"
+uuid = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+version = "0.14.7"
+
+[[deps.IntervalArithmetic]]
+deps = ["CRlibm", "FastRounding", "LinearAlgebra", "Markdown", "Random", "RecipesBase", "RoundingEmulator", "SetRounding", "StaticArrays"]
+git-tree-sha1 = "5ab7744289be503d76a944784bac3f2df7b809af"
+uuid = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
+version = "0.20.9"
+
+[[deps.IntervalSets]]
+deps = ["Dates", "Random"]
+git-tree-sha1 = "8e59ea773deee525c99a8018409f64f19fb719e6"
+uuid = "8197267c-284f-5f27-9208-e0e47529a953"
+version = "0.7.7"
+weakdeps = ["Statistics"]
+
+    [deps.IntervalSets.extensions]
+    IntervalSetsStatisticsExt = "Statistics"
+
+[[deps.IrrationalConstants]]
+git-tree-sha1 = "630b497eafcc20001bba38a4651b327dcfc491d2"
+uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
+version = "0.2.2"
+
+[[deps.Isoband]]
+deps = ["isoband_jll"]
+git-tree-sha1 = "f9b6d97355599074dc867318950adaa6f9946137"
+uuid = "f1662d9f-8043-43de-a69a-05efc1cc6ff4"
+version = "0.1.1"
+
+[[deps.IterTools]]
+git-tree-sha1 = "4ced6667f9974fc5c5943fa5e2ef1ca43ea9e450"
+uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
+version = "1.8.0"
+
+[[deps.IteratorInterfaceExtensions]]
+git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
+uuid = "82899510-4779-5014-852e-03e436cf321d"
+version = "1.0.0"
+
 [[deps.JLLWrappers]]
 deps = ["Artifacts", "Preferences"]
-git-tree-sha1 = "a7e91ef94114d5bc8952bcaa8d6ff952cf709808"
+git-tree-sha1 = "7e5d6779a1e09a36db2a7b6cff50942a0a7d0fca"
 uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
-version = "1.4.2"
+version = "1.5.0"
 
 [[deps.JSON]]
 deps = ["Dates", "Mmap", "Parsers", "Unicode"]
@@ -93,32 +689,211 @@ git-tree-sha1 = "31e996f0a15c7b280ba9f76636b3ff9e2ae58c9a"
 uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 version = "0.21.4"
 
+[[deps.JpegTurbo]]
+deps = ["CEnum", "FileIO", "ImageCore", "JpegTurbo_jll", "TOML"]
+git-tree-sha1 = "327713faef2a3e5c80f96bf38d1fa26f7a6ae29e"
+uuid = "b835a17e-a41a-41e7-81f0-2f016b05efe0"
+version = "0.1.3"
+
+[[deps.JpegTurbo_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "6f2675ef130a300a112286de91973805fcc5ffbc"
+uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8"
+version = "2.1.91+0"
+
+[[deps.KernelDensity]]
+deps = ["Distributions", "DocStringExtensions", "FFTW", "Interpolations", "StatsBase"]
+git-tree-sha1 = "90442c50e202a5cdf21a7899c66b240fdef14035"
+uuid = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b"
+version = "0.6.7"
+
+[[deps.LAME_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "f6250b16881adf048549549fba48b1161acdac8c"
+uuid = "c1c5ebd0-6772-5130-a774-d5fcae4a789d"
+version = "3.100.1+0"
+
+[[deps.LLVMOpenMP_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "f689897ccbe049adb19a065c495e75f372ecd42b"
+uuid = "1d63c593-3942-5779-bab2-d838dc0a180e"
+version = "15.0.4+0"
+
+[[deps.LZO_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "e5b909bcf985c5e2605737d2ce278ed791b89be6"
+uuid = "dd4b983a-f0e5-5f8d-a1b7-129d4a5fb1ac"
+version = "2.10.1+0"
+
+[[deps.LaTeXStrings]]
+git-tree-sha1 = "f2355693d6778a178ade15952b7ac47a4ff97996"
+uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
+version = "1.3.0"
+
+[[deps.LazyArtifacts]]
+deps = ["Artifacts", "Pkg"]
+uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
+
+[[deps.LazyModules]]
+git-tree-sha1 = "a560dd966b386ac9ae60bdd3a3d3a326062d3c3e"
+uuid = "8cdb02fc-e678-4876-92c5-9defec4f444e"
+version = "0.3.1"
+
+[[deps.LibCURL]]
+deps = ["LibCURL_jll", "MozillaCACerts_jll"]
+uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
+version = "0.6.3"
+
+[[deps.LibCURL_jll]]
+deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
+uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
+version = "7.84.0+0"
+
 [[deps.LibGit2]]
 deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
 uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
 
+[[deps.LibSSH2_jll]]
+deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
+uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
+version = "1.10.2+0"
+
 [[deps.Libdl]]
 uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
 
+[[deps.Libffi_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "0b4a5d71f3e5200a7dff793393e09dfc2d874290"
+uuid = "e9f186c6-92d2-5b65-8a66-fee21dc1b490"
+version = "3.2.2+1"
+
+[[deps.Libgcrypt_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgpg_error_jll", "Pkg"]
+git-tree-sha1 = "64613c82a59c120435c067c2b809fc61cf5166ae"
+uuid = "d4300ac3-e22c-5743-9152-c294e39db1e4"
+version = "1.8.7+0"
+
+[[deps.Libgpg_error_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "c333716e46366857753e273ce6a69ee0945a6db9"
+uuid = "7add5ba3-2f88-524e-9cd5-f83b8a55f7b8"
+version = "1.42.0+0"
+
+[[deps.Libiconv_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "f9557a255370125b405568f9767d6d195822a175"
+uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
+version = "1.17.0+0"
+
+[[deps.Libmount_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "9c30530bf0effd46e15e0fdcf2b8636e78cbbd73"
+uuid = "4b2f31a3-9ecc-558c-b454-b3730dcb73e9"
+version = "2.35.0+0"
+
+[[deps.Libuuid_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "7f3efec06033682db852f8b3bc3c1d2b0a0ab066"
+uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700"
+version = "2.36.0+0"
+
+[[deps.LightXML]]
+deps = ["Libdl", "XML2_jll"]
+git-tree-sha1 = "e129d9391168c677cd4800f5c0abb1ed8cb3794f"
+uuid = "9c8b4983-aa76-5018-a973-4c85ecc9e179"
+version = "0.9.0"
+
+[[deps.LineSearches]]
+deps = ["LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "Printf"]
+git-tree-sha1 = "7bbea35cec17305fc70a0e5b4641477dc0789d9d"
+uuid = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
+version = "7.2.0"
+
+[[deps.LinearAlgebra]]
+deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"]
+uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+
+[[deps.LinearAlgebraX]]
+deps = ["LinearAlgebra", "Mods", "Permutations", "Primes", "SimplePolynomials"]
+git-tree-sha1 = "558a338f1eeabe933f9c2d4052aa7c2c707c3d52"
+uuid = "9b3f67b0-2d00-526e-9884-9e4938f8fb88"
+version = "0.1.12"
+
 [[deps.Literate]]
 deps = ["Base64", "IOCapture", "JSON", "REPL"]
-git-tree-sha1 = "a1a0d4ff9f785a2184baca7a5c89e664f144143d"
+git-tree-sha1 = "f65966c5518ed2b9db90944d7657c312190854f1"
 uuid = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
-version = "2.14.1"
+version = "2.14.2"
+
+[[deps.LogExpFunctions]]
+deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"]
+git-tree-sha1 = "7d6dd4e9212aebaeed356de34ccf262a3cd415aa"
+uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
+version = "0.3.26"
+
+    [deps.LogExpFunctions.extensions]
+    LogExpFunctionsChainRulesCoreExt = "ChainRulesCore"
+    LogExpFunctionsChangesOfVariablesExt = "ChangesOfVariables"
+    LogExpFunctionsInverseFunctionsExt = "InverseFunctions"
+
+    [deps.LogExpFunctions.weakdeps]
+    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+    ChangesOfVariables = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
+    InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112"
 
 [[deps.Logging]]
 uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
 
 [[deps.LoggingExtras]]
 deps = ["Dates", "Logging"]
-git-tree-sha1 = "cedb76b37bc5a6c702ade66be44f831fa23c681e"
+git-tree-sha1 = "0d097476b6c381ab7906460ef1ef1638fbce1d91"
 uuid = "e6f89c97-d47a-5376-807f-9c37f3926c36"
-version = "1.0.0"
+version = "1.0.2"
+
+[[deps.MKL_jll]]
+deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"]
+git-tree-sha1 = "eb006abbd7041c28e0d16260e50a24f8f9104913"
+uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7"
+version = "2023.2.0+0"
+
+[[deps.MacroTools]]
+deps = ["Markdown", "Random"]
+git-tree-sha1 = "9ee1618cbf5240e6d4e0371d6f24065083f60c48"
+uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
+version = "0.5.11"
+
+[[deps.Makie]]
+deps = ["Animations", "Base64", "ColorBrewer", "ColorSchemes", "ColorTypes", "Colors", "Contour", "DelaunayTriangulation", "Distributions", "DocStringExtensions", "Downloads", "FFMPEG", "FileIO", "FixedPointNumbers", "Formatting", "FreeType", "FreeTypeAbstraction", "GeometryBasics", "GridLayoutBase", "ImageIO", "InteractiveUtils", "IntervalSets", "Isoband", "KernelDensity", "LaTeXStrings", "LinearAlgebra", "MacroTools", "MakieCore", "Markdown", "Match", "MathTeXEngine", "Observables", "OffsetArrays", "Packing", "PlotUtils", "PolygonOps", "PrecompileTools", "Printf", "REPL", "Random", "RelocatableFolders", "Setfield", "ShaderAbstractions", "Showoff", "SignedDistanceFields", "SparseArrays", "StableHashTraits", "Statistics", "StatsBase", "StatsFuns", "StructArrays", "TriplotBase", "UnicodeFun"]
+git-tree-sha1 = "e81675589ba7199a82443e87fc52e17eeceac2e8"
+uuid = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+version = "0.19.8"
+
+[[deps.MakieCore]]
+deps = ["Observables"]
+git-tree-sha1 = "f56b09c8b964919373d61750c6d8d4d2c602a2be"
+uuid = "20f20a25-4f0e-4fdf-b5d1-57303727442b"
+version = "0.6.5"
+
+[[deps.MappedArrays]]
+git-tree-sha1 = "2dab0221fe2b0f2cb6754eaa743cc266339f527e"
+uuid = "dbb5928d-eab1-5f90-85c2-b9b0edb7c900"
+version = "0.4.2"
 
 [[deps.Markdown]]
 deps = ["Base64"]
 uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
 
+[[deps.Match]]
+git-tree-sha1 = "1d9bc5c1a6e7ee24effb93f175c9342f9154d97f"
+uuid = "7eb4fadd-790c-5f42-8a69-bfa0b872bfbf"
+version = "1.2.0"
+
+[[deps.MathTeXEngine]]
+deps = ["AbstractTrees", "Automa", "DataStructures", "FreeTypeAbstraction", "GeometryBasics", "LaTeXStrings", "REPL", "RelocatableFolders", "Test", "UnicodeFun"]
+git-tree-sha1 = "8f52dbaa1351ce4cb847d95568cb29e62a307d93"
+uuid = "0a4f8689-d25c-4efe-a92b-7142dfc1aa53"
+version = "0.5.6"
+
 [[deps.MbedTLS]]
 deps = ["Dates", "MbedTLS_jll", "MozillaCACerts_jll", "Random", "Sockets"]
 git-tree-sha1 = "03a9b9718f5682ecb107ac9f7308991db4ce395b"
@@ -130,17 +905,96 @@ deps = ["Artifacts", "Libdl"]
 uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
 version = "2.28.2+0"
 
+[[deps.Missings]]
+deps = ["DataAPI"]
+git-tree-sha1 = "f66bdc5de519e8f8ae43bdc598782d35a25b1272"
+uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
+version = "1.1.0"
+
 [[deps.Mmap]]
 uuid = "a63ad114-7e13-5084-954f-fe012c677804"
 
+[[deps.Mods]]
+git-tree-sha1 = "61be59e4daffff43a8cec04b5e0dc773cbb5db3a"
+uuid = "7475f97c-0381-53b1-977b-4c60186c8d62"
+version = "1.3.3"
+
+[[deps.MosaicViews]]
+deps = ["MappedArrays", "OffsetArrays", "PaddedViews", "StackViews"]
+git-tree-sha1 = "7b86a5d4d70a9f5cdf2dacb3cbe6d251d1a61dbe"
+uuid = "e94cdb99-869f-56ef-bcf0-1ae2bcbe0389"
+version = "0.3.4"
+
 [[deps.MozillaCACerts_jll]]
 uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
 version = "2022.10.11"
 
+[[deps.Multisets]]
+git-tree-sha1 = "8d852646862c96e226367ad10c8af56099b4047e"
+uuid = "3b2b4ff1-bcff-5658-a3ee-dbcf1ce5ac09"
+version = "0.4.4"
+
+[[deps.NLSolversBase]]
+deps = ["DiffResults", "Distributed", "FiniteDiff", "ForwardDiff"]
+git-tree-sha1 = "a0b464d183da839699f4c79e7606d9d186ec172c"
+uuid = "d41bc354-129a-5804-8e4c-c37616107c6c"
+version = "7.8.3"
+
+[[deps.NaNMath]]
+deps = ["OpenLibm_jll"]
+git-tree-sha1 = "0877504529a3e5c3343c6f8b4c0381e57e4387e4"
+uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
+version = "1.0.2"
+
+[[deps.Netpbm]]
+deps = ["FileIO", "ImageCore", "ImageMetadata"]
+git-tree-sha1 = "d92b107dbb887293622df7697a2223f9f8176fcd"
+uuid = "f09324ee-3d7c-5217-9330-fc30815ba969"
+version = "1.1.1"
+
 [[deps.NetworkOptions]]
 uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
 version = "1.2.0"
 
+[[deps.Observables]]
+git-tree-sha1 = "6862738f9796b3edc1c09d0890afce4eca9e7e93"
+uuid = "510215fc-4207-5dde-b226-833fc4488ee2"
+version = "0.5.4"
+
+[[deps.OffsetArrays]]
+deps = ["Adapt"]
+git-tree-sha1 = "2ac17d29c523ce1cd38e27785a7d23024853a4bb"
+uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+version = "1.12.10"
+
+[[deps.Ogg_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "887579a3eb005446d514ab7aeac5d1d027658b8f"
+uuid = "e7412a2a-1a6e-54c0-be00-318e2571c051"
+version = "1.3.5+1"
+
+[[deps.OpenBLAS_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
+uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
+version = "0.3.21+4"
+
+[[deps.OpenEXR]]
+deps = ["Colors", "FileIO", "OpenEXR_jll"]
+git-tree-sha1 = "327f53360fdb54df7ecd01e96ef1983536d1e633"
+uuid = "52e1d378-f018-4a11-a4be-720524705ac7"
+version = "0.3.2"
+
+[[deps.OpenEXR_jll]]
+deps = ["Artifacts", "Imath_jll", "JLLWrappers", "Libdl", "Zlib_jll"]
+git-tree-sha1 = "a4ca623df1ae99d09bc9868b008262d0c0ac1e4f"
+uuid = "18a262bb-aa17-5467-a713-aee519bc75cb"
+version = "3.1.4+0"
+
+[[deps.OpenLibm_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
+version = "0.8.1+0"
+
 [[deps.OpenSSL]]
 deps = ["BitFlags", "Dates", "MozillaCACerts_jll", "OpenSSL_jll", "Sockets"]
 git-tree-sha1 = "51901a49222b09e3743c65b8847687ae5fc78eb2"
@@ -153,17 +1007,139 @@ git-tree-sha1 = "e78db7bd5c26fc5a6911b50a47ee302219157ea8"
 uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95"
 version = "3.0.10+0"
 
+[[deps.OpenSpecFun_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
+uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
+version = "0.5.5+0"
+
+[[deps.Optim]]
+deps = ["Compat", "FillArrays", "ForwardDiff", "LineSearches", "LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "PositiveFactorizations", "Printf", "SparseArrays", "StatsBase"]
+git-tree-sha1 = "963b004d15216f8129f6c0f7d187efa136570be0"
+uuid = "429524aa-4258-5aef-a3af-852621145aeb"
+version = "1.7.7"
+
+[[deps.Opus_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "51a08fb14ec28da2ec7a927c4337e4332c2a4720"
+uuid = "91d4177d-7536-5919-b921-800302f37372"
+version = "1.3.2+0"
+
+[[deps.OrderedCollections]]
+git-tree-sha1 = "2e73fe17cac3c62ad1aebe70d44c963c3cfdc3e3"
+uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+version = "1.6.2"
+
+[[deps.PCRE2_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "efcefdf7-47ab-520b-bdef-62a2eaa19f15"
+version = "10.42.0+0"
+
+[[deps.PDMats]]
+deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
+git-tree-sha1 = "67eae2738d63117a196f497d7db789821bce61d1"
+uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
+version = "0.11.17"
+
+[[deps.PNGFiles]]
+deps = ["Base64", "CEnum", "ImageCore", "IndirectArrays", "OffsetArrays", "libpng_jll"]
+git-tree-sha1 = "9b02b27ac477cad98114584ff964e3052f656a0f"
+uuid = "f57f5aa1-a3ce-4bc8-8ab9-96f992907883"
+version = "0.4.0"
+
+[[deps.Packing]]
+deps = ["GeometryBasics"]
+git-tree-sha1 = "ec3edfe723df33528e085e632414499f26650501"
+uuid = "19eb6ba3-879d-56ad-ad62-d5c202156566"
+version = "0.5.0"
+
+[[deps.PaddedViews]]
+deps = ["OffsetArrays"]
+git-tree-sha1 = "0fac6313486baae819364c52b4f483450a9d793f"
+uuid = "5432bcbf-9aad-5242-b902-cca2824c8663"
+version = "0.5.12"
+
+[[deps.Pango_jll]]
+deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "FriBidi_jll", "Glib_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "84a314e3926ba9ec66ac097e3635e270986b0f10"
+uuid = "36c8627f-9965-5494-a995-c6b170f724f3"
+version = "1.50.9+0"
+
+[[deps.Parameters]]
+deps = ["OrderedCollections", "UnPack"]
+git-tree-sha1 = "34c0e9ad262e5f7fc75b10a9952ca7692cfc5fbe"
+uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a"
+version = "0.12.3"
+
 [[deps.Parsers]]
 deps = ["Dates", "PrecompileTools", "UUIDs"]
 git-tree-sha1 = "716e24b21538abc91f6205fd1d8363f39b442851"
 uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
 version = "2.7.2"
 
+[[deps.Permutations]]
+deps = ["Combinatorics", "LinearAlgebra", "Random"]
+git-tree-sha1 = "6e6cab1c54ae2382bcc48866b91cf949cea703a1"
+uuid = "2ae35dd2-176d-5d53-8349-f30d82d94d4f"
+version = "0.4.16"
+
+[[deps.Pixman_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LLVMOpenMP_jll", "Libdl"]
+git-tree-sha1 = "64779bc4c9784fee475689a1752ef4d5747c5e87"
+uuid = "30392449-352a-5448-841d-b1acce4e97dc"
+version = "0.42.2+0"
+
+[[deps.Pkg]]
+deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
+uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+version = "1.9.2"
+
+[[deps.PkgVersion]]
+deps = ["Pkg"]
+git-tree-sha1 = "f9501cc0430a26bc3d156ae1b5b0c1b47af4d6da"
+uuid = "eebad327-c553-4316-9ea0-9fa01ccd7688"
+version = "0.3.3"
+
+[[deps.PlotUtils]]
+deps = ["ColorSchemes", "Colors", "Dates", "PrecompileTools", "Printf", "Random", "Reexport", "Statistics"]
+git-tree-sha1 = "f92e1315dadf8c46561fb9396e525f7200cdc227"
+uuid = "995b91a9-d308-5afd-9ec6-746e21dbc043"
+version = "1.3.5"
+
+[[deps.PolygonOps]]
+git-tree-sha1 = "77b3d3605fc1cd0b42d95eba87dfcd2bf67d5ff6"
+uuid = "647866c9-e3ac-4575-94e7-e3d426903924"
+version = "0.1.2"
+
+[[deps.Polynomials]]
+deps = ["LinearAlgebra", "RecipesBase", "Setfield"]
+git-tree-sha1 = "af8c8b863adb84abacb9a87822a778a3982901e1"
+uuid = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
+version = "4.0.2"
+
+    [deps.Polynomials.extensions]
+    PolynomialsChainRulesCoreExt = "ChainRulesCore"
+    PolynomialsFFTWExt = "FFTW"
+    PolynomialsMakieCoreExt = "MakieCore"
+    PolynomialsMutableArithmeticsExt = "MutableArithmetics"
+
+    [deps.Polynomials.weakdeps]
+    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+    FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+    MakieCore = "20f20a25-4f0e-4fdf-b5d1-57303727442b"
+    MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
+
+[[deps.PositiveFactorizations]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "17275485f373e6673f7e7f97051f703ed5b15b20"
+uuid = "85a6dd25-e78a-55b7-8502-1745935b8125"
+version = "0.2.4"
+
 [[deps.PrecompileTools]]
 deps = ["Preferences"]
-git-tree-sha1 = "9673d39decc5feece56ef3940e5dafba15ba0f81"
+git-tree-sha1 = "03b4c25b43cb84cee5c90aa9b5ea0a78fd848d2f"
 uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
-version = "1.1.2"
+version = "1.2.0"
 
 [[deps.Preferences]]
 deps = ["TOML"]
@@ -171,10 +1147,34 @@ git-tree-sha1 = "7eb1686b4f04b82f96ed7a4ea5890a4f0c7a09f1"
 uuid = "21216c6a-2e73-6563-6e65-726566657250"
 version = "1.4.0"
 
+[[deps.Primes]]
+deps = ["IntegerMathUtils"]
+git-tree-sha1 = "4c9f306e5d6603ae203c2000dd460d81a5251489"
+uuid = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae"
+version = "0.5.4"
+
 [[deps.Printf]]
 deps = ["Unicode"]
 uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 
+[[deps.ProgressMeter]]
+deps = ["Distributed", "Printf"]
+git-tree-sha1 = "ae36206463b2395804f2787ffe172f44452b538d"
+uuid = "92933f4c-e287-5a05-a399-4b506db050ca"
+version = "1.8.0"
+
+[[deps.QOI]]
+deps = ["ColorTypes", "FileIO", "FixedPointNumbers"]
+git-tree-sha1 = "18e8f4d1426e965c7b532ddd260599e1510d26ce"
+uuid = "4b34888f-f399-49d4-9bb3-47ed5cae4e65"
+version = "1.0.0"
+
+[[deps.QuadGK]]
+deps = ["DataStructures", "LinearAlgebra"]
+git-tree-sha1 = "6ec7ac8412e83d57e313393220879ede1740f9ee"
+uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+version = "2.8.2"
+
 [[deps.REPL]]
 deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
 uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
@@ -183,41 +1183,304 @@ uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
 deps = ["SHA", "Serialization"]
 uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
+[[deps.RangeArrays]]
+git-tree-sha1 = "b9039e93773ddcfc828f12aadf7115b4b4d225f5"
+uuid = "b3c3ace0-ae52-54e7-9d0b-2c1406fd6b9d"
+version = "0.3.2"
+
+[[deps.Ratios]]
+deps = ["Requires"]
+git-tree-sha1 = "1342a47bf3260ee108163042310d26f2be5ec90b"
+uuid = "c84ed2f1-dad5-54f0-aa8e-dbefe2724439"
+version = "0.4.5"
+weakdeps = ["FixedPointNumbers"]
+
+    [deps.Ratios.extensions]
+    RatiosFixedPointNumbersExt = "FixedPointNumbers"
+
+[[deps.RecipesBase]]
+deps = ["PrecompileTools"]
+git-tree-sha1 = "5c3d09cc4f31f5fc6af001c250bf1278733100ff"
+uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+version = "1.3.4"
+
 [[deps.Reexport]]
 git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
 uuid = "189a3867-3050-52da-a836-e630ba90ab69"
 version = "1.2.2"
 
+[[deps.RelocatableFolders]]
+deps = ["SHA", "Scratch"]
+git-tree-sha1 = "90bc7a7c96410424509e4263e277e43250c05691"
+uuid = "05181044-ff0b-4ac5-8273-598c1e38db00"
+version = "1.0.0"
+
+[[deps.Requires]]
+deps = ["UUIDs"]
+git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
+uuid = "ae029012-a4dd-5104-9daa-d747884805df"
+version = "1.3.0"
+
+[[deps.RingLists]]
+deps = ["Random"]
+git-tree-sha1 = "9712ebc42e91850f35272b48eb840e60c0270ec0"
+uuid = "286e9d63-9694-5540-9e3c-4e6708fa07b2"
+version = "0.2.7"
+
+[[deps.Rmath]]
+deps = ["Random", "Rmath_jll"]
+git-tree-sha1 = "f65dcb5fa46aee0cf9ed6274ccbd597adc49aa7b"
+uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
+version = "0.7.1"
+
+[[deps.Rmath_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "6ed52fdd3382cf21947b15e8870ac0ddbff736da"
+uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
+version = "0.4.0+0"
+
+[[deps.RoundingEmulator]]
+git-tree-sha1 = "40b9edad2e5287e05bd413a38f61a8ff55b9557b"
+uuid = "5eaf0fd0-dfba-4ccb-bf02-d820a40db705"
+version = "0.2.1"
+
 [[deps.SHA]]
 uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 version = "0.7.0"
 
+[[deps.Scratch]]
+deps = ["Dates"]
+git-tree-sha1 = "30449ee12237627992a99d5e30ae63e4d78cd24a"
+uuid = "6c6a2e73-6563-6170-7368-637461726353"
+version = "1.2.0"
+
 [[deps.Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 
+[[deps.SetRounding]]
+git-tree-sha1 = "d7a25e439d07a17b7cdf97eecee504c50fedf5f6"
+uuid = "3cc68bcd-71a2-5612-b932-767ffbe40ab0"
+version = "0.2.1"
+
+[[deps.Setfield]]
+deps = ["ConstructionBase", "Future", "MacroTools", "StaticArraysCore"]
+git-tree-sha1 = "e2cc6d8c88613c05e1defb55170bf5ff211fbeac"
+uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46"
+version = "1.1.1"
+
+[[deps.ShaderAbstractions]]
+deps = ["ColorTypes", "FixedPointNumbers", "GeometryBasics", "LinearAlgebra", "Observables", "StaticArrays", "StructArrays", "Tables"]
+git-tree-sha1 = "0d15c3e7b2003f4451714f08ffec2b77badc2dc4"
+uuid = "65257c39-d410-5151-9873-9b3e5be5013e"
+version = "0.3.0"
+
+[[deps.SharedArrays]]
+deps = ["Distributed", "Mmap", "Random", "Serialization"]
+uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
+
+[[deps.Showoff]]
+deps = ["Dates", "Grisu"]
+git-tree-sha1 = "91eddf657aca81df9ae6ceb20b959ae5653ad1de"
+uuid = "992d4aef-0814-514b-bc4d-f2e9a6c4116f"
+version = "1.0.3"
+
+[[deps.SignedDistanceFields]]
+deps = ["Random", "Statistics", "Test"]
+git-tree-sha1 = "d263a08ec505853a5ff1c1ebde2070419e3f28e9"
+uuid = "73760f76-fbc4-59ce-8f25-708e95d2df96"
+version = "0.4.0"
+
 [[deps.SimpleBufferStream]]
 git-tree-sha1 = "874e8867b33a00e784c8a7e4b60afe9e037b74e1"
 uuid = "777ac1f9-54b0-4bf8-805c-2214025038e7"
 version = "1.1.0"
 
+[[deps.SimpleGraphs]]
+deps = ["AbstractLattices", "Combinatorics", "DataStructures", "IterTools", "LightXML", "LinearAlgebra", "LinearAlgebraX", "Optim", "Primes", "Random", "RingLists", "SimplePartitions", "SimplePolynomials", "SimpleRandom", "SparseArrays", "Statistics"]
+git-tree-sha1 = "b608903049d11cc557c45e03b3a53e9260579c19"
+uuid = "55797a34-41de-5266-9ec1-32ac4eb504d3"
+version = "0.8.4"
+
+[[deps.SimplePartitions]]
+deps = ["AbstractLattices", "DataStructures", "Permutations"]
+git-tree-sha1 = "dcc02923a53f316ab97da8ef3136e80b4543dbf1"
+uuid = "ec83eff0-a5b5-5643-ae32-5cbf6eedec9d"
+version = "0.3.0"
+
+[[deps.SimplePolynomials]]
+deps = ["Mods", "Multisets", "Polynomials", "Primes"]
+git-tree-sha1 = "d537c31cf9995236166e3e9afc424a5a1c59ff9d"
+uuid = "cc47b68c-3164-5771-a705-2bc0097375a0"
+version = "0.2.14"
+
+[[deps.SimpleRandom]]
+deps = ["Distributions", "LinearAlgebra", "Random"]
+git-tree-sha1 = "3a6fb395e37afab81aeea85bae48a4db5cd7244a"
+uuid = "a6525b86-64cd-54fa-8f65-62fc48bdc0e8"
+version = "0.3.1"
+
+[[deps.SimpleTraits]]
+deps = ["InteractiveUtils", "MacroTools"]
+git-tree-sha1 = "5d7e3f4e11935503d3ecaf7186eac40602e7d231"
+uuid = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
+version = "0.9.4"
+
+[[deps.Sixel]]
+deps = ["Dates", "FileIO", "ImageCore", "IndirectArrays", "OffsetArrays", "REPL", "libsixel_jll"]
+git-tree-sha1 = "2da10356e31327c7096832eb9cd86307a50b1eb6"
+uuid = "45858cf5-a6b0-47a3-bbea-62219f50df47"
+version = "0.1.3"
+
 [[deps.Sockets]]
 uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
 
+[[deps.SortingAlgorithms]]
+deps = ["DataStructures"]
+git-tree-sha1 = "c60ec5c62180f27efea3ba2908480f8055e17cee"
+uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
+version = "1.1.1"
+
+[[deps.SparseArrays]]
+deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"]
+uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+
+[[deps.SpecialFunctions]]
+deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
+git-tree-sha1 = "e2cfc4012a19088254b3950b85c3c1d8882d864d"
+uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
+version = "2.3.1"
+weakdeps = ["ChainRulesCore"]
+
+    [deps.SpecialFunctions.extensions]
+    SpecialFunctionsChainRulesCoreExt = "ChainRulesCore"
+
+[[deps.StableHashTraits]]
+deps = ["CRC32c", "Compat", "Dates", "SHA", "Tables", "TupleTools", "UUIDs"]
+git-tree-sha1 = "0b8b801b8f03a329a4e86b44c5e8a7d7f4fe10a3"
+uuid = "c5dd0088-6c3f-4803-b00e-f31a60c170fa"
+version = "0.3.1"
+
+[[deps.StackViews]]
+deps = ["OffsetArrays"]
+git-tree-sha1 = "46e589465204cd0c08b4bd97385e4fa79a0c770c"
+uuid = "cae243ae-269e-4f55-b966-ac2d0dc13c15"
+version = "0.1.1"
+
+[[deps.StaticArrays]]
+deps = ["LinearAlgebra", "Random", "StaticArraysCore"]
+git-tree-sha1 = "9cabadf6e7cd2349b6cf49f1915ad2028d65e881"
+uuid = "90137ffa-7385-5640-81b9-e52037218182"
+version = "1.6.2"
+weakdeps = ["Statistics"]
+
+    [deps.StaticArrays.extensions]
+    StaticArraysStatisticsExt = "Statistics"
+
+[[deps.StaticArraysCore]]
+git-tree-sha1 = "36b3d696ce6366023a0ea192b4cd442268995a0d"
+uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
+version = "1.4.2"
+
+[[deps.Statistics]]
+deps = ["LinearAlgebra", "SparseArrays"]
+uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+version = "1.9.0"
+
+[[deps.StatsAPI]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "45a7769a04a3cf80da1c1c7c60caf932e6f4c9f7"
+uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
+version = "1.6.0"
+
+[[deps.StatsBase]]
+deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
+git-tree-sha1 = "75ebe04c5bed70b91614d684259b661c9e6274a4"
+uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+version = "0.34.0"
+
+[[deps.StatsFuns]]
+deps = ["HypergeometricFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
+git-tree-sha1 = "f625d686d5a88bcd2b15cd81f18f98186fdc0c9a"
+uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
+version = "1.3.0"
+
+    [deps.StatsFuns.extensions]
+    StatsFunsChainRulesCoreExt = "ChainRulesCore"
+    StatsFunsInverseFunctionsExt = "InverseFunctions"
+
+    [deps.StatsFuns.weakdeps]
+    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+    InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112"
+
+[[deps.StructArrays]]
+deps = ["Adapt", "DataAPI", "GPUArraysCore", "StaticArraysCore", "Tables"]
+git-tree-sha1 = "521a0e828e98bb69042fec1809c1b5a680eb7389"
+uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
+version = "0.6.15"
+
+[[deps.SuiteSparse]]
+deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
+uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
+
+[[deps.SuiteSparse_jll]]
+deps = ["Artifacts", "Libdl", "Pkg", "libblastrampoline_jll"]
+uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c"
+version = "5.10.1+6"
+
 [[deps.TOML]]
 deps = ["Dates"]
 uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 version = "1.0.3"
 
+[[deps.TableTraits]]
+deps = ["IteratorInterfaceExtensions"]
+git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
+uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
+version = "1.0.1"
+
+[[deps.Tables]]
+deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "OrderedCollections", "TableTraits", "Test"]
+git-tree-sha1 = "1544b926975372da01227b382066ab70e574a3ec"
+uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
+version = "1.10.1"
+
+[[deps.Tar]]
+deps = ["ArgTools", "SHA"]
+uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
+version = "1.10.0"
+
+[[deps.TensorCore]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "1feb45f88d133a655e001435632f019a9a1bcdb6"
+uuid = "62fd8b95-f654-4bbd-a8a5-9c27f68ccd50"
+version = "0.1.1"
+
 [[deps.Test]]
 deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
 uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
+[[deps.TiffImages]]
+deps = ["ColorTypes", "DataStructures", "DocStringExtensions", "FileIO", "FixedPointNumbers", "IndirectArrays", "Inflate", "Mmap", "OffsetArrays", "PkgVersion", "ProgressMeter", "UUIDs"]
+git-tree-sha1 = "8621f5c499a8aa4aa970b1ae381aae0ef1576966"
+uuid = "731e570b-9d59-4bfa-96dc-6df516fadf69"
+version = "0.6.4"
+
 [[deps.TranscodingStreams]]
 deps = ["Random", "Test"]
 git-tree-sha1 = "9a6ae7ed916312b41236fcef7e0af564ef934769"
 uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
 version = "0.9.13"
 
+[[deps.TriplotBase]]
+git-tree-sha1 = "4d4ed7f294cda19382ff7de4c137d24d16adc89b"
+uuid = "981d1d27-644d-49a2-9326-4793e63143c3"
+version = "0.1.0"
+
+[[deps.TupleTools]]
+git-tree-sha1 = "3c712976c47707ff893cf6ba4354aa14db1d8938"
+uuid = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6"
+version = "1.3.0"
+
 [[deps.URIs]]
 git-tree-sha1 = "b7a5e99f24892b6824a954199a45e9ffcc1c70f0"
 uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
@@ -227,15 +1490,156 @@ version = "1.5.0"
 deps = ["Random", "SHA"]
 uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
 
+[[deps.UnPack]]
+git-tree-sha1 = "387c1f73762231e86e0c9c5443ce3b4a0a9a0c2b"
+uuid = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
+version = "1.0.2"
+
 [[deps.Unicode]]
 uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
 
+[[deps.UnicodeFun]]
+deps = ["REPL"]
+git-tree-sha1 = "53915e50200959667e78a92a418594b428dffddf"
+uuid = "1cfade01-22cf-5700-b092-accc4b62d6e1"
+version = "0.4.1"
+
+[[deps.WoodburyMatrices]]
+deps = ["LinearAlgebra", "SparseArrays"]
+git-tree-sha1 = "de67fa59e33ad156a590055375a30b23c40299d3"
+uuid = "efce3f68-66dc-5838-9240-27a6d6f5f9b6"
+version = "0.5.5"
+
+[[deps.XML2_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Zlib_jll"]
+git-tree-sha1 = "04a51d15436a572301b5abbb9d099713327e9fc4"
+uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
+version = "2.10.4+0"
+
+[[deps.XSLT_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgcrypt_jll", "Libgpg_error_jll", "Libiconv_jll", "Pkg", "XML2_jll", "Zlib_jll"]
+git-tree-sha1 = "91844873c4085240b95e795f692c4cec4d805f8a"
+uuid = "aed1982a-8fda-507f-9586-7b0439959a61"
+version = "1.1.34+0"
+
+[[deps.Xorg_libX11_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libxcb_jll", "Xorg_xtrans_jll"]
+git-tree-sha1 = "afead5aba5aa507ad5a3bf01f58f82c8d1403495"
+uuid = "4f6342f7-b3d2-589e-9d20-edeb45f2b2bc"
+version = "1.8.6+0"
+
+[[deps.Xorg_libXau_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "6035850dcc70518ca32f012e46015b9beeda49d8"
+uuid = "0c0b7dd1-d40b-584c-a123-a41640f87eec"
+version = "1.0.11+0"
+
+[[deps.Xorg_libXdmcp_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "34d526d318358a859d7de23da945578e8e8727b7"
+uuid = "a3789734-cfe1-5b06-b2d0-1dd0d9d62d05"
+version = "1.1.4+0"
+
+[[deps.Xorg_libXext_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
+git-tree-sha1 = "b7c0aa8c376b31e4852b360222848637f481f8c3"
+uuid = "1082639a-0dae-5f34-9b06-72781eeb8cb3"
+version = "1.3.4+4"
+
+[[deps.Xorg_libXrender_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
+git-tree-sha1 = "19560f30fd49f4d4efbe7002a1037f8c43d43b96"
+uuid = "ea2f1a96-1ddc-540d-b46f-429655e07cfa"
+version = "0.9.10+4"
+
+[[deps.Xorg_libpthread_stubs_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "8fdda4c692503d44d04a0603d9ac0982054635f9"
+uuid = "14d82f49-176c-5ed1-bb49-ad3f5cbd8c74"
+version = "0.1.1+0"
+
+[[deps.Xorg_libxcb_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "XSLT_jll", "Xorg_libXau_jll", "Xorg_libXdmcp_jll", "Xorg_libpthread_stubs_jll"]
+git-tree-sha1 = "b4bfde5d5b652e22b9c790ad00af08b6d042b97d"
+uuid = "c7cfdc94-dc32-55de-ac96-5a1b8d977c5b"
+version = "1.15.0+0"
+
+[[deps.Xorg_xtrans_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "e92a1a012a10506618f10b7047e478403a046c77"
+uuid = "c5fb5394-a638-5e4d-96e5-b29de1b5cf10"
+version = "1.5.0+0"
+
 [[deps.Zlib_jll]]
 deps = ["Libdl"]
 uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
 version = "1.2.13+0"
 
+[[deps.isoband_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "51b5eeb3f98367157a7a12a1fb0aa5328946c03c"
+uuid = "9a68df92-36a6-505f-a73e-abb412b6bfb4"
+version = "0.2.3+0"
+
+[[deps.libaom_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "3a2ea60308f0996d26f1e5354e10c24e9ef905d4"
+uuid = "a4ae2306-e953-59d6-aa16-d00cac43593b"
+version = "3.4.0+0"
+
+[[deps.libass_jll]]
+deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
+git-tree-sha1 = "5982a94fcba20f02f42ace44b9894ee2b140fe47"
+uuid = "0ac62f75-1d6f-5e53-bd7c-93b484bb37c0"
+version = "0.15.1+0"
+
+[[deps.libblastrampoline_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
+version = "5.8.0+0"
+
+[[deps.libfdk_aac_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "daacc84a041563f965be61859a36e17c4e4fcd55"
+uuid = "f638f0a6-7fb0-5443-88ba-1cc74229b280"
+version = "2.0.2+0"
+
+[[deps.libpng_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
+git-tree-sha1 = "94d180a6d2b5e55e447e2d27a29ed04fe79eb30c"
+uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f"
+version = "1.6.38+0"
+
+[[deps.libsixel_jll]]
+deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Pkg", "libpng_jll"]
+git-tree-sha1 = "d4f63314c8aa1e48cd22aa0c17ed76cd1ae48c3c"
+uuid = "075b6546-f08a-558a-be8f-8157d0f608a5"
+version = "1.10.3+0"
+
+[[deps.libvorbis_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Ogg_jll", "Pkg"]
+git-tree-sha1 = "b910cb81ef3fe6e78bf6acee440bda86fd6ae00c"
+uuid = "f27f6e37-5d2b-51aa-960f-b287f2bc3b7a"
+version = "1.3.7+1"
+
+[[deps.nghttp2_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
+version = "1.48.0+0"
+
 [[deps.p7zip_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
 version = "17.4.0+0"
+
+[[deps.x264_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "4fea590b89e6ec504593146bf8b988b2c00922b2"
+uuid = "1270edf5-f2f9-52d2-97e9-ab00b5d0237a"
+version = "2021.5.5+0"
+
+[[deps.x265_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "ee567a171cce03570d77ad3a43e90218e38937a9"
+uuid = "dfaa095f-4041-5dcd-9319-2fabd8486b76"
+version = "3.5.0+0"
diff --git a/docs/Project.toml b/docs/Project.toml
index 377e2b2..e11571d 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,4 +1,5 @@
 [deps]
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 DataDeps = "124859b0-ceae-595e-8997-d05f6a7a8dfe"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
diff --git a/docs/make.jl b/docs/make.jl
index 86d64a9..a75c5ff 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -14,15 +14,21 @@ ENV["DATADEPS_ALWAYS_ACCEPT"] = "true"
 #####
 
 const EXAMPLES_DIR = joinpath(@__DIR__, "..", "examples")
-const OUTPUT_DIR   = joinpath(@__DIR__, "src/literated")
+const OUTPUT_DIR   = joinpath(@__DIR__, "src/generated")
 
-to_be_literated = []
+example_scripts = [
+    "freezing_bucket.jl",
+]
 
-for file in to_be_literated
-    filepath = joinpath(EXAMPLES_DIR, file)
+for filename in example_scripts
+    filepath = joinpath(EXAMPLES_DIR, filename)
     Literate.markdown(filepath, OUTPUT_DIR; flavor = Literate.DocumenterFlavor())
 end
 
+example_pages = [
+    "Freezing bucket" => "generated/freezing_bucket.md",
+]
+
 #####
 ##### Build and deploy docs
 #####
@@ -35,6 +41,7 @@ format = Documenter.HTML(
 
 pages = [
     "Home" => "index.md",
+    "Examples" => example_pages,
 
     "Library" => [ 
         "Contents"       => "library/outline.md",
diff --git a/examples/diffusive_ice_column_model.jl b/examples/diffusive_ice_column_model.jl
index bc6ff05..18e399b 100644
--- a/examples/diffusive_ice_column_model.jl
+++ b/examples/diffusive_ice_column_model.jl
@@ -1,4 +1,4 @@
-using ClimaSeaIce.EulerianThermodynamicSeaIceModels: EulerianThermodynamicSeaIceModel, MolecularDiffusivity
+using ClimaSeaIce.EnthalpyMethodSeaIceModels: EnthalpyMethodSeaIceModel, MolecularDiffusivity
 using Oceananigans
 using Oceananigans.Units
 using Oceananigans.Operators: Δzᶜᶜᶠ
@@ -6,7 +6,7 @@ using Oceananigans.Grids: znode
 using GLMakie
 
 #####
-##### Set up a EulerianThermodynamicSeaIceModel
+##### Set up a EnthalpyMethodSeaIceModel
 #####
 
 # Build a grid with 10 cm resolution
@@ -49,7 +49,7 @@ top_T_bc = ValueBoundaryCondition(air_ice_temperature)
 bottom_T_bc = ValueBoundaryCondition(ice_ocean_temperature)
 T_bcs = FieldBoundaryConditions(top=top_T_bc, bottom=bottom_T_bc)
 
-model = EulerianThermodynamicSeaIceModel(; grid, closure, boundary_conditions=(; T=T_bcs))
+model = EnthalpyMethodSeaIceModel(; grid, closure, boundary_conditions=(; T=T_bcs))
 
 # Initialize and run
 set!(model, T=initial_ice_ocean_temperature)
diff --git a/examples/freezing_bucket.jl b/examples/freezing_bucket.jl
new file mode 100644
index 0000000..24e7bd0
--- /dev/null
+++ b/examples/freezing_bucket.jl
@@ -0,0 +1,69 @@
+using Oceananigans
+using Oceananigans.Units
+using ClimaSeaIce
+
+# # Setting up a model of a bucket of water freezing from the top down
+#
+# We'll model a slab of ice in a bucket with one grid point, which
+# requires only a zero-dimensional grid.
+
+grid = RectilinearGrid(size=(), topology=(Flat, Flat, Flat))
+
+## Set the temperature at the top of the bucket to ``Tᵤ = -10ᵒC``.
+model = SlabSeaIceModel(grid; top_thermal_boundary_condition=PrescribedTemperature(-10))
+
+## We'll freeze the bucket for 10 straight days
+simulation = Simulation(model, Δt=10minute, stop_time=10days)
+
+## Initialize the ice thickness to 1 cm
+set!(model, h=0.01)
+
+# # Collecting data and running the simulation
+#
+# Before running the simulation, we set up a `Callback` to create
+# a timeseries of the ice thickness saved at every time step.
+
+## Container to hold the data
+timeseries = []
+
+## Callback function to collect the data from the `sim`ulation
+function accumulate_timeseries(sim)
+    h = sim.model.ice_thickness
+    push!(timeseries, (time(sim), first(h)))
+end
+
+## Add the callback to `simulation`
+simulation.callbacks[:save] = Callback(accumulate_timeseries)
+
+# Now we're ready to hit the Big Red Button (it should run pretty quick):
+run!(simulation)
+
+# # Visualizing the result
+#
+# It'd be a shame to run such a "cool" simulation without looking at the
+# results. We'll visualize it with Makie.
+
+using CairoMakie
+
+# `timeseries` is a `Vector` of `Tuple`. So we have to do a bit of processing
+# to build `Vector`s of time `t` and thickness `h`. It's not much work though:
+t = map(ts -> ts[1], timeseries)
+h = map(ts -> ts[2], timeseries)
+
+# Just for fun, we also compute the velocity of the ice-water interface:
+dhdt = @. (h[2:end] - h[1:end-1]) / simulation.Δt
+
+# All that's left, really, is to put those `lines!` in an `Axis`:
+set_theme!(Theme(fontsize=24, linewidth=4))
+
+fig = Figure(resolution=(1200, 600))
+
+axh = Axis(fig[1, 1], xlabel="Time (days)", ylabel="Ice thickness (cm)")
+axd = Axis(fig[1, 2], xlabel="Ice thickness (cm)", ylabel="Freezing rate (μm s⁻¹)")
+
+lines!(axh, t ./ day, 1e2 .* h)
+lines!(axd, 1e2 .* h[1:end-1], 1e6 .* dhdt)
+
+current_figure() # hide
+fig
+
diff --git a/examples/ice_ocean_column_model.jl b/examples/ice_ocean_column_model.jl
deleted file mode 100644
index 303ba09..0000000
--- a/examples/ice_ocean_column_model.jl
+++ /dev/null
@@ -1,112 +0,0 @@
-using ClimaSeaIce.EulerianThermodynamicSeaIceModels: EulerianThermodynamicSeaIceModel, MolecularDiffusivity, temperature_flux
-using Oceananigans
-using Oceananigans.Units
-using Oceananigans.Operators: Δzᶜᶜᶠ
-using Oceananigans.TurbulenceClosures: CATKEVerticalDiffusivity
-using Oceananigans.Grids: znode
-using GLMakie
-
-#####
-##### Set up an Ocean column model
-#####
-
-ocean_grid = RectilinearGrid(size=20, z=(-100, 0), topology=(Flat, Flat, Bounded))
-equation_of_state = LinearEquationOfState(thermal_expansion=2e-4, haline_contraction=8e-5)
-
-# Set up boundary conditions ready for coupling (on the CPU)
-top_T_flux = zeros(1, 1)
-top_S_flux = zeros(1, 1)
-T_bcs = FieldBoundaryConditions(top=FluxBoundaryCondition(top_T_flux))
-S_bcs = FieldBoundaryConditions(top=FluxBoundaryCondition(top_S_flux))
-
-ocean_model = HydrostaticFreeSurfaceModel(grid = ocean_grid,
-                                          closure = CATKEVerticalDiffusivity(),
-                                          tracers = (:T, :S, :e),
-                                          buoyancy = SeawaterBuoyancy(; equation_of_state),
-                                          boundary_conditions = (; T=T_bcs, S=S_bcs))
-
-
-# Ocean initial condition
-α = equation_of_state.thermal_expansion
-g = ocean_model.buoyancy.model.gravitational_acceleration
-N² = 1e-6
-dTdz = - α * g * N²
-T₀ = 0.1 # Temperature at surface
-Tᵢ(x, y, z) = T₀ + dTdz * z
-Sᵢ(x, y, z) = 35
-set!(ocean_model, T=Tᵢ, S=Sᵢ)
-
-#####
-##### Set up a EulerianThermodynamicSeaIceModel
-#####
-
-# Build a grid with 10 cm resolution
-ice_grid = RectilinearGrid(size=20, z=(-1, 0), topology=(Flat, Flat, Bounded))
-
-# Set up a simple problem and build the ice model
-ice_closure = MolecularDiffusivity(ice_grid, κ_ice=1e-5, κ_water=1e-6)
-
-#####
-##### Create temperature boundary conditions
-#####
-
-# Atmopshere-ice surface temperature parameters
-initial_air_ice_temperature = -5
-air_ice_temperature_fluctuation_amplitude = 5
-air_ice_tendency = -0.5 / day # ᵒC s⁻¹
-
-parameters = (Tᵢ = initial_air_ice_temperature,
-              Tₐ = air_ice_temperature_fluctuation_amplitude,
-              dTdt = air_ice_tendency)
-
-@inline air_ice_temperature(x, y, t, p) = p.dTdt * t + p.Tₐ * sin(2π*t / day) + p.Tᵢ
-top_T_bc = ValueBoundaryCondition(air_ice_temperature; parameters)
-T_bcs = FieldBoundaryConditions(top=top_T_bc)
-
-ice_model = EulerianThermodynamicSeaIceModel(grid = ice_grid,
-                                     closure = ice_closure,
-                                     boundary_conditions=(; T=T_bcs))
-
-
-κ = 1e-5
-Δz = Δzᶜᶜᶠ(1, 1, 1, ice_grid)
-Δt = 0.1 * Δz^2 / κ
-
-ice_simulation = Simulation(ice_model; Δt, stop_iteration=10)
-
-@inline is_ocean(i, j, k, grid, ϕ) = @inbounds ϕ[i, j, k] == 1
-
-function fix_liquid_temperature!(ϕ_ice_model, T_ice_model, grid, T_ocean_model) 
-    Nz = size(grid, 3)
-
-    i = j = 1
-    @inbounds for k = 1:Nz
-        T_ice_model[i, j, k] = ifelse(is_ocean(i, j, k, grid, ϕ_ice_model),
-                                      T_ocean_model[i, j],
-                                      T_ice_model[i, j, k])
-    end
-
-    return nothing
-end
-
-function set_liquid_state!(ice_simulation)
-    ice_model = ice_simulation.model
-
-    # Extract ocean temperature at top grid point
-    ocean_Nz = size(ocean_model.grid, 3)
-    ocean_temperature = interior(ocean_model.tracers.T, :, :, ocean_Nz)
-
-    # Eventually this will launch a kernel on the GPU
-    ϕ = ice_model.state.ϕ
-    T = ice_model.state.T
-    fix_liquid_temperature!(ϕ, T, ice_model.grid, ocean_temperature)
-
-
-    return nothing
-end
-
-ice_simulation.callbacks[:ocean_state] = Callback(set_liquid_state!)
-
-#=
-ocean_simulation = Simulation(ocean_model; Δt)
-=#
diff --git a/examples/idealized_boundary_conditions.jl b/examples/idealized_boundary_conditions.jl
deleted file mode 100644
index 37bdd31..0000000
--- a/examples/idealized_boundary_conditions.jl
+++ /dev/null
@@ -1,35 +0,0 @@
-using Oceananigans.Units
-using GLMakie
-
-# This file creates idealized ice-interface temperatures
-# for driving an ice column model.
-
-# Information about the diurnal temperature fluctuations and cooling trends at the atmosphere boundary
-top_T_initial = -5
-top_T_amplitude = 5
-top_T_slope = -0.5 / day # ᵒC s⁻¹
-
-# Information about ocean cooling
-bottom_T_initial = 0
-bottom_T_slope = -0.5 / day # ᵒC s⁻¹
-
-# Information about time steps
-dt = 1000
-tf = 864000
-time = 0:dt:tf
-
-# Calculate BCs
-top_T_BC = @. top_T_slope * time + top_T_amplitude * sin(2π*time/day) + top_T_initial
-bottom_T_BC = @. bottom_T_slope * time + bottom_T_initial
-
-# Plot boundary conditions
-
-set_theme!(Theme(fontsize=24, linewidth=3))
-fig = Figure()
-ax = Axis(fig[1, 1], title="Boundary Conditions", xlabel="Time (s)", ylabel="Temperature (ᵒC)")
-lines!(ax, time, top_T_BC, label="Air-ice surface temperature")
-lines!(ax, time, bottom_T_BC, label="Ice-ocean temperature")
-axislegend(ax)
-     
-display(fig)
-
diff --git a/examples/melting_in_march.jl b/examples/melting_in_march.jl
new file mode 100644
index 0000000..1bc7444
--- /dev/null
+++ b/examples/melting_in_march.jl
@@ -0,0 +1,73 @@
+using Oceananigans
+using Oceananigans.Units
+using ClimaSeaIce
+using GLMakie
+
+# Generate a 0D grid for a single column slab model 
+grid = RectilinearGrid(size=(), topology=(Flat, Flat, Flat))
+
+# Build a model of an ice slab that has internal conductive fluxes
+# and that emits radiation from its top surface.
+solar_insolation = -500 # W m⁻²
+outgoing_radiation = RadiativeEmission()
+
+# Define a FluxFunction representing a sensible heat flux
+parameters = (
+    transfer_coefficient     = 1e-3,  # Unitless
+    atmosphere_density       = 1.225, # kg m⁻³
+    atmosphere_heat_capacity = 1004,  # 
+    atmosphere_temperature   = -5,    # ᵒC
+    atmosphere_wind_speed    = 5      # m s⁻¹
+)
+
+@inline function sensible_heat_flux(i, j, grid, Tᵤ, clock, fields, parameters)
+    Cₛ = parameters.transfer_coefficient
+    ρₐ = parameters.atmosphere_density
+    cₐ = parameters.atmosphere_heat_capacity
+    Tₐ = parameters.atmosphere_temperature
+    uₐ = parameters.atmosphere_wind_speed
+
+    # Flux is positive (cooling by fluxing heat up away from upper surface)
+    # when Tₐ < Tᵤ:
+    return Cₛ * ρₐ * cₐ * uₐ * (Tᵤ - Tₐ)
+end
+
+aerodynamic_flux = FluxFunction(sensible_heat_flux; parameters)
+
+top_thermal_flux = (outgoing_radiation, solar_insolation, aerodynamic_flux)
+model = SlabSeaIceModel(grid; top_thermal_flux)
+set!(model, h=1)
+
+simulation = Simulation(model, Δt=10minute, stop_time=1day)
+
+# Accumulate data
+timeseries = []
+
+function accumulate_timeseries(sim)
+    T = model.top_temperature
+    h = model.ice_thickness
+    push!(timeseries, (time(sim), first(h), first(T)))
+end
+
+simulation.callbacks[:save] = Callback(accumulate_timeseries)
+
+run!(simulation)
+
+# Extract and visualize data
+
+t = [datum[1] for datum in timeseries]
+h = [datum[2] for datum in timeseries]
+T = [datum[3] for datum in timeseries]
+
+set_theme!(Theme(fontsize=24, linewidth=4))
+
+fig = Figure(resolution=(1000, 800))
+
+axT = Axis(fig[1, 1], xlabel="Time (days)", ylabel="Top temperature (ᵒC)")
+axh = Axis(fig[2, 1], xlabel="Time (days)", ylabel="Ice thickness (m)")
+
+lines!(axT, t / day, T)
+lines!(axh, t / day, h)
+
+display(fig)
+
diff --git a/examples/perpetual_night.jl b/examples/perpetual_night.jl
new file mode 100644
index 0000000..ef7a850
--- /dev/null
+++ b/examples/perpetual_night.jl
@@ -0,0 +1,47 @@
+using Oceananigans
+using Oceananigans.Units
+using ClimaSeaIce
+using ClimaSeaIce.ThermalBoundaryConditions: RadiativeEmission
+using GLMakie
+
+# Generate a zero-dimensional grid for a single column slab model 
+grid = RectilinearGrid(size=(), topology=(Flat, Flat, Flat))
+
+# Build a model of an ice slab that has internal conductive fluxes
+# and that emits radiation from its top surface.
+model = SlabSeaIceModel(grid; top_thermal_flux=RadiativeEmission())
+set!(model, h=0.01)
+
+simulation = Simulation(model, Δt=1hour, stop_time=40days)
+
+# Accumulate data
+timeseries = []
+
+function accumulate_timeseries(sim)
+    T = model.top_temperature
+    h = model.ice_thickness
+    push!(timeseries, (time(sim), first(h), first(T)))
+end
+
+simulation.callbacks[:save] = Callback(accumulate_timeseries)
+
+run!(simulation)
+
+# Extract and visualize data
+
+t = [datum[1] for datum in timeseries]
+h = [datum[2] for datum in timeseries]
+T = [datum[3] for datum in timeseries]
+
+set_theme!(Theme(fontsize=24, linewidth=4))
+
+fig = Figure(resolution=(1000, 800))
+
+axT = Axis(fig[1, 1], xlabel="Time (days)", ylabel="Top temperature (ᵒC)")
+axh = Axis(fig[2, 1], xlabel="Time (days)", ylabel="Ice thickness (m)")
+
+lines!(axT, t / day, T)
+lines!(axh, t / day, h)
+
+display(fig)
+
diff --git a/src/ClimaSeaIce.jl b/src/ClimaSeaIce.jl
index 252fb18..3ba44d9 100644
--- a/src/ClimaSeaIce.jl
+++ b/src/ClimaSeaIce.jl
@@ -1,8 +1,140 @@
 module ClimaSeaIce
 
-include("EulerianThermodynamicSeaIceModels.jl")
-include("PrognosticThicknessThermodynamicsModels.jl")
+import Oceananigans.TimeSteppers: time_step!
 
-using .EulerianThermodynamicSeaIceModels: EulerianThermodynamicSeaIceModel
+export
+    MeltingConstrainedFluxBalance,
+    PrescribedTemperature,
+    RadiativeEmission,
+    ConductiveFlux,
+    FluxFunction,
+    SlabSeaIceModel
+
+#####
+##### A bit of thermodynamics to start the day
+#####
+
+struct LinearLiquidus{FT}
+    freshwater_melting_temperature :: FT
+    slope :: FT
+end
+
+"""
+    LinearLiquidus(FT=Float64,
+                   slope = 0.054, # psu / ᵒC
+                   freshwater_melting_temperature = 0) # ᵒC
+
+Return a linear model for the dependence of the melting temperature of
+saltwater on salinity,
+
+```math
+Tₘ(S) = T₀ - m S ,
+```
+
+where ``Tₘ(S)`` is the melting temperature as a function of salinity ``S``,
+``T₀`` is the melting temperature of freshwater, and ``m`` is the ratio
+between the melting temperature and salinity (in other words the linear model
+should be thought of as defining ``m`` and could be written ``m ≡ (T₀ - Tₘ) / S``.
+The signs are arranged so that ``m > 0`` for saltwater).
+
+The defaults assume that salinity is given in practical salinity units `psu` and
+temperature is in degrees Celsius.
+
+Note: the function `melting_temperature(liquidus, salinity)` returns the
+melting temperature given `salinity`.
+"""
+function LinearLiquidus(FT::DataType=Float64;
+                        slope = 0.054, # psu / ᵒC
+                        freshwater_melting_temperature = 0) # ᵒC
+
+    return LinearLiquidus(convert(FT, freshwater_melting_temperature),
+                          convert(FT, slope))
+end
+
+@inline function melting_temperature(liquidus::LinearLiquidus, salinity)
+    return liquidus.freshwater_melting_temperature - liquidus.slope * salinity
+end
+
+struct PhaseTransitions{FT, L}
+    ice_density :: FT
+    ice_heat_capacity :: FT
+    liquid_density :: FT
+    liquid_heat_capacity :: FT
+    reference_latent_heat :: FT
+    reference_temperature :: FT
+    liquidus :: L
+end
+
+"""
+    PhaseTransitions(FT=Float64,
+                     ice_density           = 917,   # kg m⁻³
+                     ice_heat_capacity     = 2000,  # J / (kg ᵒC)
+                     liquid_density        = 999.8, # kg m⁻³
+                     liquid_heat_capacity  = 4186,  # J / (kg ᵒC)
+                     reference_latent_heat = 334e3  # J kg⁻³
+                     liquidus = LinearLiquidus(FT)) # default assumes psu, ᵒC
+
+Return a representation of transitions between the solid and liquid phases
+of salty water: in other words, the freezing and melting of sea ice.
+
+The latent heat of fusion ``ℒ(T)`` (more simply just "latent heat") is
+a function of temperature ``T`` via
+
+```math
+ρᵢ ℒ(T) = ρᵢ ℒ₀ + (ρℓ cℓ - ρᵢ cᵢ) * (T - T₀)    
+```
+
+where ``ρᵢ`` is the `ice_density`, ``ρℓ`` is the liquid density,
+``cᵢ`` is the heat capacity of ice, and ``cℓ`` is the heat capacity of
+liquid, and ``T₀`` is a reference temperature, all of which are assumed constant.
+
+The default `liquidus` assumes that salinity has practical salinity units (psu)
+and that temperature is degrees Celsius.
+"""
+@inline function PhaseTransitions(FT=Float64,
+                                  ice_density           = 917,    # kg m⁻³
+                                  ice_heat_capacity     = 2000,   # J / (kg ᵒC)
+                                  liquid_density        = 999.8,  # kg m⁻³
+                                  liquid_heat_capacity  = 4186,   # J / (kg ᵒC)
+                                  reference_latent_heat = 334e3,  # J kg⁻³
+                                  reference_temperature = 0,      # ᵒC
+                                  liquidus = LinearLiquidus(FT))
+
+    return PhaseTransitions(convert(FT, ice_density),
+                            convert(FT, ice_heat_capacity),
+                            convert(FT, liquid_density),
+                            convert(FT, liquid_heat_capacity),
+                            convert(FT, reference_latent_heat),
+                            convert(FT, reference_temperature),
+                            liquidus)
+end
+
+@inline function latent_heat(thermo::PhaseTransitions, T)
+    T₀ = thermo.reference_temperature    
+    ℒ₀ = thermo.reference_latent_heat
+    ρᵢ = thermo.ice_density
+    ρℓ = thermo.liquid_density
+    cᵢ = thermo.ice_heat_capacity
+    cℓ = thermo.liquid_heat_capacity
+
+    return ρᵢ * ℒ₀ + (ρℓ * cℓ - ρᵢ * cᵢ) * (T - T₀)
+end
+
+struct ForwardEulerTimestepper end
+
+include("ThermalBoundaryConditions/ThermalBoundaryConditions.jl")
+
+using .ThermalBoundaryConditions:
+    IceWaterThermalEquilibrium,
+    MeltingConstrainedFluxBalance,
+    RadiativeEmission,
+    FluxFunction,
+    PrescribedTemperature
+
+include("EnthalpyMethodSeaIceModels.jl")
+include("SlabSeaIceModels.jl")
+
+using .EnthalpyMethodSeaIceModels: EnthalpyMethodSeaIceModel
+using .SlabSeaIceModels: SlabSeaIceModel, ConductiveFlux
 
 end # module
diff --git a/src/EulerianThermodynamicSeaIceModels.jl b/src/EnthalpyMethodSeaIceModels.jl
similarity index 79%
rename from src/EulerianThermodynamicSeaIceModels.jl
rename to src/EnthalpyMethodSeaIceModels.jl
index 6bebf50..78f6a61 100644
--- a/src/EulerianThermodynamicSeaIceModels.jl
+++ b/src/EnthalpyMethodSeaIceModels.jl
@@ -1,4 +1,4 @@
-module EulerianThermodynamicSeaIceModels
+module EnthalpyMethodSeaIceModels
 
 using Oceananigans.BoundaryConditions:
     fill_halo_regions!,
@@ -21,14 +21,14 @@ import Oceananigans.Fields: set!
 import Oceananigans.TimeSteppers: time_step!, update_state!
 import Oceananigans.Simulations: reset!
 
-mutable struct EulerianThermodynamicSeaIceModel{Grid,
-                                                Tim,
-                                                Clk,
-                                                Clo,
-                                                State,
-                                                Cp,
-                                                Fu,
-                                                Tend} <: AbstractModel{Nothing}
+mutable struct EnthalpyMethodSeaIceModel{Grid,
+                                         Tim,
+                                         Clk,
+                                         Clo,
+                                         State,
+                                         Cp,
+                                         Fu,
+                                         Tend} <: AbstractModel{Nothing}
     grid :: Grid
     timestepper :: Tim # unused placeholder for now
     clock :: Clk
@@ -40,27 +40,27 @@ mutable struct EulerianThermodynamicSeaIceModel{Grid,
     tendencies :: Tend
 end
 
-const ETSIM = EulerianThermodynamicSeaIceModel
+const ETSIM = EnthalpyMethodSeaIceModel
 
 function Base.show(io::IO, model::ETSIM)
     clock = model.clock
-    print(io, "EulerianThermodynamicSeaIceModel(t=", prettytime(clock.time), ", iteration=", clock.iteration, ")", '\n')
+    print(io, "EnthalpyMethodSeaIceModel(t=", prettytime(clock.time), ", iteration=", clock.iteration, ")", '\n')
     print(io, "    grid: ", summary(model.grid))
 end
 
 const reference_density = 999.8 # kg m⁻³
 
 """
-    EulerianThermodynamicSeaIceModel(; grid, kw...)
+    EnthalpyMethodSeaIceModel(; grid, kw...)
 
 Return a thermodynamic model for ice sandwiched between an atmosphere and ocean on an Eulerian grid.
 """
-function EulerianThermodynamicSeaIceModel(; grid,
-                                          closure = default_closure(grid),
-                                          ice_heat_capacity = 2090.0 / reference_density,
-                                          water_heat_capacity = 3991.0 / reference_density,
-                                          fusion_enthalpy = 3.3e5 / reference_density,
-                                          boundary_conditions = NamedTuple())
+function EnthalpyMethodSeaIceModel(; grid,
+                                   closure = default_closure(grid),
+                                   ice_heat_capacity = 2090.0 / reference_density,
+                                   water_heat_capacity = 3991.0 / reference_density,
+                                   fusion_enthalpy = 3.3e5 / reference_density,
+                                   boundary_conditions = NamedTuple())
 
     # Prognostic fields: temperature, enthalpy, porosity
     field_names = (:T, :H, :ϕ)
@@ -74,7 +74,7 @@ function EulerianThermodynamicSeaIceModel(; grid,
     tendencies = (; H=CenterField(grid))
     clock = Clock{eltype(grid)}(0, 0, 1)
 
-    return EulerianThermodynamicSeaIceModel(grid,
+    return EnthalpyMethodSeaIceModel(grid,
                                     nothing,
                                     clock,
                                     closure,
@@ -99,7 +99,7 @@ function set!(model::ETSIM; T=nothing, H=nothing)
 
     if setting_enthalpy
         set!(model.state.H, H)
-        update_temperature!(model)
+        compute_temperature!(model)
     end
 
     return nothing
@@ -116,7 +116,7 @@ prognostic_fields(model::ETSIM) = (; model.state.H)
 ##### Time-stepping
 #####
 
-function update_porosity!(model)
+function compute_porosity!(model)
     T = model.state.T
     ϕ = model.state.ϕ
     grid = model.grid
@@ -137,7 +137,7 @@ end
     end
 end
 
-function update_temperature!(model)
+function compute_temperature!(model)
     H = model.state.H
     T = model.state.T
     c = model.ice_heat_capacity
@@ -171,13 +171,13 @@ function update_state!(model::ETSIM)
     grid = model.grid
     arch = grid.architecture
     args = (model.clock, fields(model))
-    update_temperature!(model)
-    update_porosity!(model)
-    update_diffusivity!(model.closure, model)
+    compute_temperature!(model)
+    compute_porosity!(model)
+    compute_diffusivity!(model.closure, model)
     return nothing
 end
 
-function time_step!(model, Δt; callbacks=nothing)
+function time_step!(model::ETSIM, Δt; callbacks=nothing)
     grid = model.grid
     arch = grid.architecture
     Ψ = model.state
@@ -219,13 +219,12 @@ struct MolecularDiffusivity{C}
     κ :: C
 end
 
-
 function MolecularDiffusivity(grid; κ_ice=1e-5, κ_water=1e-6)
     κ = CenterField(grid)
     return MolecularDiffusivity(κ_ice, κ_water, κ)
 end
 
-function update_diffusivity!(closure::MolecularDiffusivity, model)
+function compute_diffusivity!(closure::MolecularDiffusivity, model)
     κ = closure.κ
     κ_ice = closure.κ_ice
     κ_water = closure.κ_water
diff --git a/src/PrognosticThicknessThermodynamicsModels.jl b/src/PrognosticThicknessThermodynamicsModels.jl
deleted file mode 100644
index 467aba2..0000000
--- a/src/PrognosticThicknessThermodynamicsModels.jl
+++ /dev/null
@@ -1,15 +0,0 @@
-module PrognosticThicknessThermodynamicsModels
-
-abstract type AbstractPrognosticThicknessThermodynamicsModel end
-
-mutable struct ZeroLayerThermodynamicSeaIceModel{G, H, A, L, KS, KI} <: AbstractPrognosticThicknessThermodynamicsModel
-    horizontal_grid :: G
-    thicknesses :: H
-    previous_thicknesses :: H
-    albedo :: A
-    latent_heat_of_fusion :: L
-    snow_diffusivity :: KS
-    ice_diffusivity :: KI
-end
-
-end # module
diff --git a/src/SlabSeaIceModels.jl b/src/SlabSeaIceModels.jl
new file mode 100644
index 0000000..bea3ee0
--- /dev/null
+++ b/src/SlabSeaIceModels.jl
@@ -0,0 +1,256 @@
+module SlabSeaIceModels
+
+using ClimaSeaIce:
+    PhaseTransitions,
+    latent_heat,
+    ForwardEulerTimestepper,
+    melting_temperature
+
+using ClimaSeaIce.ThermalBoundaryConditions:
+    MeltingConstrainedFluxBalance,
+    IceWaterThermalEquilibrium,
+    PrescribedTemperature,
+    FluxFunction,
+    SurfaceTemperatureDependent,
+    bottom_temperature,
+    top_temperature,
+    top_flux_imbalance,
+    bottom_flux_imbalance
+
+# using RootSolvers: find_zero
+
+using Oceananigans.Fields: Field, Center, ZeroField, ConstantField
+using Oceananigans.TimeSteppers: Clock, tick!
+
+# Simulations interface
+import Oceananigans: fields, prognostic_fields
+import Oceananigans.Fields: field, set!
+import Oceananigans.TimeSteppers: time_step!, update_state!
+import Oceananigans.Simulations: reset!, initialize!
+import Oceananigans.OutputWriters: default_included_properties
+
+# TODO: move to Oceananigans
+field(loc, a::Number, grid) = ConstantField(a)
+
+struct SlabSeaIceModel{GR, CL, TS, IT, ST, IS, STF, TBC, CF, P}
+    grid :: GR
+    clock :: CL
+    timestepper :: TS
+    # State
+    ice_thickness :: IT
+    top_temperature :: ST
+    ice_salinity :: IS
+    # Boundary conditions
+    external_thermal_fluxes :: STF
+    thermal_boundary_conditions :: TBC
+    # Internal flux
+    internal_thermal_flux :: CF
+    # Melting and freezing stuff
+    phase_transitions :: P
+end
+
+const SSIM = SlabSeaIceModel
+
+initialize!(::SSIM) = nothing
+default_included_properties(::SSIM) = tuple(:grid)
+
+fields(model::SSIM) = (h = model.ice_thickness,
+                       Tᵤ = model.top_temperature,
+                       Sᵢ = model.ice_salinity)
+
+# TODO: make this correct
+prognostic_fields(model::SSIM) = fields(model)
+
+struct ConductiveFlux{K}
+    conductivity :: K
+end
+
+ConductiveFlux(FT::DataType=Float64; conductivity) = ConductiveFlux(convert(FT, conductivity))
+
+@inline function slab_internal_thermal_flux(conductive_flux::ConductiveFlux,
+                                            top_temperature,
+                                            bottom_temperature,
+                                            ice_thickness)
+
+    k = conductive_flux.conductivity
+    Tu = top_temperature
+    Tb = bottom_temperature
+    h = ice_thickness
+
+    return - k * (Tu - Tb) / h
+end
+
+@inline function slab_internal_thermal_flux(i, j, grid,
+                                            top_temperature::Number,
+                                            clock, fields, parameters)
+    flux = parameters.flux
+    bottom_bc = parameters.bottom_thermal_boundary_condition
+    liquidus = parameters.liquidus
+    Tu = top_temperature
+    Tb = bottom_temperature(i, j, grid, bottom_bc, liquidus)
+    h = @inbounds fields.h[i, j, 1]
+    return slab_internal_thermal_flux(flux, Tu, Tb, h)
+end
+
+"""
+    SlabSeaIceModel(grid; kw...)
+
+Pretty simple model for sea ice.
+"""
+function SlabSeaIceModel(grid;
+                         clock                             = Clock{eltype(grid)}(0, 0, 1),
+                         ice_thickness                     = Field{Center, Center, Nothing}(grid),
+                         ice_salinity                      = 0, # psu
+                         top_temperature                   = nothing,
+                         top_thermal_flux                  = 0,
+                         bottom_thermal_flux               = 0,
+                         top_thermal_boundary_condition    = MeltingConstrainedFluxBalance(),
+                         bottom_thermal_boundary_condition = IceWaterThermalEquilibrium(),
+                         internal_thermal_flux             = ConductiveFlux(eltype(grid), conductivity=2),
+                         phase_transitions                 = PhaseTransitions(eltype(grid)))
+
+    # Only one time-stepper is supported currently
+    timestepper = ForwardEulerTimestepper()
+    FT = eltype(grid)
+
+    # TODO: pass `clock` into `field`, so functions can be time-dependent?
+
+    # Wrap ice_salinity in a field
+    ice_salinity = field((Center, Center, Nothing), ice_salinity, grid)
+
+    # Construct default top temperature if one is not provided
+    if isnothing(top_temperature)
+        # Check top boundary condition
+        if top_thermal_boundary_condition isa PrescribedTemperature  
+            top_temperature = top_thermal_boundary_condition.temperature 
+        else # build the default
+            top_temperature = Field{Center, Center, Nothing}(grid)
+        end
+    end
+
+    # Convert to `field` (does nothing if it's already a Field)
+    top_temperature = field((Center, Center, Nothing), top_temperature, grid)
+
+    # Construct an internal thermal flux function that captures the liquidus and
+    # bottom boundary condition.
+    parameters = (flux = internal_thermal_flux,
+                  liquidus = phase_transitions.liquidus,
+                  bottom_thermal_boundary_condition = bottom_thermal_boundary_condition)
+
+    internal_thermal_flux_function = FluxFunction(slab_internal_thermal_flux;
+                                                  parameters,
+                                                  top_temperature_dependent=true)
+
+    # Package the external fluxes and boundary conditions
+    external_thermal_fluxes = (top = top_thermal_flux,    
+                               bottom = bottom_thermal_flux) 
+
+    thermal_boundary_conditions = (top = top_thermal_boundary_condition,
+                                   bottom = bottom_thermal_boundary_condition)
+
+    return SlabSeaIceModel(grid,
+                           clock,
+                           timestepper,
+                           ice_thickness,
+                           top_temperature,
+                           ice_salinity,
+                           external_thermal_fluxes,
+                           thermal_boundary_conditions,
+                           internal_thermal_flux_function,
+                           phase_transitions)
+end
+
+function set!(model::SSIM; h=nothing)
+    !isnothing(h) && set!(model.ice_thickness, h)
+    return nothing
+end
+
+function time_step!(model::SSIM, Δt; callbacks=nothing)
+    grid = model.grid
+    Nx, Ny, Nz = size(grid)
+
+    phase_transitions = model.phase_transitions
+    liquidus = phase_transitions.liquidus
+
+    h = model.ice_thickness
+    Tu = model.top_temperature
+    model_fields = fields(model)
+    ice_salinity = model.ice_salinity
+    clock = model.clock
+    bottom_thermal_bc = model.thermal_boundary_conditions.bottom
+    top_thermal_bc = model.thermal_boundary_conditions.top
+
+    Qi = model.internal_thermal_flux
+    Qb = model.external_thermal_fluxes.bottom
+    Qs = model.external_thermal_fluxes.top
+
+    for i = 1:Nx, j=1:Ny
+
+        @inbounds begin
+            Tuⁿ = Tu[i, j, 1]
+            S = ice_salinity[i, j, 1]
+        end
+
+        # 1. Update thickness with ForwardEuler step
+        # 1a. Calculate residual fluxes across the top top and bottom
+        δQs = top_flux_imbalance(i, j, grid, top_thermal_bc, Tuⁿ, Qi, Qs, clock, model_fields)
+        δQb =  bottom_flux_imbalance(i, j, grid, bottom_thermal_bc,  Tuⁿ, Qi, Qb, clock, model_fields)
+
+        # 1b. Impose an implicit flux balance if Tu ≤ Tm.
+        Tₘ = melting_temperature(liquidus, S) 
+        δQs = ifelse(Tuⁿ <= Tₘ, zero(grid), δQs)
+
+        # 1c. Calculate the latent heat of fusion at the top and bottom
+        Tbⁿ = bottom_temperature(i, j, grid, bottom_thermal_bc, liquidus)
+        ℒ_Tu = latent_heat(phase_transitions, Tuⁿ)
+        ℒ_Tb = latent_heat(phase_transitions, Tbⁿ)
+
+        # 1d. Increment the thickness
+        @inbounds begin
+            h⁺ = h[i, j, 1] + Δt * (δQb / ℒ_Tb + δQs / ℒ_Tu)
+            h⁺ = max(zero(grid), h⁺)
+            h[i, j, 1] = h⁺
+        end
+
+        if !isa(top_thermal_bc, PrescribedTemperature)
+            # 2. Update top temperature
+            Tu⁺ = top_temperature(i, j, grid, top_thermal_bc, Tuⁿ, Qi, Qs, clock, model_fields)
+
+            # 3. Set top temperature (could skip this for `PrescribedTemperature`)
+            @inbounds Tu[i, j, 1] = Tu⁺ 
+        end
+    end
+
+    tick!(model.clock, Δt)
+
+    return nothing
+end
+
+function update_state!(model::SSIM)
+
+    top_thermal_bc = model.thermal_boundary_conditions.top
+
+    if !isa(top_thermal_bc, PrescribedTemperature)
+        grid = model.grid
+        Nx, Ny, Nz = size(grid)
+
+        Tu = model.top_temperature
+        model_fields = fields(model)
+        clock = model.clock
+        top_thermal_bc = model.thermal_boundary_conditions.top
+        Qi = model.internal_thermal_flux
+        Qs = model.external_thermal_fluxes.top
+
+        # Update top temperature
+        for i = 1:Nx, j=1:Ny
+            Tu⁻ = @inbounds Tu[i, j, 1]
+            Tuⁿ = top_temperature(i, j, grid, top_thermal_bc, Tu⁻, Qi, Qs, clock, model_fields)
+            @inbounds Tu[i, j, 1] = Tuⁿ
+        end
+
+    end
+
+    return nothing
+end
+
+end # module
diff --git a/src/ThermalBoundaryConditions/ThermalBoundaryConditions.jl b/src/ThermalBoundaryConditions/ThermalBoundaryConditions.jl
new file mode 100644
index 0000000..b64e095
--- /dev/null
+++ b/src/ThermalBoundaryConditions/ThermalBoundaryConditions.jl
@@ -0,0 +1,17 @@
+module ThermalBoundaryConditions
+
+"""
+    PrescribedTemperature()
+
+Thermal boundary condition indicating that temperature is prescribed on the boundary.
+"""
+struct PrescribedTemperature{T}
+    temperature :: T
+end
+
+include("bottom_thermal_boundary_conditions.jl")
+include("top_thermal_boundary_conditions.jl")
+include("boundary_fluxes.jl")
+
+end
+
diff --git a/src/ThermalBoundaryConditions/bottom_thermal_boundary_conditions.jl b/src/ThermalBoundaryConditions/bottom_thermal_boundary_conditions.jl
new file mode 100644
index 0000000..20a69e7
--- /dev/null
+++ b/src/ThermalBoundaryConditions/bottom_thermal_boundary_conditions.jl
@@ -0,0 +1,54 @@
+using ClimaSeaIce: melting_temperature
+
+#####
+##### Bottom thermal boundary conditions
+#####
+
+struct IceWaterThermalEquilibrium{S}
+    salinity :: S
+end
+
+"""
+    IceWaterThermalEquilibrium(; external_fluxes::Tuple=tuple(), salinity=0)
+
+Represents an ice-water interface in thermal equilibrium, such that
+the bottom temperature ``T_b`` is equal to the melting temperature,
+
+```math
+T_b = Tₘ(S)
+```
+
+where ``S`` is the `salinity` at the ice-water boundary.
+
+Both freezing and melting may occur at an ice-water boundary.
+The rate of change of the ice-water interface, ``d h_b / dt``,
+is related to the sum of `fluxes` into the ice-water interface,
+
+"""
+IceWaterThermalEquilibrium(; salinity=0) = IceWaterThermalEquilibrium(salinity)
+
+@inline bottom_temperature(i, j, grid, bc::PrescribedTemperature, args...) = @inbounds bc.temperature[i, j]
+@inline bottom_temperature(i, j, grid, bc::PrescribedTemperature{<:Number}, args...) = bc.temperature
+
+@inline function bottom_temperature(i, j, grid, bc::IceWaterThermalEquilibrium, liquidus)
+    Sₒ = @inbounds bc.salinity[i, j, 1]
+    return melting_temperature(liquidus, Sₒ)
+end
+
+@inline function bottom_flux_imbalance(i, j, grid, bottom_thermal_bc, top_temperature,
+                                       internal_fluxes, external_fluxes, clock, model_fields)
+
+    #          
+    #   ice        ↑   Qi ≡ internal_fluxes. Example: Qi = - k ∂z T
+    #            |⎴⎴⎴|
+    # ----------------------- ↕ hᵇ → ℒ ∂t hᵇ = δQ _given_ T = Tₘ
+    #            |⎵⎵⎵|
+    #   water      ↑   Qx ≡ external_fluxes
+    #        
+
+    Qi = getflux(internal_fluxes, i, j, grid, top_temperature, clock, model_fields)
+    Qx = getflux(external_fluxes, i, j, grid, top_temperature, clock, model_fields)
+
+    return Qi - Qx
+end
+
diff --git a/src/ThermalBoundaryConditions/boundary_fluxes.jl b/src/ThermalBoundaryConditions/boundary_fluxes.jl
new file mode 100644
index 0000000..cb40150
--- /dev/null
+++ b/src/ThermalBoundaryConditions/boundary_fluxes.jl
@@ -0,0 +1,96 @@
+#####
+##### Fluxes
+#####
+
+# Flux extractor function
+# getflux(flux, i, j, grid, T_top, clock, model_fields)
+
+@inline getflux(flux::Number, i, j, grid, args...) = flux
+@inline getflux(flux::AbstractArray{<:Any, 2}, i, j, grid, args...) = @inbounds flux[i, j]
+
+# Make sure to check the size on flux construction
+@inline getflux(flux::AbstractArray{<:Any, 3}, i, j, grid, args...) = @inbounds flux[i, j, 1]
+
+# Tuples
+@inline getflux(fluxes::Tuple{}, i, j, grid, args...) = zero(grid)
+@inline getflux(fluxes::Tuple{<:Any}, args...) = getflux(fluxes[1], args...)
+@inline getflux(fluxes::Tuple{<:Any, <:Any}, args...) = getflux(fluxes[1], args...) + getflux(fluxes[2], args...)
+@inline getflux(fluxes::Tuple, args...) = getflux(fluxes[1], args...) + _getflux(fluxes[2:end], args...)
+
+@inline _getflux(fluxes::Tuple{}, i, j, grid, args...) = zero(grid)
+@inline _getflux(fluxes::Tuple{<:Any}, grid, args...) = getflux(fluxes[1], args...)
+@inline _getflux(fluxes::Tuple, args...) = getflux(fluxes[1], args...) + getflux(fluxes[2:end], args...)
+
+struct SurfaceTemperatureDependent end
+
+struct FluxFunction{P, T, F}
+    func :: F
+    parameters :: P
+    function FluxFunction{SurfaceTemperatureDependence}(func, parameters=nothing) where SurfaceTemperatureDependence
+        T = SurfaceTemperatureDependence
+        P = typeof(parameters)
+        F = typeof(func)
+        new{P, T, F}(func, parameters)
+    end
+end
+
+"""
+    FluxFunction(func; parameters=nothing, top_temperature_dependent=false)
+
+Return `FluxFunction` representing a flux across an air-ice, air-snow, or ice-water interface.
+The flux is computed by `func` with the signature
+
+```julia
+flux = func(i, j, grid, clock, top_temperature, model_fields)
+```
+
+if `isnothing(parameters)`, or
+
+```julia
+flux = func(i, j, grid, clock, top_temperature, model_fields, parameters)
+```
+
+if `!isnothing(parameters)`. If `func` is `top_temperature_dependent`, then it will be recomputed
+during a diagnostic solve for the top temperature.
+"""
+function FluxFunction(func; parameters=nothing, top_temperature_dependent=false)
+    T = top_temperature_dependent ? SurfaceTemperatureDependent() : Nothing
+    return FluxFunction{T}(func, parameters)
+end
+
+@inline getflux(flux::FluxFunction{<:Nothing}, args...) = flux.func(args...)
+@inline getflux(flux::FluxFunction, args...) = flux.func(args..., flux.parameters)
+
+#####
+##### Implementations
+#####
+
+
+struct RadiativeEmission{FT}
+    emissivity :: FT
+    stefan_boltzmann_constant :: FT
+    reference_temperature :: FT
+end
+
+"""
+    RadiativeEmission(FT=Float64; kw...)
+
+Returns a flux representing radiative emission from a surface.
+"""
+function RadiativeEmission(FT=Float64;
+                           emissivity = 1,
+                           stefan_boltzmann_constant = 5.67e-8,
+                           reference_temperature = 273.15)
+
+    return RadiativeEmission(convert(FT, emissivity),
+                             convert(FT, stefan_boltzmann_constant), 
+                             convert(FT, reference_temperature))
+end
+
+@inline function getflux(emission::RadiativeEmission, i, j, grid, T, clock, fields)
+    ϵ = emission.emissivity
+    σ = emission.stefan_boltzmann_constant
+    Tᵣ = emission.reference_temperature
+    return ϵ * σ * (T + Tᵣ)^4
+end
+
diff --git a/src/ThermalBoundaryConditions/top_thermal_boundary_conditions.jl b/src/ThermalBoundaryConditions/top_thermal_boundary_conditions.jl
new file mode 100644
index 0000000..73f12c8
--- /dev/null
+++ b/src/ThermalBoundaryConditions/top_thermal_boundary_conditions.jl
@@ -0,0 +1,103 @@
+#####
+##### Surface thermal boundary conditions
+#####
+
+struct MeltingConstrainedFluxBalance{STS}
+    top_temperature_solver :: STS
+end
+
+struct LinearizedSurfaceTemperatureSolver end
+struct NonlinearSurfaceTemperatureSolver end
+
+"""
+    MeltingConstrainedFluxBalance(top_temperature_solver=NonlinearSurfaceTemperatureSolver())
+
+Return a boundary condition that determines the top (or "upper surface") temperature ``Tᵤ``
+to equilibrate the top heat fluxes,
+
+```math
+Qₓ₁(Tᵤ) + Qₓ₂(Tᵤ) + ⋯ - Qᵢ = Σᴺ Qₓₙ(Tᵤ) - Qᵢ = δQ ,
+```
+
+where ``Qᵢ`` is the intrinsic flux into the top from within the ice
+(typically, a conductive flux), ``Qₓₙ`` represent external fluxes into the
+air above the ice, ``δQ`` is the residual flux, and ``Tᵤ`` is the (upper)
+top temperature. ``Tᵤ`` is evaluated under the constraint that
+
+```math
+Tᵤ ≤ Tₘ(S),
+```
+
+where ``Tₘ(S)`` is the melting temperature, which is a function of the
+ice salinity at the top, ``S``. When ``Tᵤ < Tₘ(S)``, the top is
+frozen and ``δQ = 0``. When the constraint operates, such that ``Tᵤ = Tₘ(S)``,
+the top is melting and the residual flux is non-zero.
+
+```math
+δQ ≡ Σᴺ Qₙ(Tₘ) - Qᵢ(Tₘ).
+```
+
+The residual flux is consumed by the cost of transforming ice into liquid water,
+and is related to the rate of change of ice thickness, ``h``, by
+
+```math
+d/dt hₛ = δQ / ℒ(Tᵤ)
+```
+
+where ``ℒ(Tᵤ)`` is the latent heat, equal to the different between
+the higher internal energy of liquid water and the lower internal energy of solid ice,
+at the temperature ``Tᵤ``.
+
+"""
+MeltingConstrainedFluxBalance() = MeltingConstrainedFluxBalance(NonlinearSurfaceTemperatureSolver())
+
+#####
+##### Flux imbalance and temperature
+#####
+
+@inline function top_flux_imbalance(i, j, grid, top_thermal_bc, top_temperature,
+                                    internal_fluxes, external_fluxes, clock, model_fields)
+
+    # Schematic of the Stefan condition at an upper top
+    #        
+    #  air       ↑   Qx ≡ external_fluxes. Example: Qx = σ T⁴
+    #          |⎴⎴⎴|
+    # ----------------------- ↕ hₛ(t) → ℒ ∂t hₛ = δQ | T ≤ Tₘ
+    #          |⎵⎵⎵|
+    #  ice       ↑   Qi ≡ internal_fluxes. Example Qi = - k ∂z T
+    #      
+
+    Qi = getflux(internal_fluxes, i, j, grid, top_temperature, clock, model_fields)
+    Qx = getflux(external_fluxes, i, j, grid, top_temperature, clock, model_fields)
+
+    # The imbalance is defined such that
+    # negative imbalance => heat accumulates (out > in) ⟹  growth
+    return Qx - Qi
+end
+
+@inline top_temperature(i, j, grid, ::PrescribedTemperature, Tu, args...) = Tu
+
+using RootSolvers: SecantMethod, find_zero, CompactSolution
+
+@inline function top_temperature(i, j, grid, top_thermal_bc,
+                                 current_top_temperature,
+                                 internal_fluxes, external_fluxes,
+                                 clock, model_fields)
+
+    Tu = @inbounds current_top_temperature[i, j, 1]
+    T₁ = Tu
+    T₂ = Tu - 1
+    FT = eltype(grid)
+    method = SecantMethod{FT}(T₁, T₂)
+    solution_type = CompactSolution()
+
+    flux_balance(T) = getflux(external_fluxes, i, j, grid, T, clock, model_fields) -
+                      getflux(internal_fluxes, i, j, grid, T, clock, model_fields)
+                      
+    solution = find_zero(flux_balance, method, solution_type)
+    
+    return solution.root
+end
+
+
+
diff --git a/test/runtests.jl b/test/runtests.jl
index 01d305c..59657e8 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,15 +1,15 @@
 using Test
 
-using ClimaSeaIce
 using Oceananigans
-using ClimaSeaIce.EulerianThermodynamicSeaIceModels: EulerianThermodynamicSeaIceModel, MolecularDiffusivity
+using ClimaSeaIce
+using ClimaSeaIce.EnthalpyMethodSeaIceModels: EnthalpyMethodSeaIceModel, MolecularDiffusivity
 
 κ = 1e-5
 grid = RectilinearGrid(size=3, z=(-3, 0), topology=(Flat, Flat, Bounded))
 closure = MolecularDiffusivity(grid, κ_ice=κ, κ_water=κ)
-model = EulerianThermodynamicSeaIceModel(; grid, closure)
+model = EnthalpyMethodSeaIceModel(; grid, closure)
 
-@test model isa EulerianThermodynamicSeaIceModel
+@test model isa EnthalpyMethodSeaIceModel
 
 # Test that it runs
 simulation = Simulation(model; Δt = 0.1 / κ, stop_iteration=3)