diff --git a/LICENSE b/LICENSE
index b5cf1c11..a2bdcc3a 100644
--- a/LICENSE
+++ b/LICENSE
@@ -1,6 +1,6 @@
 MIT License
 
-Copyright (c) 2022 Olav Møyner
+Copyright (c) 2022 Olav Møyner, SINTEF Digital and Contributors
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal
diff --git a/Project.toml b/Project.toml
index 137bf434..b9303555 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "JutulDarcy"
 uuid = "82210473-ab04-4dce-b31b-11573c4f8e0a"
 authors = ["Olav Møyner <olav.moyner@gmail.com>"]
-version = "0.2.32"
+version = "0.2.33"
 
 [deps]
 AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c"
@@ -54,7 +54,7 @@ DelimitedFiles = "1.6"
 DocStringExtensions = "0.9.3"
 ForwardDiff = "0.10.30"
 GLMakie = "0.10.0"
-GeoEnergyIO = "1.1.11"
+GeoEnergyIO = "1.1.12"
 HYPRE = "1.4.0"
 Jutul = "0.2.37"
 Krylov = "0.9.1"
diff --git a/docs/make.jl b/docs/make.jl
index 3b314d76..74363ddb 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -7,7 +7,13 @@ using DocumenterCitations
 using DocumenterVitepress
 ##
 cd(@__DIR__)
-function build_jutul_darcy_docs(build_format = nothing; build_examples = true, build_validation_examples = build_examples, build_notebooks = build_examples, clean = true)
+function build_jutul_darcy_docs(build_format = nothing;
+        build_examples = true,
+        build_validation_examples = build_examples,
+        build_notebooks = true,
+        clean = true,
+        deploy = true
+    )
     DocMeta.setdocmeta!(JutulDarcy, :DocTestSetup, :(using JutulDarcy; using Jutul); recursive=true)
     DocMeta.setdocmeta!(Jutul, :DocTestSetup, :(using Jutul); recursive=true)
     bib = CitationBibliography(joinpath(@__DIR__, "src", "refs.bib"))
@@ -26,6 +32,7 @@ function build_jutul_darcy_docs(build_format = nothing; build_examples = true, b
         "Intro: Compositional flow" => "co2_brine_2d_vertical",
         "Quarter-five-spot with variation" => "five_spot_ensemble",
         "Gravity circulation with CPR preconditioner" => "two_phase_unstable_gravity",
+        "CO2 property calculations" => "co2_props",
         "CO2 injection in saline aquifer" => "co2_sloped",
         "Compositional with five components" => "compositional_5components",
         "Parameter matching of Buckley-Leverett" => "optimize_simple_bl",
@@ -44,7 +51,7 @@ function build_jutul_darcy_docs(build_format = nothing; build_examples = true, b
         return content*"\n\n # ## Example on GitHub\n "*
         "# If you would like to run this example yourself, it can be downloaded from "*
         "the JutulDarcy.jl GitHub repository [as a script](https://github.com/sintefmath/JutulDarcy.jl/blob/main/examples/$pth.jl), "*
-        "or as a [Notebook](https://github.com/sintefmath/JutulDarcy.jl/blob/gh-pages/dev/$pth.ipynb)"
+        "or as a [Jupyter Notebook](https://github.com/sintefmath/JutulDarcy.jl/blob/gh-pages/dev/final_site/notebooks/$pth.ipynb)"
     end
     if clean
         for (ex, pth) in examples
@@ -81,25 +88,20 @@ function build_jutul_darcy_docs(build_format = nothing; build_examples = true, b
             upd(content) = update_footer(content, pth)
             Literate.markdown(in_pth, out_dir, preprocess = upd)
         end
-        if build_notebooks
-            Literate.notebook(in_pth, notebook_dir, execute = false)
-        end
     end
     ## Docs
     if isnothing(build_format)
-        if false
-            build_format = Documenter.HTML(;
-                prettyurls=get(ENV, "CI", "false") == "true",
-                canonical="https://sintefmath.github.io/JutulDarcy.jl",
-                edit_link="main",
-                size_threshold_ignore = ["ref/jutul.md", "docstrings.md"],
-                assets=String["assets/citations.css"],
-            )
-        else
-            build_format = DocumenterVitepress.MarkdownVitepress(
-                repo = "https://github.com/sintefmath/JutulDarcy.jl",
-            )
-        end
+        # Old Documenter code in case we want to go back.
+        # build_format = Documenter.HTML(;
+        #     prettyurls=get(ENV, "CI", "false") == "true",
+        #     canonical="https://sintefmath.github.io/JutulDarcy.jl",
+        #     edit_link="main",
+        #     size_threshold_ignore = ["ref/jutul.md", "docstrings.md"],
+        #     assets=String["assets/citations.css"],
+        # )
+        build_format = DocumenterVitepress.MarkdownVitepress(
+            repo = "https://github.com/sintefmath/JutulDarcy.jl",
+        )
     end
     makedocs(;
         modules=[JutulDarcy, Jutul],
@@ -113,10 +115,9 @@ function build_jutul_darcy_docs(build_format = nothing; build_examples = true, b
             "Introduction" => [
                 "JutulDarcy.jl" => "index.md",
                 "Getting started" =>"man/intro.md",
-                "References" => "extras/refs.md",
-                "FAQ" => "extras/faq.md",
-                "Jutul functions" => "ref/jutul.md"
-                ],
+                "Your first JutulDarcy.jl simulation" => "man/first_ex.md",
+                "FAQ" => "extras/faq.md"
+            ],
             "Manual" => [
                     "man/highlevel.md",
                     "man/basics/input_files.md",
@@ -128,24 +129,38 @@ function build_jutul_darcy_docs(build_format = nothing; build_examples = true, b
                     "man/basics/secondary.md",
                     "man/basics/parameters.md",
                     "man/basics/plotting.md",
+                    "man/basics/utilities.md",
                     ],
+            "Further reading" => [
+                "man/advanced/mpi.md",
+                "man/advanced/compiled.md",
+                "Jutul functions" => "ref/jutul.md",
+                "Bibliography" => "extras/refs.md"
+            ],
             "Examples: Introduction" => intros_markdown,
             "Examples: Usage" => examples_markdown,
-            "Examples: Validation" => validation_markdown,
-            "Advanced usage" => [
-                "man/advanced/mpi.md",
-                "man/advanced/compiled.md"
-            ]
+            "Examples: Validation" => validation_markdown
         ],
     )
-
-    deploydocs(;
-        repo="github.com/sintefmath/JutulDarcy.jl.git",
-        devbranch="main",
-        target = "build", # this is where Vitepress stores its output
-        branch = "gh-pages",
-        push_preview = true
-    )
+    if build_notebooks
+        # Subfolder of final site build folder
+        notebook_dir = joinpath(@__DIR__, "build", "final_site", "notebooks")
+        mkpath(notebook_dir)
+        for (ex, pth) in examples
+            in_pth = example_path(pth)
+            @info "$ex Writing notebook to $notebook_dir"
+            Literate.notebook(in_pth, notebook_dir, execute = false)
+        end
+    end
+    if deploy
+        deploydocs(;
+            repo="github.com/sintefmath/JutulDarcy.jl.git",
+            devbranch="main",
+            target = "build", # this is where Vitepress stores its output
+            branch = "gh-pages",
+            push_preview = true
+        )
+    end
 end
 ##
 # build_jutul_darcy_docs(build_examples = false, build_validation_examples = false)
diff --git a/docs/src/extras/faq.md b/docs/src/extras/faq.md
index f6dca166..d37933cf 100644
--- a/docs/src/extras/faq.md
+++ b/docs/src/extras/faq.md
@@ -1,8 +1,173 @@
 # Frequently asked questions
 
-!!! info "Note about units"
-    JutulDarcy does currently not make us of conversion factors or explicit
-    units can in principle use any consistent unit system. Some default scaling
-    of variables assume that the magnitude pressures and velocities roughly
-    match that of strict SI (e.g. Pascals and cubic meters per second). These
-    scaling factors are primarily used when iterative linear solvers are used.
+Here are a few common questions and possible answers.
+
+## Input and output
+
+### What input formats can JutulDarcy.jl use?
+
+1. DATA files (used by Eclipse, OPM Flow, tNavigator, Echelon and others) provided that the grid is given either as a corner-point GRDECL file or in TOPS format. As with most reservoir simulators, not all features of the original format are supported, but the code will let you know when unsupported features are encountered.
+2. Cases written out from [MRST](https://www.mrst.no/) through the `jutul` module.
+3. Cases written entirely in Julia using the basic `Jutul` and `JutulDarcy` data structures, as seen in the examples of the module.
+
+### What output formats does JutulDarcy.jl have?
+
+The simulator outputs results into standard Julia data structures (e.g. Vectors and Dicts) that can easily be written out using other Julia packages, for example in CSV format. We do not currently support binary formats output by commercial simulators.
+
+Simulation results are written to disk using [JLD2](https://github.com/JuliaIO/JLD2.jl), a subset of [HDF5](https://en.wikipedia.org/wiki/Hierarchical_Data_Format) commonly used in Julia for storing objects to disk.
+
+### How do I restart an interrupted simulation?
+
+JutulDarcy keeps everything in memory by default. This is not practical for larger models. If the argument `output_path` is set to a directory, JutulDarcy writes to the `JLD2` format (variant of HDF5).
+
+```julia
+# Note: set ENV["JUTUL_OUTPUT_PATH"] in your startup.jl first!
+pth = jutul_output_path("My_test_case")
+simulate_reservoir(case, output_path = pth)
+```
+
+If a output path is set, you can restart simulations:
+
+```julia
+# Restart from the last succesfully solved step, or return output if everything is simulated
+ws, states = simulate_reservoir(case, output_path = pth, restart = true)
+# Start from the beginning, overwriting files if already present
+ws, states = simulate_reservoir(case, output_path = pth, restart = false)
+# Restart from step 10 and throw error if step 9 is not already stored on disk.
+ws, states = simulate_reservoir(case, output_path = pth, restart = 10)
+```
+
+You can restart the simulation with different options for timestepping or tolerances.
+
+### How do I decide where output is stored?
+
+`Jutul.jl` contains a system for managing output folders. It is highly recommended that you amend your `startup.jl` file to include `ENV["JUTUL_OUTPUT_PATH"]` that points to where you want output to be stored. For example, on Windows usage of the output path mechanism may look something like this:
+
+```julia
+julia> ENV["JUTUL_OUTPUT_PATH"]
+"D:/jutul_output/"
+
+julia> jutul_output_path() # Randomly generated file name
+"D:/jutul_output/jutul/jl_DwpAvQTiLo"
+
+julia> jutul_output_path("mycase")
+"D:/jutul_output/jutul/mycase"
+
+julia> jutul_output_path("mycase", subfolder = "ensemble_name")
+"D:/jutul_output/ensemble_name/mycase"
+
+julia> jutul_output_path("mycase", subfolder = missing)
+"D:/jutul_output/mycase"
+```
+
+Or equivialent on a Linux system:
+
+```julia
+julia> ENV["JUTUL_OUTPUT_PATH"]
+"/home/username/jutul_output/"
+
+julia> jutul_output_path() # Randomly generated file name
+"/home/username/jutul_output/jutul/jl_DwpAvQTiLo"
+
+julia> jutul_output_path("mycase")
+"/home/username/jutul_output/jutul/mycase"
+
+julia> jutul_output_path("mycase", subfolder = "ensemble_name")
+"/home/username/jutul_output/ensemble_name/mycase"
+
+julia> jutul_output_path("mycase", subfolder = missing)
+"/home/username/jutul_output/mycase"
+```
+
+You can also just specify a full path and keep track of output folders yourself, but using the `jutul_output_path` mechanism will make it easier to write a script that can be run on another computer with different folder structure.
+
+### How do you get out more output from a simulation?
+
+The default outputs per cell are primary variables and total masses:
+
+```julia
+reservoir_model(model).output_variables
+3-element Vector{Symbol}:
+ :Pressure
+ :Saturations
+ :TotalMasses
+```
+
+You can push variables to this list, or ask the code to output all variables:
+
+```julia
+model2, = setup_reservoir_model(domain, sys, extra_outputs = true);
+reservoir_model(model2).output_variables
+7-element Vector{Symbol}:
+ :Pressure
+ :Saturations
+ :TotalMasses
+ :PhaseMassDensities
+ :RelativePermeabilities
+ :PhaseMobilities
+ :PhaseMassMobilities
+```
+
+Can also pass `extra_outputs = [:PhaseMobilities]`  to output specific variables.
+
+## Plotting
+
+### What is required for visualization?
+
+We use the wonderful [`Makie.jl`](https://docs.makie.org/) for both 2D and 3D plots. Generally `CairoMakie` is supported for non-interactive plotting and `GLMakie` is used for interactive plotting (especially 3D). The latter requires a working graphics context, which is not directly available when the code is run over for example SSH or on a server.
+
+For more details on the backends, see the[`Makie.jl` docs](https://docs.makie.org/stable/explanations/backends/backends)
+
+## Miscellaneous
+
+### Can you add feature X or format Y?
+
+If you have a feature you'd like to have supported, please file an
+[issue](https://github.com/JuliaIO/JLD2.jl/issues) with details on the format.
+`JutulDarcy` is developed primarily through contract research, so features are
+added as needed for ongoing projects where the simulator is in use. Posting an
+issue, especially if you have a clear reference to how something should be
+implemented is still very useful. It is also possible to fund the development
+for a specific feature, or to implement the feature yourself by asking for
+pointers on how to get started.
+
+### What units does JutulDarcy.jl use?
+
+JutulDarcy uses consistent units. This typically means that all your values must
+be input in strict SI. This means pressures in Pascal, temperatures in Kelvin
+and time in seconds. Note that this is very similar to the `METRIC` type of unit
+system seen in many commercial simulators, except that units of time is not
+given in days. This also impacts permeabilities, transmissibilities and
+viscosities, which will be much smaller than in metric where days are used.
+
+Reading of input files will automatically convert data to the correct units for simulation, but care must be taken when you are writing your own code. `Jutul.jl` contains unit conversion factors to make it easier to write code:
+
+```@example
+using Jutul
+p = convert_to_si(120.0, :bar)
+```
+
+You can also extract individual units and to the setup yourself:
+
+```@example
+using Jutul
+day, stb = si_units(:day, :stb)
+# convert to m^3/s:
+rate = 100.0stb/day
+```
+
+JutulDarcy does currently not make use of conversion factors or explicit
+units can in principle use any consistent unit system. Some default scaling
+of variables assume that the magnitude pressures and velocities roughly
+match that of strict SI (e.g. Pascals and cubic meters per second). These
+scaling factors are primarily used when iterative linear solvers are used.
+
+### Why write a new reservoir simulation code in Julia?
+
+We believe that reservoir simulation should be a *library* and not necessarily an application by itself. The future of porous media simulation is deeply integrated into other workflows, and not as an application that simply writes files to disk. As a part of experimentation in differentiable and flexible solvers using automatic differentiation that started with [MRST](https://www.mrst.no), Julia was the natural next step.
+
+### What is the license of JutulDarcy.jl?
+
+The code uses the [MIT license](https://en.wikipedia.org/wiki/MIT_License), a permissive license that requires attribution, but does not place limitations on commercial use or closed-source integration.
+
+The code uses a number of dependencies that can have other licenses and we make no guarantees that the entirety of the code made available by adding `JutulDarcy.jl` to a given Julia environment is all MIT licensed.
diff --git a/docs/src/extras/refs.md b/docs/src/extras/refs.md
index 78f29bc4..e4c97813 100644
--- a/docs/src/extras/refs.md
+++ b/docs/src/extras/refs.md
@@ -1,4 +1,2 @@
-# References
-
 ```@bibliography
 ```
diff --git a/docs/src/man/advanced/compiled.md b/docs/src/man/advanced/compiled.md
index a1481348..df5d51d1 100644
--- a/docs/src/man/advanced/compiled.md
+++ b/docs/src/man/advanced/compiled.md
@@ -1,3 +1,21 @@
 # Standalone reservoir simulator
 
-It is possible build a standalone reservoir simulator using [PackageCompiler.jl](https://github.com/JuliaLang/PackageCompiler.jl). For more details and an example build file with keyword arguments, see [the JutulDarcyApps.jl repository](https://github.com/sintefmath/JutulDarcyApps.jl/tree/master/mpi_simulator).
+Scripts are interactive and useful for doing setup, simulation and post-processing in one file, but sometimes you want to run a big model unmodified from an input file:
+
+- As an alternative to a pure Julia workflow, `JutulDarcy.jl` can be compiled into a standalone reservoir simulator
+- This makes MPI simulations more ergonomic
+- Compiling the code saves time when running multiple simulations
+- The resulting executable is a standard command-line program - no Julia experience needed
+- Output is given in the same format as regular simulations, can load data by restarting a simulation from the same `output_path`
+
+This workflow uses [PackageCompiler.jl](https://github.com/JuliaLang/PackageCompiler.jl). For more details and an example build file with keyword arguments, see [the JutulDarcyApps.jl repository](https://github.com/sintefmath/JutulDarcyApps.jl/tree/master/mpi_simulator).
+
+A few things to note:
+
+- The simulator comes with a set of shared library files and will be ~500 mb
+- Binaries will match platform (compiling under Linux gives you Linux binaries)
+- The repository has a script that runs small "representative" models
+- You can input small representative models in `precompile_jutul_darcy_mpi.jl` to make sure that compilation is avoided during simulation
+- By default, the script uses the default Julia MPI binary. On a cluster, the build script may have to be modified to use the MPI type of the cluster using [MPITrampoline.jl](https://github.com/eschnett/MPItrampoline)
+
+If you get it working on a complex MPI setup, feedback on your experience and PRs are very welcome.
diff --git a/docs/src/man/advanced/mpi.md b/docs/src/man/advanced/mpi.md
index dcc75a96..7e92cd32 100644
--- a/docs/src/man/advanced/mpi.md
+++ b/docs/src/man/advanced/mpi.md
@@ -16,15 +16,40 @@ JutulDarcy also supports threads. By defualt, this only parallelizes property ev
 
 An experimental thread-parallel backend for matrices and linear algebra can be enabled by setting `backend=:csr` in the call to [`setup_reservoir_model`](@ref). This backend provides additional features such as a parallel zero-overlap ILU(0) implementation and parallel apply for AMG, but these features are still work in progress.
 
+Starting Julia with multiple threads (for example `julia --project. --threads=4`) will allow `JutulDarcy` to make use of threads to speed up calculations
+
+- The default behavior is to only speed up assembly of equations
+- The linear solver is often the most expensive part -- as mentioned above, parts can be parallelized by choosing `csr` backend when setting up the model
+- Running with a parallel preconditioner can lead to higher iteration counts since the ILU(0) preconditioner changes in parallel
+- Heavy compositional models benefit a lot from using threads
+
+Threads are easy to use and can give a bit of benefit for large models.
+
+
 ### Mixed-mode parallelism
 
 You can mix the two approaches: Adding multiple threads to each MPI process can use threads to speed up assembly and property evaluations.
 
+### Tips for parallel runs
+
+A few hints when you are looking at performance:
+
+- Reservoir simulations are memory bound, cannot expect that 10 threads = 10x performance
+- CPUs can often boost single-core performance when resources are available
+- MPI in JutulDarcy is less tested than single-process simulations, but is natural for larger models
+- There is always some cost to parallelism: If running a large ensemble with limited compute, many serial runs handled by Julia's task system is usually a better option
+- Adding the maximum number of processes does not always give the best performance. Typically you want at least 10 000 cells per process. Can be case dependent.
+
+Example: 200k cell model on laptop: 1 process 235 s -> 4 processes 145s
+
+
 ## Solving with MPI in practice
 
+There are a few adjustments needed before a script can be run in MPI.
+
 ### Setting up the environment
 
-Tou will have to set up an environment with the following packages under Julia 1.9+:
+You will have to set up an environment with the following packages under Julia 1.9+:
 `PartitionedArrays`, `MPI`, `JutulDarcy` and `HYPRE`. This is generally the best performing solver setup available, even if you are working in a shared memory environment.
 
 ### Writing the script
@@ -37,14 +62,47 @@ Write your script as usual. The following options must then be set:
 You must then run the file using the approprioate `mpiexec` as described in the MPI.jl documentation. Specialized functions will be called by `simulate_reservoir` when this is the case. We document them here, even if we recommend using the high level version of this interface:
 
 ```@docs
-setup_reservoir_simulator_parray
-simulate_reservoir_parray
+JutulDarcy.setup_reservoir_simulator_parray
+JutulDarcy.simulate_reservoir_parray
 ```
 
-### Limitations for running in MPI
+### Checklist for running in MPI
 
-!!! note
-    You should be familiar with the MPI programming model to use this feature. See [MPI.jl](https://juliaparallel.org/MPI.jl/stable/) for more details, and how MPI is handled in Julia specifically.
+- Install and load the following packages at the top of your script: `PartitionedArrays, MPI, HYPRE`
+- (Recommended): Put `MPI.Init()` at the top of your script
+- Make sure that `split_wells = true` is set in your model setup
+- Set `output_path` when running the simulation (otherwise the results will not be stored anywhere)
+- Set `mode=:mpi` when running the simulation.
+
+A few useful functions:
+
+- `MPI.install_mpiexecjl()` installs MPI that works "out of the box" with your current Julia setup
+- `MPI.mpiexec()` gives you the path to the executable and all environment variables
+
+A typical command to launch a MPI script from within Julia:
+
+```julia
+n = 5 # = 5 processes
+script_to_run = "my_script.jl"
+run(`$(mpiexec()) -n $n $(Base.julia_cmd()) --project=$(Base.active_project()) $script_to_run`)
+```
+
+Adding threads to the command will make JutulDarcy use both threads and processes
+
+:warning: **Running a script in MPI means that all parts of the script will run on each process!** If you want to do data analysis you will have to either wrap your code in `if MPI.Comm_rank(MPI.COMM_WORLD) == 0` or do the data analysis in serial (recommended).
+
+### Limitations of running in MPI
+
+MPI can be cumbersome to use when compared to a standard Julia script, and the
+current implementation relies on the model being set up on each processor before
+subdivision. This can be quite memory intensive during startup.
+
+You should be familiar with the MPI programming model to use this feature. See
+[MPI.jl](https://juliaparallel.org/MPI.jl/stable/) for more details, and how MPI
+is handled in Julia specifically.
+
+For larger models, compiling the [Standalone reservoir simulator](@ref) is
+highly recommended.
 
 !!! note
     MPI consolidates results by writing files to disk. Unless you have a plan to work with the distributed states in-memory returned by the `simulate!` call, it is best to specify a `output_path` optional argument to [`setup_reservoir_simulator`](@ref). After the simulation, that folder will contain output just as if you had run the case in serial.
diff --git a/docs/src/man/basics/forces.md b/docs/src/man/basics/forces.md
index a55ea751..6f00597d 100644
--- a/docs/src/man/basics/forces.md
+++ b/docs/src/man/basics/forces.md
@@ -1,13 +1,19 @@
 # Driving forces
 
+```@docs
+JutulDarcy.setup_reservoir_forces
+```
+
 ## Source terms
 
 ```@docs
 SourceTerm
+JutulDarcy.FlowSourceType
 ```
 
 ## Boundary conditions
 
 ```@docs
 FlowBoundaryCondition
+flow_boundary_condition
 ```
diff --git a/docs/src/man/basics/input_files.md b/docs/src/man/basics/input_files.md
index 74f5cdc3..483d1b74 100644
--- a/docs/src/man/basics/input_files.md
+++ b/docs/src/man/basics/input_files.md
@@ -28,6 +28,7 @@ If you want to directly simulate a file from disk, you can sue the high level fu
 simulate_data_file
 setup_case_from_data_file
 JutulDarcy.setup_case_from_parsed_data
+JutulDarcy.convert_co2store_to_co2_brine
 ```
 
 Reservoir simulator input files are highly complex and contain a great number of
diff --git a/docs/src/man/basics/parameters.md b/docs/src/man/basics/parameters.md
index cc15c1d0..7fb134da 100644
--- a/docs/src/man/basics/parameters.md
+++ b/docs/src/man/basics/parameters.md
@@ -32,3 +32,9 @@ JutulDarcy.EndPointScalingCoefficients
 JutulDarcy.WellIndices
 JutulDarcy.PerforationGravityDifference
 ```
+
+## Thermal
+
+```@docs
+JutulDarcy.FluidThermalConductivities
+```
\ No newline at end of file
diff --git a/docs/src/man/basics/plotting.md b/docs/src/man/basics/plotting.md
index 4206fc36..8bdb402b 100644
--- a/docs/src/man/basics/plotting.md
+++ b/docs/src/man/basics/plotting.md
@@ -6,4 +6,5 @@ plot_reservoir
 plot_well_results
 JutulDarcy.plot_reservoir_simulation_result
 Jutul.plot_cell_data
+JutulDarcy.plot_well!
 ```
diff --git a/docs/src/man/basics/secondary.md b/docs/src/man/basics/secondary.md
index 6da29c2e..b38a0350 100644
--- a/docs/src/man/basics/secondary.md
+++ b/docs/src/man/basics/secondary.md
@@ -9,9 +9,20 @@
 ```@docs
 BrooksCoreyRelativePermeabilities
 RelativePermeabilities
+JutulDarcy.TabulatedSimpleRelativePermeabilities
 JutulDarcy.ReservoirRelativePermeabilities
 ```
 
+The `ReservoirRelativePermeabilities` type also supports hysteresis for either phase.
+
+```@docs
+JutulDarcy.NoHysteresis
+JutulDarcy.CarlsonHysteresis
+JutulDarcy.KilloughHysteresis
+JutulDarcy.JargonHysteresis
+JutulDarcy.ImbibitionOnlyHysteresis
+```
+
 ```@docs
 PhaseRelativePermeability
 ```
@@ -24,16 +35,20 @@ EndPointScalingCoefficients
 
 ```@docs
 DeckPhaseViscosities
+JutulDarcy.ConstMuBTable
+JutulDarcy.MuBTable
 ```
 
 ```@docs
 PhaseMassDensities
 ```
 
-### Immiscible flow
-
 #### Phase densities
 
+```@docs
+JutulDarcy.DeckPhaseMassDensities
+```
+
 #### Shrinkage factors
 
 ```@docs
@@ -53,6 +68,7 @@ ConstantCompressibilityDensities
 
 ```@docs
 PhaseMassFractions
+JutulDarcy.KValueWrapper
 ```
 
 ## Wells
diff --git a/docs/src/man/basics/solution.md b/docs/src/man/basics/solution.md
index f1497e19..39c06075 100644
--- a/docs/src/man/basics/solution.md
+++ b/docs/src/man/basics/solution.md
@@ -69,7 +69,7 @@ This block system has several advantages:
 
 ##### Constrained Pressure Residual
 
-The CPR preconditioner [wallis-cpr, cao-cpr](@cite) [`CPRPreconditioner`](@ref) is a multi-stage physics-informed preconditioner that seeks to decouple the global pressure part of the system from the local  transport part. In the limits of incompressible flow without gravity it can be thought of as an elliptic / hyperbolic splitting.
+The CPR preconditioner [wallis-cpr, cao-cpr](@cite) [`CPRPreconditioner`](@ref) is a multi-stage physics-informed preconditioner that seeks to decouple the global pressure part of the system from the local  transport part. In the limits of incompressible flow without gravity it can be thought of as an elliptic / hyperbolic splitting. We also implement a special variant for the adjoint system that is similar to the treatment described in [adjoint_cpr](@cite).
 
 ```@docs
 CPRPreconditioner
diff --git a/docs/src/man/basics/systems.md b/docs/src/man/basics/systems.md
index e6349481..a769f18b 100644
--- a/docs/src/man/basics/systems.md
+++ b/docs/src/man/basics/systems.md
@@ -157,7 +157,12 @@ For additional details, please see Chapter 8 - Compositional Simulation with the
 
 ## Multi-phase thermal flow
 
+Thermal effects are modelled as an additional residual equation that comes in addition to the flow equations.
+
+```math
+r_i = \frac{\partial}{\partial t}\Bigl(\rho_r U_r (1-\phi) + \sum_\alpha \rho_\alpha S_\alpha U_\alpha \phi \Bigr) + \nabla \cdot \left (\sum_\alpha ( H_\alpha \rho_\alpha v_\alpha - S_\alpha \lambda_\alpha \nabla T)-\lambda_r \nabla T \right) - Q_e
+```
+
 ```@docs
-reservoir_system
 JutulDarcy.add_thermal_to_model!
 ```
diff --git a/docs/src/man/basics/utilities.md b/docs/src/man/basics/utilities.md
new file mode 100644
index 00000000..6fd0e66c
--- /dev/null
+++ b/docs/src/man/basics/utilities.md
@@ -0,0 +1,50 @@
+# Utilities
+
+This section describes various utilities that do not fit in other sections.
+
+## CO2 and brine correlations
+
+These functions are not exported, but can be found inside the `CO2Properties` submodule. The functions described here are a Julia port of the MRST module described in [salo_co2](@cite) They can be accessed by explicit import:
+
+```julia
+import JutulDarcy.CO2Properties: name_of_function
+```
+
+```@docs
+JutulDarcy.CO2Properties.compute_co2_brine_props
+JutulDarcy.CO2Properties.pvt_brine_RoweChou1970
+JutulDarcy.CO2Properties.activity_co2_DS2003
+JutulDarcy.CO2Properties.viscosity_brine_co2_mixture_IC2012
+JutulDarcy.CO2Properties.pvt_brine_BatzleWang1992
+JutulDarcy.CO2Properties.viscosity_co2_Fenghour1998
+JutulDarcy.CO2Properties.pvt_co2_RedlichKwong1949
+JutulDarcy.CO2Properties.viscosity_gas_mixture_Davidson1993
+```
+
+## Relative permeability functions
+
+```@docs
+JutulDarcy.brooks_corey_relperm
+```
+
+## CO2 inventory
+
+```@docs
+JutulDarcy.co2_inventory
+JutulDarcy.plot_co2_inventory
+```
+
+## API utilities
+
+```@docs
+reservoir_model
+JutulDarcy.reservoir_storage
+JutulDarcy.well_symbols
+```
+
+## Adjoints and gradients
+
+```@docs
+JutulDarcy.reservoir_sensitivities
+JutulDarcy.well_mismatch
+```
diff --git a/docs/src/man/basics/wells.md b/docs/src/man/basics/wells.md
index e13c1b44..ecdb6572 100644
--- a/docs/src/man/basics/wells.md
+++ b/docs/src/man/basics/wells.md
@@ -23,6 +23,8 @@ JutulDarcy.SimpleWell
 ```@docs
 WellDomain
 MultiSegmentWell
+JutulDarcy.SegmentWellBoreFrictionHB
+JutulDarcy.PotentialDropBalanceWell
 ```
 
 ## Well controls and limits
@@ -53,16 +55,30 @@ TotalRateTarget
 HistoricalReservoirVoidageTarget
 ReservoirVoidageTarget
 DisabledTarget
+JutulDarcy.TotalReservoirRateTarget
+```
+
+### Implementation of well controls
+
+```@docs
+JutulDarcy.well_target
+```
+
+### Well outputs
+
+```@docs
+JutulDarcy.print_well_result_table
 ```
 
 ### Imposing limits on wells (multiple constraints)
 
 ## Well forces
 
-### Changing perforations
+### Perforations and WI adjustments
 
 ```@docs
 PerforationMask
+JutulDarcy.Perforations
 ```
 
 ### Other forces
diff --git a/docs/src/man/first_ex.md b/docs/src/man/first_ex.md
new file mode 100644
index 00000000..c3b51a4d
--- /dev/null
+++ b/docs/src/man/first_ex.md
@@ -0,0 +1,167 @@
+# Your first JutulDarcy.jl simulation
+
+After a bit of time has passed compiling the packages, you are now ready to use JutulDarcy. There are a number of examples included in this manual, but we include a brief example here that briefly demonstrates the key concepts.
+
+## Setting up the domain
+
+We set up a simple Cartesian Mesh that is converted into a reservoir domain with
+static properties permeability and porosity values together with geometry
+information. We then use this domain to set up two wells: One vertical well for
+injection and a single perforation producer.
+
+````@example intro_ex
+using JutulDarcy, Jutul
+Darcy, bar, kg, meter, day = si_units(:darcy, :bar, :kilogram, :meter, :day)
+nx = ny = 25
+nz = 10
+cart_dims = (nx, ny, nz)
+physical_dims = (1000.0, 1000.0, 100.0).*meter
+g = CartesianMesh(cart_dims, physical_dims)
+domain = reservoir_domain(g, permeability = 0.3Darcy, porosity = 0.2)
+Injector = setup_vertical_well(domain, 1, 1, name = :Injector)
+Producer = setup_well(domain, (nx, ny, 1), name = :Producer)
+# Show the properties in the domain
+domain
+````
+
+## Setting up a fluid system
+
+We select a two-phase immicible system by declaring that the liquid and vapor
+phases are present in the model. These are assumed to have densities of 1000 and
+100 kilograms per meters cubed at reference pressure and temperature conditions.
+
+````@example intro_ex
+phases = (LiquidPhase(), VaporPhase())
+rhoLS = 1000.0kg/meter^3
+rhoGS = 100.0kg/meter^3
+reference_densities = [rhoLS, rhoGS]
+sys = ImmiscibleSystem(phases, reference_densities = reference_densities)
+````
+
+## Setting up the model
+
+We now have everything we need to set up a model. We call the reservoir model
+setup function and get out the model together with the parameters. Parameter
+represent numerical input values that are static throughout the simulation.
+These are automatically computed from the domain's geometry, permeability and
+porosity.
+
+````@example intro_ex
+model, parameters = setup_reservoir_model(domain, sys, wells = [Injector, Producer])
+model
+````
+
+The model has a set of default secondary variables (properties) that are used to compute the flow equations. We can have a look at the reservoir model to see what the defaults are for the Darcy flow part of the domain:
+
+````@example intro_ex
+reservoir_model(model)
+````
+
+The secondary variables can be swapped out, replaced and new variables can be added with arbitrary functional dependencies thanks to Jutul's flexible setup for automatic differentiation. Let us adjust the defaults by replacing the relative permeabilities with Brooks-Corey functions and the phase density functions by constant compressibilities:
+
+````@example intro_ex
+c = [1e-6, 1e-4]/bar
+density = ConstantCompressibilityDensities(
+    p_ref = 100*bar,
+    density_ref = reference_densities,
+    compressibility = c
+)
+kr = BrooksCoreyRelativePermeabilities(sys, [2.0, 3.0])
+replace_variables!(model, PhaseMassDensities = density, RelativePermeabilities = kr)
+nothing #hide
+````
+
+## Initial state: Pressure and saturations
+
+Now that we are happy with our model setup, we can designate an initial state.
+Equilibriation of reservoirs can be a complicated affair, but here we set up a
+constant pressure reservoir filled with the liquid phase. The inputs must match
+the primary variables of the model, which in this case is pressure and
+saturations in every cell.
+
+````@example intro_ex
+state0 = setup_reservoir_state(model,
+    Pressure = 120bar,
+    Saturations = [1.0, 0.0]
+)
+````
+
+## Setting up timesteps and well controls
+
+We set up reporting timesteps. These are the intervals that the simulator gives out outputs. The simulator may use shorter steps internally, but will always hit these points in the output. Here, we report every year for 25 years.
+
+````@example intro_ex
+nstep = 25
+dt = fill(365.0day, nstep);
+````
+
+We next set up a rate target with a high amount of gas injected into the model. This is not fully realistic, but will give some nice and dramatic plots for our example later on.
+
+````@example intro_ex
+pv = pore_volume(model, parameters)
+inj_rate = 1.5*sum(pv)/sum(dt)
+rate_target = TotalRateTarget(inj_rate)
+````
+
+The producer is set to operate at a fixed pressure:
+
+````@example intro_ex
+bhp_target = BottomHolePressureTarget(100bar)
+````
+
+We can finally set up forces for the model. Note that while JutulDarcy supports time-dependent forces and limits for the wells, we keep this example as simple as possible.
+
+````@example intro_ex
+I_ctrl = InjectorControl(rate_target, [0.0, 1.0], density = rhoGS)
+P_ctrl = ProducerControl(bhp_target)
+controls = Dict(:Injector => I_ctrl, :Producer => P_ctrl)
+forces = setup_reservoir_forces(model, control = controls);
+````
+
+## Simulate and analyze results
+
+We call the simulation with our initial state, our model, the timesteps, the forces and the parameters:
+
+````@example intro_ex
+wd, states, t = simulate_reservoir(state0, model, dt, parameters = parameters, forces = forces)
+````
+
+We can interactively look at the well results in the command line:
+
+````@example intro_ex
+wd(:Producer)
+````
+
+Let us look at the pressure evolution in the injector:
+
+````@example intro_ex
+wd(:Injector, :bhp)
+````
+
+If we have a plotting package available, we can visualize the results too:
+
+````@example intro_ex
+using GLMakie
+grat = wd[:Producer, :grat]
+lrat = wd[:Producer, :lrat]
+bhp = wd[:Injector, :bhp]
+fig = Figure(size = (1200, 400))
+ax = Axis(fig[1, 1],
+    title = "Injector",
+    xlabel = "Time / days",
+    ylabel = "Bottom hole pressure / bar")
+lines!(ax, t/day, bhp./bar)
+ax = Axis(fig[1, 2],
+    title = "Producer",
+    xlabel = "Time / days",
+    ylabel = "Production rate / m³/day")
+lines!(ax, t/day, abs.(grat).*day)
+lines!(ax, t/day, abs.(lrat).*day)
+fig
+````
+
+Interactive visualization of the 3D results is also possible if GLMakie is loaded:
+
+````@example intro_ex
+plot_reservoir(model, states, key = :Saturations, step = 3)
+````
diff --git a/docs/src/man/highlevel.md b/docs/src/man/highlevel.md
index 6fbe7473..2b180ebb 100644
--- a/docs/src/man/highlevel.md
+++ b/docs/src/man/highlevel.md
@@ -2,6 +2,15 @@
 
 ## Setup
 
+The basic outline of building a reservoir simulation problem consists of:
+
+1. Making a mesh
+2. Converting the mesh into a reservoir, adding properties
+3. Add any number of wells
+4. Setup a physical system and setup a reservoir model
+5. Set up timesteps, well controls and other forces
+6. Simulate!
+
 ### Meshes
 
 JutulDarcy can use meshes that supported by Jutul. This includes the Cartesian and Unstructured meshes as well as any meshes in the more general [Meshes.jl](https://github.com/JuliaGeometry/Meshes.jl) package.
diff --git a/docs/src/man/intro.md b/docs/src/man/intro.md
index f61e3168..4251e877 100644
--- a/docs/src/man/intro.md
+++ b/docs/src/man/intro.md
@@ -19,7 +19,7 @@ using Pkg
 Pkg.add("JutulDarcy")
 ```
 
-You can then run any of the examples in the [`examples`](https://github.com/sintefmath/JutulDarcy.jl/tree/main/examples) directory by including them. The examples are sorted by complexity. We suggest you start with [Gravity segregation example](Intro: Gravity segregation example).
+You can then run any of the examples in the [`examples`](https://github.com/sintefmath/JutulDarcy.jl/tree/main/examples) directory by including them. The examples are sorted by complexity. We suggest you start with [Intro: Gravity segregation](Intro: Gravity segregation example).
 
 To generate a folder that contains the examples locally, you can run the following code to create a folder `jutuldarcy_examples` in your current working directory:
 
@@ -35,6 +35,10 @@ using JutulDarcy
 generate_jutuldarcy_examples("/home/username/")
 ```
 
+```@docs
+generate_jutuldarcy_examples
+```
+
 ### Adding additional packages
 
 We also rely heavily on the Jutul base package for some functionality, so we recommend that you also install it together with JutulDarcy:
@@ -56,160 +60,3 @@ Pkg.add("Optim") # Optimization library
 Pkg.add("HYPRE") # Better linear solver
 Pkg.add("GeoEnergyIO") # Parsing input files
 ```
-
-## A short motivational example with two wells
-
-After a bit of time has passed compiling the packages, you are now ready to use JutulDarcy. There are a number of examples included in this manual, but we include a brief example here that briefly demonstrates the key concepts.
-
-### Setting up the domain
-
-We set up a simple Cartesian Mesh that is converted into a reservoir domain with permeability and porosity values. We then use this domain to set up two wells: One vertical well for injection and a single perforation producer:
-
-````@example intro_ex
-using JutulDarcy, Jutul
-Darcy, bar, kg, meter, day = si_units(:darcy, :bar, :kilogram, :meter, :day)
-nx = ny = 25
-nz = 10
-cart_dims = (nx, ny, nz)
-physical_dims = (1000.0, 1000.0, 100.0).*meter
-g = CartesianMesh(cart_dims, physical_dims)
-domain = reservoir_domain(g, permeability = 0.3Darcy, porosity = 0.2)
-Injector = setup_vertical_well(domain, 1, 1, name = :Injector)
-Producer = setup_well(domain, (nx, ny, 1), name = :Producer)
-# Show the properties in the domain
-domain
-````
-
-### Setting up a fluid system
-
-We select a two-phase immicible system by declaring that the liquid and vapor phases are present in the model. These are assumed to have fixed densitys of 1000 and 100 kilograms per meters cubed at some reference pressure and temperature conditions.
-
-````@example intro_ex
-phases = (LiquidPhase(), VaporPhase())
-rhoLS = 1000.0kg/meter^3
-rhoGS = 100.0kg/meter^3
-reference_densities = [rhoLS, rhoGS]
-sys = ImmiscibleSystem(phases, reference_densities = reference_densities)
-````
-
-### Setting up the model
-
-We now have everything we need to set up a model. We call the setup function and get out the model together with the parameters - numerical input values that are static throughout the simulation. These are automatically computed from the domain's geometry, permeability and porosity.
-
-````@example intro_ex
-model, parameters = setup_reservoir_model(domain, sys, wells = [Injector, Producer])
-model
-````
-
-The model has a set of default secondary variables (properties) that are used to compute the flow equations. We can have a look at the reservoir model to see what the defaults are for the Darcy flow part of the domain:
-
-````@example intro_ex
-reservoir_model(model)
-````
-
-The secondary variables can be swapped out, replaced and new variables can be added with arbitrary functional dependencies thanks to Jutul's flexible setup for automatic differentiation. Let us adjust the defaults by replacing the relative permeabilities with Brooks-Corey functions and the phase density functions by constant compressibilities:
-
-````@example intro_ex
-c = [1e-6, 1e-4]/bar
-density = ConstantCompressibilityDensities(
-    p_ref = 100*bar,
-    density_ref = reference_densities,
-    compressibility = c
-)
-kr = BrooksCoreyRelativePermeabilities(sys, [2.0, 3.0])
-replace_variables!(model, PhaseMassDensities = density, PhaseRelativePermeability = kr)
-nothing #hide
-````
-
-### Initial state
-
-Now that we are happy with our model setup, we can designate an initial state. For simplicity, we set up a constant pressure reservoir filled with the liquid phase.
-
-````@example intro_ex
-state0 = setup_reservoir_state(model,
-    Pressure = 120bar,
-    Saturations = [1.0, 0.0]
-)
-````
-
-### Setting up timesteps and well controls
-
-We set up reporting timesteps. These are the intervals that the simulator gives out outputs. The simulator may use shorter steps internally, but will always hit these points in the output.
-
-````@example intro_ex
-nstep = 25
-dt = fill(365.0day, nstep)
-nothing #hide
-````
-
-We next set up a rate target with a high amount of gas injected into the model. This is not fully realistic, but will give some nice and dramatic plots for our example later on.
-
-````@example intro_ex
-pv = pore_volume(model, parameters)
-inj_rate = 1.5*sum(pv)/sum(dt)
-rate_target = TotalRateTarget(inj_rate)
-````
-
-The producer is set to operate at a fixed pressure:
-
-````@example intro_ex
-bhp_target = BottomHolePressureTarget(100bar)
-````
-
-We can finally set up forces for the model. Note that while JutulDarcy supports time-dependent forces and limits for the wells, we keep this example as simple as possible.
-
-````@example intro_ex
-I_ctrl = InjectorControl(rate_target, [0.0, 1.0], density = rhoGS)
-P_ctrl = ProducerControl(bhp_target)
-controls = Dict(:Injector => I_ctrl, :Producer => P_ctrl)
-forces = setup_reservoir_forces(model, control = controls)
-nothing #hide
-````
-
-### Simulate and analyze results
-
-We call the simulation with our initial state, our model, the timesteps, the forces and the parameters:
-
-````@example intro_ex
-wd, states, t = simulate_reservoir(state0, model, dt, parameters = parameters, forces = forces)
-````
-
-We can interactively look at the well results in the command line:
-
-````@example intro_ex
-wd(:Producer)
-````
-
-Let us look at the pressure evolution in the injector:
-
-````@example intro_ex
-wd(:Injector, :bhp)
-````
-
-If we have a plotting package available, we can visualize the results too:
-
-````@example intro_ex
-using GLMakie
-grat = wd[:Producer, :grat]
-lrat = wd[:Producer, :lrat]
-bhp = wd[:Injector, :bhp]
-fig = Figure(size = (1200, 400))
-ax = Axis(fig[1, 1],
-    title = "Injector",
-    xlabel = "Time / days",
-    ylabel = "Bottom hole pressure / bar")
-lines!(ax, t/day, bhp./bar)
-ax = Axis(fig[1, 2],
-    title = "Producer",
-    xlabel = "Time / days",
-    ylabel = "Production rate / m³/day")
-lines!(ax, t/day, abs.(grat).*day)
-lines!(ax, t/day, abs.(lrat).*day)
-fig
-````
-
-Interactive visualization of the 3D results is also possible if GLMakie is loaded:
-
-````@example intro_ex
-plot_reservoir(model, states, key = :Saturations, step = 3)
-````
diff --git a/docs/src/ref/jutul.md b/docs/src/ref/jutul.md
index febabc42..72886ea7 100644
--- a/docs/src/ref/jutul.md
+++ b/docs/src/ref/jutul.md
@@ -1,6 +1,8 @@
 # Documentation from Jutul.jl
 
-JutulDarcy.jl builds upon Jutul.jl. You can use JutulDarcy.jl without knowing the inner workings of Jutul.jl, but if you want to dive under the hood these docstrings may be helpful.
+`JutulDarcy.jl` builds upon `Jutul.jl`. You can use `JutulDarcy.jl` without
+knowing the inner workings of `Jutul.jl`, but if you want to dive under the hood
+these docstrings may be helpful.
 
 ```@autodocs
 Modules = [Jutul]
diff --git a/docs/src/refs.bib b/docs/src/refs.bib
index 1e697d95..6f6d00a5 100644
--- a/docs/src/refs.bib
+++ b/docs/src/refs.bib
@@ -25,6 +25,17 @@ @inproceedings{cao-cpr
   organization={Society of Petroleum Engineers}
 }
 
+@article{moyner_tchelepi_2018,
+  title={A mass-conservative sequential implicit multiscale method for isothermal equation-of-state compositional problems},
+  author={M{\o}yner, Olav and Tchelepi, Hamdi A},
+  journal={SPE Journal},
+  volume={23},
+  number={06},
+  pages={2376--2393},
+  year={2018},
+  publisher={SPE}
+}
+
 @inproceedings{wallis-cpr,
   title={Incomplete Gaussian elimination as a preconditioning for generalized conjugate gradient acceleration},
   author={Wallis, John R},
@@ -33,3 +44,84 @@ @inproceedings{wallis-cpr
   year={1983},
   organization={Society of Petroleum Engineers}
 }
+
+@inproceedings{hypre,
+  title={hypre: A library of high performance preconditioners},
+  author={Falgout, Robert D and Yang, Ulrike Meier},
+  booktitle={International Conference on computational science},
+  pages={632--641},
+  year={2002},
+  organization={Springer}
+}
+
+@article{boomeramg,
+  title={BoomerAMG: A parallel algebraic multigrid solver and preconditioner},
+  author={Yang, Ulrike Meier and others},
+  journal={Applied Numerical Mathematics},
+  volume={41},
+  number={1},
+  pages={155--177},
+  year={2002},
+  publisher={Elsevier}
+}
+
+@article{amgcl,
+  title={AMGCL: An efficient, flexible, and extensible algebraic multigrid implementation},
+  author={Demidov, Denis},
+  journal={Lobachevskii Journal of Mathematics},
+  volume={40},
+  pages={535--546},
+  year={2019},
+  publisher={Springer}
+}
+
+@article{krylov_jl,
+  author  = {Montoison, Alexis and Orban, Dominique},
+  title   = {{Krylov.jl: A Julia basket of hand-picked Krylov methods}},
+  journal = {Journal of Open Source Software},
+  volume  = {8},
+  number  = {89},
+  pages   = {5187},
+  year    = {2023},
+  doi     = {10.21105/joss.05187}
+}
+
+@inproceedings{spe9,
+  title={Ninth {SPE} comparative solution project: a reexamination of black-oil simulation},
+  author={Killough, JE},
+  booktitle={SPE Reservoir Simulation Conference},
+  pages={SPE--29110},
+  year={1995},
+  organization={SPE}
+}
+@inproceedings{olympus,
+  title={Overview of the olympus field development optimization challenge},
+  author={Fonseca, RM and Della Rossa, E and Emerick, AA and Hanea, RG and Jansen, JD},
+  booktitle={ECMOR XVI-16th European Conference on the Mathematics of Oil Recovery},
+  volume={2018},
+  number={1},
+  pages={1--10},
+  year={2018},
+  organization={European Association of Geoscientists \& Engineers}
+}
+
+@article{adjoint_cpr,
+  title={Adaptation of the CPR preconditioner for efficient solution of the adjoint equation},
+  author={Han, Choongyong and Wallis, John and Sarma, Pallav and Li, Gary and Schrader, Mark L and Chen, Wen},
+  journal={{SPE} Journal},
+  volume={18},
+  number={02},
+  pages={207--213},
+  year={2013},
+  publisher={OnePetro}
+}
+
+@article{salo_co2,
+  title={Three-dimensional simulation of geologic carbon dioxide sequestration using MRST},
+  author={Sal{\'o}-Salgado, Llu{\'\i}s and M{\o}yner, Olav and Lie, Knut-Andreas and Juanes, Ruben},
+  journal={Advances in Geo-Energy Research},
+  volume={14},
+  number={1},
+  pages={34--48},
+  year={2024}
+}
diff --git a/examples/co2_props.jl b/examples/co2_props.jl
new file mode 100644
index 00000000..c4fc304c
--- /dev/null
+++ b/examples/co2_props.jl
@@ -0,0 +1,136 @@
+# # Plotting CO2 properties used in compositional model
+# JutulDarcy includes a set of pre-generated tables for simulation of CO2
+# storage in saline aquifers, as well as functions for calculating PVT and
+# solubility properties of CO2-brine mixtures with varying degrees of salinity.
+#
+# For more details on the property calculations used herein, please see the
+# paper [Three-dimensional simulation of geologic carbon dioxide sequestration
+# using MRST by L. Saló et al (2024)](https://doi.org/10.46690/ager.2024.10.06).
+#
+# This example demonstrates how to plot the properties, and does a basic
+# comparison of how the properties change when the aqueous phase has salts
+# added.
+
+using Jutul, JutulDarcy, GLMakie
+import JutulDarcy.CO2Properties: co2_brine_property_tables
+
+# ## Define plotting functions
+# We define a plotting function that varies temperature, and one that compares
+# the same property with and without salts for a given temperature.
+
+T = [20.0, 40.0, 60.0, 80.0, 100.0]
+p = range(5, 300, 100)
+
+function plot_property(prop, tab, component, title = "")
+    trans = x -> x
+    if prop == :viscosity
+        u = "(Pa s)"
+    elseif prop == :density
+        u = "(kg/m^3)"
+    elseif prop == :K
+        u = " (log10)"
+        trans = log10
+    else
+        u = ""
+    end
+    is_h2o = component == :H2O
+    fig = Figure(size = (800, 600))
+    ax = Axis(fig[1, 1], title = title, xlabel = "Pressure (bar)", ylabel = "$component $prop $u")
+
+    F = tab[prop]
+    if prop == :K
+        F = F.K
+    end
+    for (i, T_i) in enumerate(T)
+        val = F.(convert_to_si(p, :bar), T_i + 273.15)
+        if is_h2o
+            val = map(first, val)
+        else
+            val = map(last, val)
+        end
+        lines!(ax, p, trans.(val), label = "T=$(T_i) °C",
+            colormap = :hot,
+            color = i,
+            colorrange = (1, length(T)+3),
+            linewidth = 2,
+        )
+    end
+    axislegend(position = :lt, orientation = :horizontal)
+    fig
+end
+
+function plot_brine_comparison(prop, T, tab1, tab2, component, title = "")
+    trans = x -> x
+    if prop == :viscosity
+        u = "(Pa s)"
+    elseif prop == :density
+        u = "(kg/m^3)"
+    elseif prop == :K
+        u = ""
+        u = " (log10)"
+        trans = log10
+    else
+        u = ""
+    end
+    is_h2o = component == :H2O
+    fig = Figure(size = (800, 600))
+    ax = Axis(fig[1, 1], title = title, xlabel = "Pressure (bar)", ylabel = "$component $prop $u")
+
+    pbar = convert_to_si.(p, :bar)
+
+    F1 = tab1[prop]
+    F2 = tab2[prop]
+    if prop == :K
+        F1 = F1.K
+        F2 = F2.K
+    end
+    val1 = F1.(pbar, T + 273.15)
+    val2 = F2.(pbar, T + 273.15)
+
+    if is_h2o
+        val1 = map(first, val1)
+        val2 = map(first, val2)
+    else
+        val1 = map(last, val1)
+        val2 = map(last, val2)
+    end
+    lines!(ax, trans.(val1), label = "Water",
+        linewidth = 2,
+    )
+    lines!(ax, trans.(val2), label = "Brine",
+        linewidth = 2,
+    )
+    axislegend(position = :lt, orientation = :horizontal)
+    fig
+end
+# ## Generate tables
+# We get tables for a wide pressure and temperature range. The first set of
+# tables assumes no salts, and the second set of tables uses specified molar
+# fractions of salts in the aqueous phase.
+tab_water = co2_brine_property_tables()
+tab_brine = co2_brine_property_tables(
+    salt_names = ["NaCl", "KCl"],
+    salt_mole_fractions = [0.05, 0.01]
+)
+
+# ## Plot aqueous mass density
+plot_property(:density, tab_water, :H2O, "Pure phase density, no salts")
+# ## Plot gaseous mass density
+plot_property(:density, tab_water, :CO2, "Pure phase density, no salts")
+# ## Plot aqueous viscosity
+plot_property(:viscosity, tab_water, :H2O, "Pure phase viscosity, no salts")
+# ## Plot gaseous viscosity
+plot_property(:viscosity, tab_water, :CO2, "Pure phase viscosity, no salts")
+# ## Plot K-value of CO2
+# The K value defines the ratio between liquid and vapor phases at equilibrium
+# conditions and is closely related to solubility. For a component we relate the
+# liquid mass fraction ``x_i`` to the vapor mole fraction ``y_i`` of that
+# component by the K-value ``K_i``:
+#
+# ``y_i = K_i(p, T) x_i``
+#
+plot_property(:K, tab_water, :CO2, "K value of CO2 component in aqueous phase, no salts")
+# ## Compare viscosity with and without salts
+plot_brine_comparison(:density, 30.0, tab_water, tab_brine, :H2O, "Pure phase viscosity with salts")
+# ## Compare K value with and without salts
+plot_brine_comparison(:K, 30.0, tab_water, tab_brine, :CO2, "K value of CO2 component in aqueous phase")
diff --git a/examples/intro_example.jl b/examples/intro_example.jl
deleted file mode 100644
index 0aefc9e3..00000000
--- a/examples/intro_example.jl
+++ /dev/null
@@ -1,79 +0,0 @@
-
-
-using JutulDarcy, Jutul
-Darcy, bar, kg, meter, day = si_units(:darcy, :bar, :kilogram, :meter, :day)
-nx = ny = 25
-nz = 10
-cart_dims = (nx, ny, nz)
-physical_dims = (1000.0, 1000.0, 100.0).*meter
-g = CartesianMesh(cart_dims, physical_dims)
-domain = reservoir_domain(g, permeability = 0.3Darcy, porosity = 0.2)
-Injector = setup_vertical_well(domain, 1, 1, name = :Injector);
-Producer = setup_well(domain, (nx, ny, 1), name = :Producer);
-
-
-phases = (LiquidPhase(), VaporPhase())
-rhoLS = 1000.0kg/meter^3
-rhoGS = 100.0kg/meter^3
-reference_densities = [rhoLS, rhoGS]
-sys = ImmiscibleSystem(phases, reference_densities = reference_densities)
-
-model, parameters = setup_reservoir_model(domain, sys, wells = [Injector, Producer])
-
-
-c = [1e-6, 1e-4]/bar
-density = ConstantCompressibilityDensities(
-    p_ref = 100*bar,
-    density_ref = reference_densities,
-    compressibility = c
-)
-kr = BrooksCoreyRelativePermeabilities(sys, [2.0, 3.0])
-replace_variables!(model, PhaseMassDensities = density, RelativePermeabilities = kr);
-
-state0 = setup_reservoir_state(model,
-    Pressure = 120bar,
-    Saturations = [1.0, 0.0]
-)
-# ### Set up schedule with driving forces
-
-nstep = 25
-dt = fill(365.0day, nstep)
-pv = pore_volume(model, parameters)
-inj_rate = 1.5*sum(pv)/sum(dt)
-rate_target = TotalRateTarget(inj_rate)
-I_ctrl = InjectorControl(rate_target, [0.0, 1.0], density = rhoGS)
-bhp_target = BottomHolePressureTarget(100bar)
-P_ctrl = ProducerControl(bhp_target)
-controls = Dict(:Injector => I_ctrl, :Producer => P_ctrl)
-forces = setup_reservoir_forces(model, control = controls)
-
-wd, states, t = simulate_reservoir(state0, model, dt, parameters = parameters, forces = forces)
-##
-wd(:Producer)
-##
-wd(:Injector, :bhp)
-# ### Plot the well rates
-using GLMakie
-
-grat = wd[:Producer, :grat]
-lrat = wd[:Producer, :lrat]
-bhp = wd[:Injector, :bhp]
-
-fig = Figure(size = (1200, 400))
-
-ax = Axis(fig[1, 1],
-    title = "Injector",
-    xlabel = "Time / days",
-    ylabel = "Bottom hole pressure / bar")
-lines!(ax, t/day, bhp./bar)
-
-ax = Axis(fig[1, 2],
-    title = "Producer",
-    xlabel = "Time / days",
-    ylabel = "Production rate / m³/day")
-lines!(ax, t/day, abs.(grat).*day)
-lines!(ax, t/day, abs.(lrat).*day)
-
-fig
-# ### Launch interactive plotting of reservoir values
-plot_reservoir(model, states, key = :Saturations, step = 3)
diff --git a/examples/validation_compositional.jl b/examples/validation_compositional.jl
index 2528392d..b3e6c55d 100644
--- a/examples/validation_compositional.jl
+++ b/examples/validation_compositional.jl
@@ -1,18 +1,35 @@
 # # Validation of equation-of-state compositiona simulator
-# This example solves a 1D problem and compares against existing simulators
-# (E300, AD-GPRS) to verify correctness.
+# This example solves a 1D two-phase, three component miscible displacement
+# problem and compares against existing simulators (E300, AD-GPRS) to verify
+# correctness.
 #
 # The case is loaded from an input file that can be run in other simulators. For
-# convenience, we provide solutions from the other simulators as a binary file.
+# convenience, we provide solutions from the other simulators as a binary file
+# to perform a comparison without having to run and convert results from other
+# the simulators.
+#
+# This case is a small compositional problem inspired by the examples in Voskov
+# et al (JPSE, 2012). A 1D reservoir of 1,000 meters length is discretized into
+# 1,000 cells. The model initially contains a mixture made up of 0.6 parts C10,
+# 0.1 parts CO2, and 0.3 parts C1 by moles at 150 degrees C and 75 bar pressure.
+# Wells are placed in the leftmost and rightmost cells of the domain, with the
+# leftmost well injecting pure CO$_2$ at a fixed bottom-hole pressure of 100 bar
+# and the other well producing at 50 bar. The model is isothermal and contains a
+# phase transition from the initial two-phase mixture to single-phase gas as
+# injected CO$_2$ eventually displaces the resident fluids. For further details
+# on this setup, see Møyner and Tchelepi (SPE J. 2018)
+# [moyner_tchelepi_2018](@cite).
 using JutulDarcy
 using Jutul
 using GLMakie
 dpth = JutulDarcy.GeoEnergyIO.test_input_file_path("SIMPLE_COMP")
 data_path = joinpath(dpth, "SIMPLE_COMP.DATA")
 case = setup_case_from_data_file(data_path)
-result = simulate_reservoir(case)
+result = simulate_reservoir(case, info_level = 1)
 ws, states = result;
 # ## Plot solutions and compare
+# The 1D displacement can be plotted as a line plot. We pick a step midway
+# through the simulation and plot compositions, saturations and pressure.
 cmap = :tableau_hue_circle
 ref_path = joinpath(dpth, "reference.jld2")
 ref = Jutul.JLD2.load(ref_path)
diff --git a/examples/validation_egg.jl b/examples/validation_egg.jl
index 9a63a40c..9be32053 100644
--- a/examples/validation_egg.jl
+++ b/examples/validation_egg.jl
@@ -92,8 +92,11 @@ function plot_well_comparison(response, well_names, reponse_name = "$response")
     fig
 end
 # ## Well responses and comparison
+# As the case is a two-phase model with water injection, we limit the results to
+# plots of the producer water and oil rates.
+
 # ### Water production rates
-plot_well_comparison(:wrat, prod, "Water surface rate")
-# ### Oil production rates
-plot_well_comparison(:orat, prod, "Oil surface rate")
+plot_well_comparison(:wrat, prod, "Producer water surface rate")
 
+# ### Oil production rates
+plot_well_comparison(:orat, prod, "Producer oil surface rate")
diff --git a/examples/validation_olympus_1.jl b/examples/validation_olympus_1.jl
index dea2bfa4..1af60655 100644
--- a/examples/validation_olympus_1.jl
+++ b/examples/validation_olympus_1.jl
@@ -1,6 +1,13 @@
 # # OLYMPUS_1 model
-# Model from the [ISAPP Optimization challenge](https://www.isapp2.com/optimization-challenge/reservoir-model-description.html)
+# Model from the [ISAPP Optimization
+# challenge](https://www.isapp2.com/optimization-challenge/reservoir-model-description.html)
+#
+# Two-phase complex corner-point model with primary and secondary production.
+#
+# For more details, see [olympus](@cite)
+#
 using Jutul, JutulDarcy, GLMakie, DelimitedFiles, HYPRE
+using Test #src
 olympus_dir = JutulDarcy.GeoEnergyIO.test_input_file_path("OLYMPUS_1")
 case = setup_case_from_data_file(joinpath(olympus_dir, "OLYMPUS_1.DATA"), backend = :csr)
 ws, states = simulate_reservoir(case, output_substates = true)
@@ -55,6 +62,11 @@ function plot_well_comparison(response, well_names, reponse_name = "$response";
         ref_pos = findfirst(x -> x == label_in_csv, vec(header))
         qoi = copy(well[response]).*ys
         qoi_ref = data[:, ref_pos].*ys
+
+        val = sum(qoi.*diff([0; time_jutul]))       #src
+        val_ref = sum(qoi_ref.*diff([0; time_ref])) #src
+        @test abs(val - val_ref)/abs(val) < 0.03    #src
+
         tot_rate = copy(well[:rate])
         @. qoi[tot_rate == 0] = NaN
         orat_ref = data[:, findfirst(x -> x == "$well_name:orat", vec(header))]
diff --git a/examples/validation_spe9.jl b/examples/validation_spe9.jl
index a45dd9fe..4400512e 100644
--- a/examples/validation_spe9.jl
+++ b/examples/validation_spe9.jl
@@ -3,6 +3,8 @@
 # of black-oil simulation. In SPE Reservoir Simulation Symposium,  12-15
 # February 1995, San Antonio, Texas. SPE 29110-MS]
 # (http://dx.doi.org/10.2118/29110-MS)
+# [spe9](@cite)
+
 using Jutul, JutulDarcy, GLMakie, DelimitedFiles
 spe9_dir = JutulDarcy.GeoEnergyIO.test_input_file_path("SPE9")
 case = setup_case_from_data_file(joinpath(spe9_dir, "SPE9.DATA"))
diff --git a/src/CO2Properties/generation.jl b/src/CO2Properties/generation.jl
index dd2f216a..0fc1a6c7 100644
--- a/src/CO2Properties/generation.jl
+++ b/src/CO2Properties/generation.jl
@@ -278,16 +278,18 @@ of excess Gibbs energy.
     - rhox: Scalar with density in [mol/m^3]
     - rho: Scalar with density in [kg/m^3]
 """
-function activity_co2_DS2003(T, P, m_io)
-    # Check if T, P conditions are within range
-    if P > 2000
-        jutul_message("activity_co2_DS2003", "Pressure out of tested range", color = :yellow)
-    end
-    if T < 273 || T > 533
-        jutul_message("activity_co2_DS2003", "Temperature out of tested range", color = :yellow)
-    end
-    if maximum(m_io) > 4.3
-        jutul_message("activity_co2_DS2003", "Ionic strength out of tested range", color = :yellow)
+function activity_co2_DS2003(T, P, m_io; check = true)
+    if check
+        # Check if T, P conditions are within range
+        if P > 2000
+            jutul_message("activity_co2_DS2003", "Pressure out of tested range", color = :yellow)
+        end
+        if T < 273 || T > 533
+            jutul_message("activity_co2_DS2003", "Temperature out of tested range", color = :yellow)
+        end
+        if maximum(m_io) > 4.3
+            jutul_message("activity_co2_DS2003", "Ionic strength out of tested range", color = :yellow)
+        end
     end
 
     # Interaction parameter constants
@@ -543,9 +545,10 @@ calculation) and can any subset of the following names provided:
 
 """
 function compute_co2_brine_props(p_pascal, T_K, salt_mole_fractions = Float64[], salt_names = String[];
-        check=true,
-        iterate=false,
-        maxits=15,
+        check = true,
+        iterate = false,
+        verbose = false,
+        maxits = 15,
         ionized = false
     )
     # TODO: Salt var name
@@ -709,7 +712,7 @@ function compute_co2_brine_props(p_pascal, T_K, salt_mole_fractions = Float64[],
         # 5. CO2 molality in saline solution at P, T conditions (m_co2)
         #   5.1 CO2 activity coefficient (Duan and Sun, 2003, as in Spycher
         #   & Pruess, 2005; typos in Hassanzadeh et al., 2008)
-        gamma_co2 = activity_co2_DS2003(T_K, P, m_io)
+        gamma_co2 = activity_co2_DS2003(T_K, P, m_io, check = check)
 
         #   5.2 CO2 molality
         m_co2 = m0_co2 / gamma_co2
@@ -740,7 +743,7 @@ function compute_co2_brine_props(p_pascal, T_K, salt_mole_fractions = Float64[],
             Y0_h2o = Y_h2o
             Y0_co2 = Y_co2
             if res < 1e-6
-                println("Iterations at p=$P bar: $it")
+                verbose && println("Iterations at p=$P bar: $it")
                 done = true
             elseif it == maxits
                 error("Iterative loop did not converge")
diff --git a/src/CO2Properties/setup.jl b/src/CO2Properties/setup.jl
index 5de79112..258003d8 100644
--- a/src/CO2Properties/setup.jl
+++ b/src/CO2Properties/setup.jl
@@ -3,26 +3,20 @@ function setup_reservoir_model_co2_brine(reservoir::DataDomain;
         thermal = false,
         co2_physics = :kvalue,
         co2_table_directory = missing,
+        co2_source = missing,
+        co2_density = :nist,
         extra_out = true,
         salt_names = String[],
         salt_mole_fractions = Float64[],
-        co2_source = ifelse(ismissing(co2_table_directory), :salo24, :table),
         kwarg...
     )
-    length(salt_mole_fractions) == length(salt_names) || throw(ArgumentError("salt_names ($salt_names) and salt_mole_fractions ($salt_mole_fractions) must have equal length."))
-    tables = co2_brine_property_tables(temperature, basepath = co2_table_directory)
-    if co2_source == :salo24 || length(salt_names) > 0
-        dens = tables[:density]::Jutul.BilinearInterpolant
-        visc = tables[:viscosity]::Jutul.BilinearInterpolant
-        K = tables[:K].K::Jutul.BilinearInterpolant
-        replace_co2_brine_properties!(dens, visc, K, salt_mole_fractions, salt_names)
-    elseif co2_source == :table
-        if length(salt_mole_fractions) > 0
-            jutul_message("Salts $salt_names were provided but table was also provided as $co2_table_directory. Salts will not be accounted for.", color = :yellow)
-        end
-    else
-        co2_source == :csp11 || throw(ArgumentError("co2_source argument must be either :csp11, :table or :salo24"))
-    end
+    tables = co2_brine_property_tables(temperature,
+        basepath = co2_table_directory,
+        salt_names = salt_names,
+        salt_mole_fractions = salt_mole_fractions,
+        co2_source = co2_source,
+        co2_density = co2_density
+    )
     rho = JutulDarcy.BrineCO2MixingDensities(tables[:density])
     mu = JutulDarcy.PTViscosities(tables[:viscosity])
     if thermal
@@ -90,10 +84,15 @@ function setup_reservoir_model_co2_brine(reservoir::DataDomain;
     return out
 end
 
-function replace_co2_brine_properties!(dens, visc, K, salt_mole_fractions, salt_names)
+function replace_co2_brine_properties!(dens, visc, K, salt_mole_fractions, salt_names; co2_density::Symbol, T_value)
     p = dens.X
-    T = dens.Y
-
+    if ismissing(T_value)
+        T = dens.Y
+    else
+        T = [T_value]
+    end
+    co2_density in (:rk, :nist) || throw(ArgumentError("co2_density argument must be :rk (for Redlich-Kwong) or :nist (for NIST tables)"))
+    replace_co2 = co2_density != :nist
     for (i, p_i) in enumerate(p)
         if i == 1
             p_i = p[2]
@@ -101,12 +100,24 @@ function replace_co2_brine_properties!(dens, visc, K, salt_mole_fractions, salt_
             p_i = p[end-1]
         end
         for (j, T_i) in enumerate(T)
-            if j == 1
-                T_i = T[2]
-            elseif j == length(T)
-                T_i = T[end-1]
+            if length(T) == 1
+                T_i = only(T)
+            else
+                if j == 1
+                    T_i = T[2]
+                elseif j == length(T)
+                    T_i = T[end-1]
+                end
             end
             props = compute_co2_brine_props(p_i, T_i, salt_mole_fractions, salt_names, check = false)
+            dens_calc = props[:density]
+            if replace_co2
+                dens.F[i, j] = dens_calc
+            else
+                pair_type = eltype(dens.F)
+                v_co2 = dens.F[i, j][2]
+                dens.F[i, j] = pair_type(dens_calc[1], v_co2)
+            end
             dens.F[i, j] = props[:density]
             visc.F[i, j] = props[:viscosity]
             K.F[i, j] = props[:K]
@@ -124,11 +135,47 @@ function JutulDarcy.get_phases(::Val{:co2brine})
     return (LiquidPhase(), VaporPhase())
 end
 
+"""
+    co2_brine_property_tables(T = missing;
+        basepath = missing,
+        salt_names = salt_names,
+        salt_mole_fractions = salt_mole_fractions,
+        co2_source = missing,
+        co2_density = :rk
+    )
+
+Set up brine-CO2 property tables for a wide range of p, T. If a single
+temperature value `T` is provided, an isothermal table will be generated for
+that temperature.
+
 
-function co2_brine_property_tables(T = missing; basepath = missing)
+Keyword arguments:
+
+    - `basepath`: Path to folder containing tables. This can be used to override the default tables.
+    - `salt_names`: Names of salts to in include (passed onto `compute_co2_brine_props`)
+    - `salt_mole_fractions`: Corresponding salt mole fractions (same length as `salt_names`)
+    - `co2_source = missing`: Source of CO2 values. Can be:
+        - `:table`: Use table directly.
+        - `:salo24`: Use correlations implemented in paper by Salo et al (support for salts)
+        - `:csp11`: Use values exactly as given in CSP11 benchmark (no salts)
+    - `co2_density=:rk`: How to obtain density of pure CO2 phase. Can be `:nist` for NIST tables or `:rk` for adjusted Redlich-Kwong from Spycher paper.
+
+"""
+function co2_brine_property_tables(T = missing;
+        basepath = missing,
+        salt_names = String[],
+        salt_mole_fractions = Float64[],
+        co2_source = missing,
+        co2_density = :rk
+    )
+    if ismissing(co2_source)
+        co2_source = ifelse(ismissing(basepath), :salo24, :table)
+    end
     if ismissing(basepath)
         basepath = joinpath(artifact"CO2Tables_CSP11", "csp11")
     end
+    length(salt_mole_fractions) == length(salt_names) || throw(ArgumentError("salt_names ($salt_names) and salt_mole_fractions ($salt_mole_fractions) must have equal length."))
+
     ispath(basepath) || throw(ArgumentError("basepath $basepath does not exist."))
     getpth(n) = joinpath(basepath, n)
 
@@ -138,16 +185,26 @@ function co2_brine_property_tables(T = missing; basepath = missing)
 
     prop_tab = [h2o, co2]
     K_pT = co2brine_K_values(sol, T)
-    rho_pT = co2brine_phase_property_table(prop_tab, :density, T)
-    mu_pT = co2brine_phase_property_table(prop_tab, :viscosity, T)
-    H_pT = co2brine_phase_property_table(prop_tab, :H, T)
-    C_v = co2brine_phase_property_table(prop_tab, :cv, T)
-    C_p = co2brine_phase_property_table(prop_tab, :cp, T)
-    E = co2brine_phase_property_table(prop_tab, :E, T)
+    get_tab(k::Symbol) = co2brine_phase_property_table(prop_tab, k, T)
+    rho_pT = get_tab(:density)
+    mu_pT = get_tab(:viscosity)
+    H_pT = get_tab(:H)
+    C_v = get_tab(:cv)
+    C_p = get_tab(:cp)
+    E = get_tab(:E)
+    phase_conductivity = get_tab(:phase_conductivity)
 
-    phase_conductivity = co2brine_phase_property_table(prop_tab, :phase_conductivity, T)
+    if co2_source == :table
+        if length(salt_mole_fractions) > 0
+            jutul_message("Salts $salt_names were provided but table was also provided as $base_dir. Salts will not be accounted for.", color = :yellow)
+        end
+    elseif co2_source == :salo24 || length(salt_names) > 0
+        replace_co2_brine_properties!(rho_pT, mu_pT, K_pT.K, salt_mole_fractions, salt_names, co2_density = co2_density, T_value = T)
+    else
+        co2_source == :csp11 || throw(ArgumentError("co2_source argument must be either :csp11, :table or :salo24"))
+    end
 
-    return Dict(
+    tables = Dict(
         :K => K_pT,
         :density => rho_pT,
         :viscosity => mu_pT,
@@ -160,6 +217,7 @@ function co2_brine_property_tables(T = missing; basepath = missing)
         :phase_conductivity => phase_conductivity,
         :solubility_table => sol
     )
+    return tables
 end
 
 function JutulDarcy.select_injector_mixture_spec(sys::Val{:co2brine}, name, streams, type)
diff --git a/src/facility/types.jl b/src/facility/types.jl
index c509a779..94c1f1af 100644
--- a/src/facility/types.jl
+++ b/src/facility/types.jl
@@ -59,7 +59,15 @@ cells.
 """
 struct Perforations <: JutulEntity end
 
+"""
+    WellIndices()
+
+Parameter for the connection strength between a well and the reservoir for a
+given perforation. Typical values come from a combination of Peaceman's formula,
+upscaling and/or history matching.
+"""
 struct WellIndices <: ScalarVariable end
+
 Jutul.minimum_value(::WellIndices) = 0.0
 Jutul.variable_scale(::WellIndices) = 1e-10
 
@@ -69,6 +77,12 @@ function Jutul.default_values(model, ::WellIndices)
     return vec(copy(w.perforations.WI))
 end
 
+"""
+    PerforationGravityDifference()
+
+Parameter for the height difference from the wellbore and the connected node in
+the well.
+"""
 struct PerforationGravityDifference <: ScalarVariable end
 
 Jutul.associated_entity(::PerforationGravityDifference) = Perforations()
diff --git a/src/gradients/gradients.jl b/src/gradients/gradients.jl
index a61c9d42..313205de 100644
--- a/src/gradients/gradients.jl
+++ b/src/gradients/gradients.jl
@@ -33,15 +33,26 @@ function reservoir_sensitivities(case::JutulCase, rsr::ReservoirSimResult, objec
     return reservoir_sensitivities(case, rsr.result, objective; kwarg...)
 end
 
-function reservoir_sensitivities(case::JutulCase, result::Jutul.SimResult, obj; include_parameters = false, adjoint_arg = NamedTuple(), kwarg...)
+function reservoir_sensitivities(case::JutulCase, result::Jutul.SimResult, obj;
+        include_parameters = false,
+        include_state0 = false,
+        adjoint_arg = (include_state0 = include_state0, ),
+        kwarg...)
     sens = solve_adjoint_sensitivities(case, result, obj; adjoint_arg...)
     rmodel = reservoir_model(case.model)
     if haskey(sens, :Reservoir)
         sens = sens[:Reservoir]
     end
     data_domain_with_gradients = Jutul.data_domain_to_parameters_gradient(rmodel, sens; kwarg...)
+    vars_to_add = []
     if include_parameters
-        for (k, pdef) in pairs(Jutul.get_parameters(rmodel))
+        push!(vars_to_add, Jutul.get_parameters(rmodel))
+    end
+    if include_state0
+        push!(vars_to_add, Jutul.get_primary_variables(rmodel))
+    end
+    for vars in vars_to_add
+        for (k, pdef) in pairs(vars)
             u = associated_entity(pdef)
             data_domain_with_gradients[k, u] = sens[k]
         end
diff --git a/src/input_simulation/data_input.jl b/src/input_simulation/data_input.jl
index 30ea2ce0..a7371cdc 100644
--- a/src/input_simulation/data_input.jl
+++ b/src/input_simulation/data_input.jl
@@ -69,7 +69,7 @@ function setup_case_from_parsed_data(datafile;
         normalize = true,
         convert_co2store = true,
         repair_zcorn = true,
-        process_pinch = false,
+        process_pinch = true,
         kwarg...
     )
     function msg(s)
@@ -556,7 +556,7 @@ end
 function parse_forces(model, datafile, sys, wells, controls, limits, cstep, dt, well_forces)
     bc, sources = parse_aquifer_bc(model, datafile, sys)
     if length(wells) == 0
-        return setup_reservoir_forces(mode, bc = bc, sources = sources)
+        return setup_reservoir_forces(model, bc = bc, sources = sources)
     end
     forces = []
     @assert length(controls) == length(limits) == length(well_forces)
@@ -804,10 +804,10 @@ function initialize_numerical_aquifers!(init, rmodel, aquifers)
     return init
 end
 
-function parse_reservoir(data_file; zcorn_depths = true, repair_zcorn = true, process_pinch = false)
+function parse_reservoir(data_file; zcorn_depths = true, repair_zcorn = true, process_pinch = true)
     grid = data_file["GRID"]
     cartdims = grid["cartDims"]
-    G = mesh_from_grid_section(grid, missing, repair_zcorn, process_pinch)
+    G = mesh_from_grid_section(grid; repair_zcorn = repair_zcorn, process_pinch = process_pinch)
 
     # Handle numerical aquifers
     aqunum = get(grid, "AQUNUM", missing)
@@ -1057,6 +1057,7 @@ function parse_reservoir(data_file; zcorn_depths = true, repair_zcorn = true, pr
         @. domain[:cell_centroids][3, :] = z
     end
 
+    # pinch_face = get_mesh_entity_tag(G, Faces(), :cpgrid_connection_type, :pinched, throw = false)
     if !isnothing(aquifers)
         domain[:numerical_aquifers, nothing] = aquifers
         vol = domain[:volumes]
@@ -2057,7 +2058,7 @@ function apply_weltarg!(controls, limits, wk)
     well, ctype, value = wk
     ctype = lowercase(ctype)
     ctrl = controls[well]
-    @assert !(ctrl isa DisabledControl) "Cannot use WELTARG on disabled well."
+    @assert !(ctrl isa DisabledControl) "Cannot use WELTARG $ctype = $value on disabled well $well."
     is_injector = ctrl isa InjectorControl
     limit = limits[well]
     limit = OrderedDict{Symbol, Float64}(pairs(limit))
diff --git a/src/multiphase.jl b/src/multiphase.jl
index 3dced3e4..a173b896 100644
--- a/src/multiphase.jl
+++ b/src/multiphase.jl
@@ -129,6 +129,12 @@ end
 
 Jutul.parameter_is_differentiable(::Saturations, model) = false
 
+"""
+    ConnateWater()
+
+Parameter for connate water per cell. Used in some three-phase relative
+permeability evaluators.
+"""
 struct ConnateWater <: ScalarVariable end
 
 function Jutul.default_values(model, ::ConnateWater)
diff --git a/src/thermal/thermal.jl b/src/thermal/thermal.jl
index df6bcec3..34a1dfb8 100644
--- a/src/thermal/thermal.jl
+++ b/src/thermal/thermal.jl
@@ -225,6 +225,14 @@ function Jutul.default_parameter_values(data_domain, model, param::RockThermalCo
     return T
 end
 
+"""
+    add_thermal_to_model!(model::MultiModel)
+
+Add energy conservation equation and thermal primary variable together with
+standard set of parameters to existing flow model. Note that more complex models
+require additional customization after this function call to get correct
+results.
+"""
 function add_thermal_to_model!(model::MultiModel)
     for (k, m) in pairs(model.models)
         if m.system isa MultiPhaseSystem
diff --git a/src/types.jl b/src/types.jl
index d5de5c50..3866c277 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -368,6 +368,9 @@ end
 abstract type PorousMediumDomain <: JutulMesh end
 abstract type ReservoirGrid <: PorousMediumDomain end
 
+"""
+Abstract supertype for all well domains.
+"""
 abstract type WellDomain <: PorousMediumDomain
     # Wells are not porous themselves per se, but they are discretizing
     # part of a porous medium.
diff --git a/src/variables/relperm/endscale.jl b/src/variables/relperm/endscale.jl
index 0ae463d0..acb18816 100644
--- a/src/variables/relperm/endscale.jl
+++ b/src/variables/relperm/endscale.jl
@@ -79,6 +79,13 @@ struct EndPointScalingCoefficients{phases} <: VectorVariables
     drainage::Bool
 end
 
+"""
+    EndPointScalingCoefficients(phases::Symbol; drainage = true)
+
+Parameters that represent end-point scaling for a relative permeability phase
+(drainage or imbibition) as parametrized by four points (connate, critical,
+saturation at which maximum kr first occurs and maximum kr).
+"""
 function EndPointScalingCoefficients(phases::Symbol; drainage = true)
     return EndPointScalingCoefficients{phases}(drainage)
 end