diff --git a/docs/source/api/constants_api.md b/docs/source/api/constants_api.md
index 1b992418..0565da39 100644
--- a/docs/source/api/constants_api.md
+++ b/docs/source/api/constants_api.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # The {mod}`~pyrealm.constants` module
diff --git a/docs/source/api/core_api.md b/docs/source/api/core_api.md
index 1324720e..528f10d3 100644
--- a/docs/source/api/core_api.md
+++ b/docs/source/api/core_api.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # The {mod}`~pyrealm.core` module
diff --git a/docs/source/api/demography_api.md b/docs/source/api/demography_api.md
index 9290a0b0..db9dd62e 100644
--- a/docs/source/api/demography_api.md
+++ b/docs/source/api/demography_api.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # The {mod}`~pyrealm.demography` module
diff --git a/docs/source/api/pmodel_api.md b/docs/source/api/pmodel_api.md
index 9bb4f669..b81a90d0 100644
--- a/docs/source/api/pmodel_api.md
+++ b/docs/source/api/pmodel_api.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # The {mod}`~pyrealm.pmodel` module
diff --git a/docs/source/api/splash_api.md b/docs/source/api/splash_api.md
index fe55c923..3cb0f3ac 100644
--- a/docs/source/api/splash_api.md
+++ b/docs/source/api/splash_api.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # The {mod}`~pyrealm.splash` module
diff --git a/docs/source/api/tmodel_api.md b/docs/source/api/tmodel_api.md
index b7fa9312..de3c0158 100644
--- a/docs/source/api/tmodel_api.md
+++ b/docs/source/api/tmodel_api.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # The {mod}`~pyrealm.tmodel` module
diff --git a/docs/source/development/code_qa_and_typing.md b/docs/source/development/code_qa_and_typing.md
index 848f7a5c..ed2c797d 100644
--- a/docs/source/development/code_qa_and_typing.md
+++ b/docs/source/development/code_qa_and_typing.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # Code quality and static typing
diff --git a/docs/source/development/code_testing.md b/docs/source/development/code_testing.md
index 43829394..3d76d292 100644
--- a/docs/source/development/code_testing.md
+++ b/docs/source/development/code_testing.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # Package testing and profiling
diff --git a/docs/source/development/documentation.md b/docs/source/development/documentation.md
index 0ce6b2f8..5552f055 100644
--- a/docs/source/development/documentation.md
+++ b/docs/source/development/documentation.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # Documentation
@@ -60,6 +71,16 @@ kernelspec:
   display_name: Python 3
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 ```
 
diff --git a/docs/source/development/github_actions.md b/docs/source/development/github_actions.md
index d26a7452..a0815478 100644
--- a/docs/source/development/github_actions.md
+++ b/docs/source/development/github_actions.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # GitHub Actions
diff --git a/docs/source/development/overview.md b/docs/source/development/overview.md
index 837b1efe..8138127f 100644
--- a/docs/source/development/overview.md
+++ b/docs/source/development/overview.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # Developing `pyrealm`
diff --git a/docs/source/development/profiling_and_benchmarking.md b/docs/source/development/profiling_and_benchmarking.md
index d705d43f..4570cf26 100644
--- a/docs/source/development/profiling_and_benchmarking.md
+++ b/docs/source/development/profiling_and_benchmarking.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # Profiling and benchmarking
diff --git a/docs/source/development/pyrealm_build_data.md b/docs/source/development/pyrealm_build_data.md
index b59ec2b2..b94a4d87 100644
--- a/docs/source/development/pyrealm_build_data.md
+++ b/docs/source/development/pyrealm_build_data.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # The `pyrealm_build_data` package
diff --git a/docs/source/development/release_process.md b/docs/source/development/release_process.md
index e9ff6b00..27dd9021 100644
--- a/docs/source/development/release_process.md
+++ b/docs/source/development/release_process.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # Package release process
diff --git a/docs/source/index.md b/docs/source/index.md
index 43aef59b..32ec1e69 100644
--- a/docs/source/index.md
+++ b/docs/source/index.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # The `pyrealm` package
diff --git a/docs/source/users/constants.md b/docs/source/users/constants.md
index 899076b3..91600360 100644
--- a/docs/source/users/constants.md
+++ b/docs/source/users/constants.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # Package constants
@@ -35,7 +46,7 @@ of default model constants. The core API and details for each class can be seen
 
 These can be used to generate the default set of model parameters:
 
-```{code-cell}
+```{code-cell} ipython3
 from pyrealm.constants import CoreConst, TModelTraits
 
 core_const = CoreConst()
@@ -47,7 +58,7 @@ print(tmodel_const)
 
 And individual values can be altered using the parameter arguments:
 
-```{code-cell}
+```{code-cell} ipython3
 # Estimate processes under the moon's gravity...
 core_const_moon = CoreConst(k_G=1.62)
 # ... allowing a much greater maximum height
@@ -61,7 +72,7 @@ In order to ensure that a set of parameters cannot change while models are being
 instances of these parameter classes are **frozen**. You cannot  edit an existing
 instance and will need to create a new instance to use different parameters.
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [raises-exception]
 
 core_const_moon.k_G = 9.80665
@@ -75,7 +86,7 @@ export of parameter settings to dictionaries and to JSON formatted files. The co
 shows these methods working. First, a trait definition in a JSON file is read into a
 dictionary:
 
-```{code-cell}
+```{code-cell} ipython3
 import json
 import pprint
 
@@ -88,7 +99,7 @@ the {meth}`~pyrealm.constants.base.ConstantsClass.from_dict` method. The
 {meth}`~pyrealm.constants.base.ConstantsClass.from_json` method allows this to
 be done more directly and the resulting instances are identical.
 
-```{code-cell}
+```{code-cell} ipython3
 traits1 = TModelTraits.from_dict(trt_dict)
 traits2 = TModelTraits.from_json("../files/traits.json")
 
diff --git a/docs/source/users/demography/canopy.md b/docs/source/users/demography/canopy.md
index f9c0688c..2ec8508d 100644
--- a/docs/source/users/demography/canopy.md
+++ b/docs/source/users/demography/canopy.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # Canopy model
@@ -30,7 +41,7 @@ The key variables in calculating the canopy model are the crown projected area $
 and leaf projected area $\tilde{A}_{cp}(z)$, which are calculated for a stem
 of a given size using the  [crown model](./crown.md).
 
-```{code-cell}
+```{code-cell} ipython3
 from matplotlib import pyplot as plt
 import numpy as np
 import pandas as pd
@@ -85,7 +96,7 @@ and the required canopy gap fraction.
 
 The code below creates a simple community:
 
-```{code-cell}
+```{code-cell} ipython3
 # Two PFTs
 # - a shorter understory tree with a columnar canopy and no crown gaps
 # - a taller canopy tree with a top heavy canopy and more crown gaps
@@ -115,14 +126,14 @@ community = Community(
 
 We can then look at the expected allometries for the stems in each cohort:
 
-```{code-cell}
+```{code-cell} ipython3
 print("H = ", community.stem_allometry.stem_height)
 print("Ac = ", community.stem_allometry.crown_area)
 ```
 
 We can now calculate the canopy model for the community:
 
-```{code-cell}
+```{code-cell} ipython3
 canopy = Canopy(community=community, canopy_gap_fraction=2 / 32)
 ```
 
@@ -133,14 +144,14 @@ heights for each stem in the three cohorts.
 There are four canopy layers, with the top two very close together because of the
 large crown area in the two stems in the cohort of `tall` trees.
 
-```{code-cell}
+```{code-cell} ipython3
 canopy.layer_heights
 ```
 
 The `stem_crown_area` attribute then provides the crown area of each stem found in each
 layer.
 
-```{code-cell}
+```{code-cell} ipython3
 canopy.stem_crown_area
 ```
 
@@ -149,7 +160,7 @@ the first two layers are taken up entirely by the two stems in the cohort of lar
 trees. We can confirm that the calculation is correct by calculating the total crown area
 across the cohorts at each height:
 
-```{code-cell}
+```{code-cell} ipython3
 np.sum(canopy.stem_crown_area * community.cohort_data["n_individuals"], axis=1)
 ```
 
@@ -163,7 +174,7 @@ identical to the projected crown area for the first two cohorts because the crow
 fraction $f_g$ is zero for this PFT. The projected leaf area is however displaced
 towards the ground in the last cohort, because the `tall` PFT has a large gap fraction.
 
-```{code-cell}
+```{code-cell} ipython3
 canopy.stem_leaf_area
 ```
 
@@ -174,7 +185,7 @@ community crown and leaf area profile across a range of height values. For each
 we calculate the sum of the product of stem projected area and the number of
 individuals in each cohort.
 
-```{code-cell}
+```{code-cell} ipython3
 # Set of vertical height to calculate crown profiles
 at_z = np.linspace(0, 26, num=261)[:, None]
 
@@ -198,7 +209,7 @@ superimpose the calculated $z^*_l$ values and the cumulative canopy area for eac
 to confirm that the calculated values coincide with the profile. Note here that the
 total area at each closed layer height is omitting the community gap fraction.
 
-```{code-cell}
+```{code-cell} ipython3
 fig, ax = plt.subplots(ncols=1)
 
 # Calculate the crown area at which each canopy layer closes.
@@ -256,7 +267,7 @@ $f_{abs} = 1 - e ^ {-kL}$,
 where $k$ is the light extinction coefficient ($k$) and $L$ is the leaf area index
 (LAI). The LAI can be calculated for each stem and layer:
 
-```{code-cell}
+```{code-cell} ipython3
 # LAI = Acp_within_layer / canopy_area
 # print(LAI)
 ```
@@ -264,7 +275,7 @@ where $k$ is the light extinction coefficient ($k$) and $L$ is the leaf area ind
 This can be used to calculate the LAI of individual stems but also the LAI of each layer
 in the canopy:
 
-```{code-cell}
+```{code-cell} ipython3
 # LAI_stem = LAI.sum(axis=0)
 # LAI_layer = LAI.sum(axis=1)
 
@@ -275,7 +286,7 @@ in the canopy:
 The layer LAI values can now be used to calculate the light transmission of each layer and
 hence the cumulative light extinction profile through the canopy.
 
-```{code-cell}
+```{code-cell} ipython3
 # f_abs = 1 - np.exp(-pft.traits.par_ext * LAI_layer)
 # ext = np.cumprod(f_abs)
 
@@ -287,7 +298,7 @@ One issue that needs to be resolved is that the T Model implementation in `pyrea
 follows the original implementation of the T Model in having LAI as a fixed trait of
 a given plant functional type, so is constant for all stems of that PFT.
 
-```{code-cell}
+```{code-cell} ipython3
 # print("f_abs = ", (1 - np.exp(-pft.traits.par_ext * pft.traits.lai)))
 ```
 
diff --git a/docs/source/users/demography/community.md b/docs/source/users/demography/community.md
index 0961b8e2..8bee4431 100644
--- a/docs/source/users/demography/community.md
+++ b/docs/source/users/demography/community.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # Plant Communities
@@ -19,7 +30,7 @@ This area of `pyrealm` is in active development.
 
 :::
 
-```{code-cell}
+```{code-cell} ipython3
 from matplotlib import pyplot as plt
 import numpy as np
 import pandas as pd
@@ -28,7 +39,7 @@ from pyrealm.demography.flora import PlantFunctionalType, Flora
 from pyrealm.demography.community import Community
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 short_pft = PlantFunctionalType(
     name="short", h_max=15, m=1.5, n=1.5, f_g=0, ca_ratio=380
 )
@@ -52,10 +63,10 @@ community = Community(
 )
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 community
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 
 ```
diff --git a/docs/source/users/demography/crown.md b/docs/source/users/demography/crown.md
index dcf41bc7..6f542632 100644
--- a/docs/source/users/demography/crown.md
+++ b/docs/source/users/demography/crown.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # The tree crown model
@@ -20,7 +31,7 @@ notes and initial demonstration code.
 
 :::
 
-```{code-cell}
+```{code-cell} ipython3
 from matplotlib import pyplot as plt
 import numpy as np
 import pandas as pd
@@ -146,7 +157,7 @@ The {class}`~pyrealm.demography.flora.PlantFunctionalType` class is typically
 used to set specific PFTs, but the functions to calculate $q_m$ and $p_{zm}$
 are used directly below to provides a demonstration of the impacts of each trait.
 
-```{code-cell}
+```{code-cell} ipython3
 # Set a range of values for m and n traits
 m = n = np.arange(1.0, 5, 0.1)
 
@@ -155,7 +166,7 @@ q_m = calculate_crown_q_m(m=m, n=n[:, None])
 z_max_prop = calculate_crown_z_max_proportion(m=m, n=n[:, None])
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(10.9, 4))
 
 # Plot q_m as a function of m and n
@@ -191,7 +202,7 @@ profiles for PFTs. It requires:
 The code below creates a set of PFTS with differing crown trait values and then creates
 a `Flora` object using the PFTs.
 
-```{code-cell}
+```{code-cell} ipython3
 # A PFT with a small crown area and equal m and n values
 narrow_pft = PlantFunctionalType(name="narrow", h_max=20, m=1.5, n=1.5, ca_ratio=20)
 # A PFT with an intermediate crown area  and m < n
@@ -206,7 +217,7 @@ flora
 
 The Flora object can also be used to show a table of canopy variables:
 
-```{code-cell}
+```{code-cell} ipython3
 # TODO - add a Flora.to_pandas() method
 flora_data = pd.DataFrame({k: getattr(flora, k) for k in flora.trait_attrs})
 flora_data[["name", "ca_ratio", "m", "n", "f_g", "q_m", "z_max_prop"]]
@@ -217,7 +228,7 @@ The T Model uses DBH to define stem size - here the the code is being used to
 back-calculate the required DBH values to give three stems with similar heights
 near the maximum height for each PFT.
 
-```{code-cell}
+```{code-cell} ipython3
 # Generate the expected stem allometries at similar heights for each PFT
 stem_height = np.array([19, 17, 15])
 stem_dbh = calculate_dbh_from_height(
@@ -226,14 +237,14 @@ stem_dbh = calculate_dbh_from_height(
 stem_dbh
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 # Calculate the stem allometries
 allometry = StemAllometry(stem_traits=flora, at_dbh=stem_dbh)
 ```
 
 We can again use {mod}`pandas` to get a table of those allometric predictions:
 
-```{code-cell}
+```{code-cell} ipython3
 pd.DataFrame({k: getattr(allometry, k) for k in allometry.allometry_attrs})
 ```
 
@@ -243,7 +254,7 @@ that is with a shape `(N, 1)`.
 
 We can then calculate the crown profiles.
 
-```{code-cell}
+```{code-cell} ipython3
 # Create a set of vertical heights as a column array.
 z = np.linspace(-1, 20.0, num=211)[:, None]
 
@@ -259,7 +270,7 @@ above calculated at each height $z$:
 * The projected crown area
 * The projected leaf area
 
-```{code-cell}
+```{code-cell} ipython3
 crown_profiles
 ```
 
@@ -276,7 +287,7 @@ stem. For each stem:
 Note that the equation for the relative radius $q(z)$ does define values where
 $z <0$ or $z > H$.
 
-```{code-cell}
+```{code-cell} ipython3
 fig, ax = plt.subplots(ncols=1)
 
 # Find the maximum of the actual and relative maximum crown widths
@@ -318,7 +329,7 @@ We can also use the `CanopyProfile` class with a single row of heights to calcul
 the crown profile at the expected $z_max$ and show that this matches the expected
 crown area from the T Model allometry.
 
-```{code-cell}
+```{code-cell} ipython3
 # Calculate the crown profile across those heights for each PFT
 z_max = flora.z_max_prop * stem_height
 profile_at_zmax = CrownProfile(stem_traits=flora, stem_allometry=allometry, z=z_max)
@@ -334,7 +345,7 @@ using the PFTs defined above because they have very different crown areas, so th
 below generates new profiles for a new set of PFTs that have similar crown area ratios
 but different shapes and gap fractions.
 
-```{code-cell}
+```{code-cell} ipython3
 no_gaps_pft = PlantFunctionalType(
     name="no_gaps", h_max=20, m=1.5, n=1.5, f_g=0, ca_ratio=380
 )
@@ -370,7 +381,7 @@ lines) change with height along the stem.
   lines are identical, but as `f_g` increases, more of the leaf area is displaced down
   within the crown.
 
-```{code-cell}
+```{code-cell} ipython3
 fig, ax = plt.subplots(ncols=1)
 
 for pft_idx, offset, colour in zip((0, 1, 2), (0, 5, 10), ("r", "g", "b")):
@@ -390,7 +401,7 @@ We can also generate predictions for a single PFT with varying crown gap fractio
 the plot below, note that all leaf area is above $z_{max}$ when $f_g=1$ and all leaf
 area is *below*
 
-```{code-cell}
+```{code-cell} ipython3
 fig, ax = plt.subplots(ncols=1)
 
 # Loop over f_g values
@@ -434,6 +445,6 @@ ax.set_xlabel(r"Projected leaf area ($\tilde{A}_{cp}(z)$, m2)")
 ax.legend(frameon=False)
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 
 ```
diff --git a/docs/source/users/demography/flora.md b/docs/source/users/demography/flora.md
index 4b9bf3f7..8db0615d 100644
--- a/docs/source/users/demography/flora.md
+++ b/docs/source/users/demography/flora.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # Plant Functional Types and Traits
@@ -23,7 +34,7 @@ notes and initial demonstration code.
 This page introduces the main components of the {mod}`~pyrealm.demography` module that
 describe plant functional types (PFTs) and their traits.
 
-```{code-cell}
+```{code-cell} ipython3
 from matplotlib import pyplot as plt
 import numpy as np
 import pandas as pd
@@ -100,7 +111,7 @@ their maximum height.
 Note that the `q_m` and `z_max_prop` traits are calculated from the `m` and `n` traits
 and cannot be set directly.
 
-```{code-cell}
+```{code-cell} ipython3
 short_pft = PlantFunctionalType(name="short", h_max=10)
 medium_pft = PlantFunctionalType(name="medium", h_max=20)
 tall_pft = PlantFunctionalType(name="tall", h_max=30)
@@ -108,7 +119,7 @@ tall_pft = PlantFunctionalType(name="tall", h_max=30)
 
 The traits values set for a PFT instance can then be shown:
 
-```{code-cell}
+```{code-cell} ipython3
 short_pft
 ```
 
@@ -129,13 +140,13 @@ that will be used in a demographic simulation. It can be created directly by pro
 the list of {class}`~pyrealm.demography.flora.PlantFunctionalType` instances. The only
 requirement is that each PFT instance uses a different name.
 
-```{code-cell}
+```{code-cell} ipython3
 flora = Flora([short_pft, medium_pft, tall_pft])
 
 flora
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 pd.DataFrame({k: getattr(flora, k) for k in flora.trait_attrs})
 ```
 
@@ -153,7 +164,7 @@ within {class}`~pyrealm.demography.community.Community` objects.
 A `StemTraits` instance can be created directly by providing arrays for each trait, but is
 more easily created from a `Flora` object by providing a list of PFT names:
 
-```{code-cell}
+```{code-cell} ipython3
 # Get stem traits for a range of stems
 stem_pfts = ["short", "short", "short", "medium", "medium", "tall"]
 stem_traits = flora.get_stem_traits(pft_names=stem_pfts)
diff --git a/docs/source/users/demography/module_overview.md b/docs/source/users/demography/module_overview.md
index e2aaceab..45d7f248 100644
--- a/docs/source/users/demography/module_overview.md
+++ b/docs/source/users/demography/module_overview.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # The demography module
diff --git a/docs/source/users/demography/t_model.md b/docs/source/users/demography/t_model.md
index 23510daa..cd38532c 100644
--- a/docs/source/users/demography/t_model.md
+++ b/docs/source/users/demography/t_model.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # The T Model module
@@ -19,7 +30,7 @@ The T Model {cite}`Li:2014bc` provides a model of both:
   (PFT), and
 * a carbon allocation model, given stem allometry and potential GPP.
 
-```{code-cell}
+```{code-cell} ipython3
 from matplotlib import pyplot as plt
 import numpy as np
 import pandas as pd
@@ -31,7 +42,7 @@ from pyrealm.demography.t_model_functions import StemAllocation, StemAllometry
 To generate predictions under the T Model, we need a Flora object providing the
 [trait values](./flora.md) for each of the PFTsto be modelled:
 
-```{code-cell}
+```{code-cell} ipython3
 # Three PFTS
 short_pft = PlantFunctionalType(name="short", h_max=10)
 medium_pft = PlantFunctionalType(name="medium", h_max=20)
@@ -61,14 +72,14 @@ the predictions of the T Model for:
 The DBH input can be a scalar array or a one dimensional array providing a single value
 for each PFT. This then calculates a single estimate at the given size for each stem.
 
-```{code-cell}
+```{code-cell} ipython3
 # Calculate a single prediction
 single_allometry = StemAllometry(stem_traits=flora, at_dbh=np.array([0.1, 0.1, 0.1]))
 ```
 
 We can display those predictions as a `pandas.DataFrame`:
 
-```{code-cell}
+```{code-cell} ipython3
 pd.DataFrame(
     {k: getattr(single_allometry, k) for k in single_allometry.allometry_attrs}
 )
@@ -79,7 +90,7 @@ predictions are made at each DBH value for each PFT and the allometry attributes
 predictions arranged with each PFT as a column and each DBH prediction as a row. This
 makes them convenient to plot using `matplotlib`.
 
-```{code-cell}
+```{code-cell} ipython3
 # Column array of DBH values from 0 to 1.6 metres
 dbh_col = np.arange(0, 1.6, 0.01)[:, None]
 # Get the predictions
@@ -89,7 +100,7 @@ allometries = StemAllometry(stem_traits=flora, at_dbh=dbh_col)
 The code below shows how to use the returned allometries to generate a plot of the
 scaling relationships across all of the PFTs in a `Flora` instance.
 
-```{code-cell}
+```{code-cell} ipython3
 fig, axes = plt.subplots(ncols=2, nrows=4, sharex=True, figsize=(10, 10))
 
 plot_details = [
@@ -118,14 +129,14 @@ The T Model also predicts how potential GPP will be allocated to respiration, tu
 and growth for stems with a given PFT and allometry. Again, a single value can be
 provided to get a single estimate of the allocation model for each stem:
 
-```{code-cell}
+```{code-cell} ipython3
 single_allocation = StemAllocation(
     stem_traits=flora, stem_allometry=single_allometry, at_potential_gpp=np.array([55])
 )
 single_allocation
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 pd.DataFrame(
     {k: getattr(single_allocation, k) for k in single_allocation.allocation_attrs}
 )
@@ -136,14 +147,14 @@ allocation per stem. In the first example, the code takes the allometric predict
 from above and calculates the GPP allocation for stems of varying size with the same
 potential GPP:
 
-```{code-cell}
+```{code-cell} ipython3
 potential_gpp = np.repeat(5, dbh_col.size)[:, None]
 allocation = StemAllocation(
     stem_traits=flora, stem_allometry=allometries, at_potential_gpp=potential_gpp
 )
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 fig, axes = plt.subplots(ncols=2, nrows=5, sharex=True, figsize=(10, 12))
 
 plot_details = [
@@ -175,7 +186,7 @@ fig.delaxes(axes[-1])
 An alternative calculation is to make allocation predictions for varying potential GPP
 for constant allometries:
 
-```{code-cell}
+```{code-cell} ipython3
 # Column array of DBH values from 0 to 1.6 metres
 dbh_constant = np.repeat(0.2, 50)[:, None]
 # Get the allometric predictions
@@ -189,7 +200,7 @@ allocation_2 = StemAllocation(
 )
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 fig, axes = plt.subplots(ncols=2, nrows=5, sharex=True, figsize=(10, 12))
 
 axes = axes.flatten()
@@ -206,6 +217,6 @@ for ax, (var, ylab) in zip(axes, plot_details):
 fig.delaxes(axes[-1])
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 
 ```
diff --git a/docs/source/users/hygro.md b/docs/source/users/hygro.md
index a235875e..bfd2f7fd 100644
--- a/docs/source/users/hygro.md
+++ b/docs/source/users/hygro.md
@@ -5,15 +5,26 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # Hygrometric functions
 
-```{code-cell}
+```{code-cell} ipython3
 # This code loads required packages and then creates a representative range of
 # values of the core variables to use in function plots.
 #
@@ -72,7 +83,7 @@ and returns kPa, so if you are using VP to prepare input data for
 
 ## Saturated vapour pressure
 
-```{code-cell}
+```{code-cell} ipython3
 # Create a sequence of air temperatures and calculate the saturated vapour pressure
 vp_sat = hygro.calc_vp_sat(ta_1d)
 
@@ -85,7 +96,7 @@ pyplot.show()
 
 ## Vapour pressure to VPD
 
-```{code-cell}
+```{code-cell} ipython3
 vpd = hygro.convert_vp_to_vpd(vp_2d, ta_2d.transpose())
 
 # Plot vpd
@@ -100,7 +111,7 @@ pyplot.show()
 
 ## Relative humidity to VPD
 
-```{code-cell}
+```{code-cell} ipython3
 vpd = hygro.convert_rh_to_vpd(rh_2d, ta_2d.transpose())
 
 # Plot vpd
@@ -117,7 +128,7 @@ pyplot.show()
 
 ## Specific humidity to VPD
 
-```{code-cell}
+```{code-cell} ipython3
 # Create a sequence of air temperatures and calculate the saturated vapour pressure
 vpd1 = hygro.convert_sh_to_vpd(sh_1d, ta=20, patm=101.325)
 vpd2 = hygro.convert_sh_to_vpd(sh_1d, ta=30, patm=101.325)
diff --git a/docs/source/users/pmodel/c3c4model.md b/docs/source/users/pmodel/c3c4model.md
index ac4ba204..9630134f 100644
--- a/docs/source/users/pmodel/c3c4model.md
+++ b/docs/source/users/pmodel/c3c4model.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # C3 / C4 Competition
@@ -41,7 +52,7 @@ expected C4 fraction ($F_4$) in a community (see
 correction term for the estimated percentage tree cover and the plot below shows how
 $F_4$ changes with $A_4$, given differing estimates of tree cover.
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 import numpy as np
@@ -87,7 +98,7 @@ The plot below shows how $h$ varies with the expected GPP from C3 plants alone.
 dashed line shows the C3 GPP estimate above which canopy closure leads to complete
 shading of C4 plants.
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 # Just use the competition model to predict h across a GPP gradient
@@ -123,7 +134,7 @@ cover of 0.5.
 
 ### Code
 
-```{code-cell}
+```{code-cell} ipython3
 # Use a simple temperature sequence to generate a range of optimal chi values
 n_pts = 51
 tc_1d = np.linspace(-10, 45, n_pts)
@@ -201,7 +212,7 @@ Panel F
 : The contributions of plants using the C3 and C4 pathways to predicted
   $\delta\ce{^{13}C}$ .
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 # Generate the plots
diff --git a/docs/source/users/pmodel/isotopic_discrimination.md b/docs/source/users/pmodel/isotopic_discrimination.md
index 574b2f67..99bf30b8 100644
--- a/docs/source/users/pmodel/isotopic_discrimination.md
+++ b/docs/source/users/pmodel/isotopic_discrimination.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # Isotopic discrimination
@@ -31,7 +42,7 @@ values of $\chi$. The sequence of $\chi$ values used is created by using the P M
 estimate $\chi$ across a temperature gradient, giving the range of $\chi$ values shown
 below for C3 and C4 plants.
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 import numpy as np
@@ -81,12 +92,12 @@ The calculations differ between C3 and C4 plants, and this is set by the selecti
 the `method_optchi` argument used for the {class}`~pyrealm.pmodel.pmodel.PModel`
 instance.
 
-```{code-cell}
+```{code-cell} ipython3
 carb_c3 = CalcCarbonIsotopes(mod_c3, d13CO2=-8.4, D14CO2=19.2)
 carb_c3.summarize()
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 carb_c4 = CalcCarbonIsotopes(mod_c4, d13CO2=-8.4, D14CO2=19.2)
 carb_c4.summarize()
 ```
@@ -95,7 +106,7 @@ The plots below show how the calculated values alter with $\chi$. The difference
 direction of these relationships between C3 and C4 pathways creates a predictable
 isotopic signature of relative contributions of the two pathways.
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 # Create side by side subplots
diff --git a/docs/source/users/pmodel/module_overview.md b/docs/source/users/pmodel/module_overview.md
index 04c29687..951a560c 100644
--- a/docs/source/users/pmodel/module_overview.md
+++ b/docs/source/users/pmodel/module_overview.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # The P Model module
diff --git a/docs/source/users/pmodel/pmodel_details/envt_variation_outputs.md b/docs/source/users/pmodel/pmodel_details/envt_variation_outputs.md
index e7364837..fc69263a 100644
--- a/docs/source/users/pmodel/pmodel_details/envt_variation_outputs.md
+++ b/docs/source/users/pmodel/pmodel_details/envt_variation_outputs.md
@@ -5,15 +5,26 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # P Model predictions
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 from itertools import product
@@ -115,7 +126,7 @@ environmental variables:
 
 All of the pairwise plots share the same legend:
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 fig, ax = pyplot.subplots(1, 1, figsize=(6, 1.2))
@@ -180,7 +191,7 @@ absorbed irradiance.
 Light use efficiency measures conversion efficiency of moles of absorbed irradiance into
 grams of Carbon ($\mathrm{g\,C}\; \mathrm{mol}^{-1}$ photons).
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 plot_fun("lue", r"LUE ($\mathrm{g\,C}\; \mathrm{mol}^{-1}$ photons).")
@@ -192,7 +203,7 @@ The intrinsic water-use efficiency is ratio of net photosynthetic CO2
 assimilation to stomatal conductance, and captures the cost of assimilation per
 unit of water, in units of $\mu\mathrm{mol}\;\mathrm{mol}^{-1}$.
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 plot_fun("iwue", r"IWUE ($\mu\mathrm{mol}\;\mathrm{mol}^{-1}$)")
@@ -249,7 +260,7 @@ calculated using ``fapar=1, ppfd=1``, which are the default values to
 
 ### Gross primary productivity (``gpp``, GPP)
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 plot_fun("gpp", r"GPP   ($\mu\mathrm{g\,C}\,\mathrm{m}^{-2}\,\mathrm{s}^{-1}$)")
@@ -257,7 +268,7 @@ plot_fun("gpp", r"GPP   ($\mu\mathrm{g\,C}\,\mathrm{m}^{-2}\,\mathrm{s}^{-1}$)")
 
 ### Dark respiration (``rd``)
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 plot_fun("rd", r"$r_d$   ($\mu\mathrm{mol}\,\mathrm{m}^{-2}\,\mathrm{s}^{-1}$)")
@@ -265,7 +276,7 @@ plot_fun("rd", r"$r_d$   ($\mu\mathrm{mol}\,\mathrm{m}^{-2}\,\mathrm{s}^{-1}$)")
 
 ### Maximum rate of carboxylation (``vcmax``)
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 plot_fun("vcmax", r"$v_{cmax}$   ($\mu\mathrm{mol}\,\mathrm{m}^{-2}\,\mathrm{s}^{-1}$)")
@@ -273,7 +284,7 @@ plot_fun("vcmax", r"$v_{cmax}$   ($\mu\mathrm{mol}\,\mathrm{m}^{-2}\,\mathrm{s}^
 
 ### Maximum rate of carboxylation at standard temperature (``vcmax25``)
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 plot_fun(
@@ -283,7 +294,7 @@ plot_fun(
 
 ### Maximum rate of electron transport. (``jmax``)
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 plot_fun("jmax", r"$J_{max}$   ($\mu\mathrm{mol}\,\mathrm{m}^{-2}\,\mathrm{s}^{-1}$)")
@@ -298,7 +309,7 @@ instability in estimates of $g_s$. The
 {meth}`~pyrealm.pmodel.pmodel.PModel.estimate_productivity` method will set $g_s$ to be
 undefined (`np.nan`) when VPD is zero or when $c_a - c_i = 0$.
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 plot_fun("gs", r"$g_s$   ($\mu\mathrm{mol}\,\mathrm{m}^{-2}\,\mathrm{s}^{-1}$)")
@@ -311,7 +322,7 @@ below show how each variable changes, for a constant environment with `tc` of 20
 `patm` of 101325 Pa, `vpd` of 1000 Pa and $\ce{CO2}$ of 400 ppm, when absorbed
 irradiance changes from 0 to 2000 $\mu\text{mol}\,\mathrm{m}^{-2}\,\text{s}^{-1}$.
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 # Calculate the photosynthetic environment
diff --git a/docs/source/users/pmodel/pmodel_details/extreme_values.md b/docs/source/users/pmodel/pmodel_details/extreme_values.md
index 7119ac7f..b4de861e 100644
--- a/docs/source/users/pmodel/pmodel_details/extreme_values.md
+++ b/docs/source/users/pmodel/pmodel_details/extreme_values.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # Extreme forcing values
@@ -51,9 +62,8 @@ settings, the roots of these quadratics are:
 
 Note that the default values for C3 photosynthesis give **non-zero values below 0°C**.
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
-:trusted: true
 
 from matplotlib import pyplot
 import numpy as np
@@ -91,9 +101,8 @@ The photorespiratory compensation point ($\Gamma^*$) varies with as a function o
 temperature and atmospheric pressure, and behaves smoothly with extreme inputs. Note
 that again, $\Gamma^_$ has non-zero values for sub-zero temperatures.
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
-:trusted: true
 
 # Calculate gammastar at different pressures
 tc_1d = np.linspace(-80, 100, n_pts)
@@ -117,9 +126,8 @@ pyplot.show()
 The Michaelis-Menten coefficient for photosynthesis ($K_{mm}$)  also varies with
 temperature and atmospheric pressure and again behaves smoothly with extreme values.
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
-:trusted: true
 
 fig, ax = pyplot.subplots(1, 1)
 
@@ -142,9 +150,8 @@ The density ($\rho$) and viscosity ($\mu$) of water both vary with temperature a
 atmospheric pressure. Looking at the density of water, there is a serious numerical
 issue with low temperatures arising from the equations for the density of water.
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
-:trusted: true
 
 fig, ax = pyplot.subplots(1, 1)
 
@@ -167,9 +174,8 @@ Zooming in, the behaviour of this function is not reliable at extreme low temper
 leading to unstable estimates of $\eta^*$ and the P Model should not be used to make
 predictions below about -30 °C.
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
-:trusted: true
 
 fig, ax = pyplot.subplots(1, 1)
 
diff --git a/docs/source/users/pmodel/pmodel_details/jmax_limitation.md b/docs/source/users/pmodel/pmodel_details/jmax_limitation.md
index 5e07802f..b32d481e 100644
--- a/docs/source/users/pmodel/pmodel_details/jmax_limitation.md
+++ b/docs/source/users/pmodel/pmodel_details/jmax_limitation.md
@@ -5,15 +5,26 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # $J_{max}$ limitation
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 from matplotlib import pyplot
@@ -50,7 +61,7 @@ The plot below shows the effects of each method on the light use efficienct acro
 temperature gradient. The other forcing variables are fixed ($P=101325.0 , \ce{CO2}= 400
 \text{ppm}, \text{VPD}=820$) and $\phi_0$ is also fixed ($\phi_0=0.08$).
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 # Calculate variation in m_jlim with temperature
diff --git a/docs/source/users/pmodel/pmodel_details/optimal_chi.md b/docs/source/users/pmodel/pmodel_details/optimal_chi.md
index fb22b49e..3fdac455 100644
--- a/docs/source/users/pmodel/pmodel_details/optimal_chi.md
+++ b/docs/source/users/pmodel/pmodel_details/optimal_chi.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # Optimal $\chi$ and leaf $\ce{CO2}$
@@ -63,7 +74,7 @@ for use within a P Model.
   - {class}`~pyrealm.pmodel.optimal_chi.OptimalChiC4NoGammaRootzoneStress`
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 from itertools import product
@@ -196,7 +207,7 @@ def plot_opt_chi(mod):
 This **C3 method** follows the approach detailed in {cite:t}`Prentice:2014bc`, see
 {class}`~pyrealm.pmodel.optimal_chi.OptimalChiPrentice14` for details.
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 # Run the P Model and plot predictions
@@ -210,7 +221,7 @@ This **C4 method** follows the approach detailed in {cite:t}`Prentice:2014bc`, b
 a C4 specific version of the unit cost ratio ($\beta$). It also sets $m_j = m_c = 1$.
 See {class}`~pyrealm.pmodel.optimal_chi.OptimalChiC4` for details.
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 # Run the P Model and plot predictions
@@ -228,7 +239,7 @@ and also also sets $m_j = 1$, but $m_c$ is calculated as in
 {class}`~pyrealm.pmodel.optimal_chi.OptimalChiPrentice14`. See
 {meth}`~pyrealm.pmodel.optimal_chi.OptimalChiC4NoGamma` for details.
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 # Run the P Model and plot predictions
@@ -263,7 +274,7 @@ The calculation details are provided in the description of the
 {class}`~pyrealm.pmodel.optimal_chi.OptimalChiLavergne20C3` method, but the
 variation in $\beta$ with $\theta$ is shown below.
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 # Theta is required for the calculation of beta
@@ -291,7 +302,7 @@ The plots below show the impacts on optimal $\chi$ across a temperature gradient
 values of VPD and soil moisture, with constant atmospheric pressure (101325 Pa) and CO2
 (280 ppm).
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 # Environments with high and low soil moisture
@@ -356,7 +367,7 @@ pyplot.tight_layout()
 The plots below illustrate the impact of temperature and  $\theta$ on  $m_j$ and $m_c$,
 again with constant atmospheric pressure (101325 Pa) and CO2 (280 ppm).
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 fig, ((ax1, ax3), (ax2, ax4)) = pyplot.subplots(2, 2, figsize=(10, 10), sharey=True)
@@ -454,7 +465,7 @@ but the variation in $\beta$ with rootzone stress is shown below.
 * {class}`~pyrealm.pmodel.optimal_chi.OptimalChiC4RootzoneStress`
 * {class}`~pyrealm.pmodel.optimal_chi.OptimalChiC4NoGammaRootzoneStress`
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 from pyrealm.pmodel.optimal_chi import (
@@ -500,7 +511,7 @@ The plots below show the impacts on optimal $\chi$ across a temperature gradient
 values of VPD and rootzone stress, with constant atmospheric pressure (101325 Pa) and CO2
 (280 ppm).
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 # Environments with high and low rootzone stress
diff --git a/docs/source/users/pmodel/pmodel_details/photosynthetic_environment.md b/docs/source/users/pmodel/pmodel_details/photosynthetic_environment.md
index cedd27fe..9e2c902e 100644
--- a/docs/source/users/pmodel/pmodel_details/photosynthetic_environment.md
+++ b/docs/source/users/pmodel/pmodel_details/photosynthetic_environment.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # Photosynthetic environment
@@ -26,7 +37,7 @@ The descriptions below show the typical ranges of these values under common
 environmental inputs along with links to the more detailed documentation of
 the key functions.
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 # This code loads required packages and then creates a representative range of
@@ -68,7 +79,7 @@ Details: {func}`pyrealm.pmodel.functions.calc_gammastar`
 The photorespiratory compensation point ($\Gamma^*$) varies with as a function
 of temperature and atmospheric pressure:
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 # Calculate gammastar
@@ -91,7 +102,7 @@ Details: {func}`pyrealm.pmodel.functions.calc_kmm`
 The Michaelis-Menten coefficient for photosynthesis ($K_{mm}$) also varies with
 temperature and atmospheric pressure:
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 # Calculate K_mm
@@ -118,7 +129,7 @@ pressure ($\eta^*$).
 
 The figure shows how $\eta^*$ varies with temperature and pressure.
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 # Calculate the viscosity under the range of values and the standard
@@ -146,7 +157,7 @@ Details: {func}`pyrealm.pmodel.functions.calc_co2_to_ca`
 The partial pressure of $\ce{CO2}$ is a function of the atmospheric concentration of
 $\ce{CO2}$ in parts per million and the atmospheric pressure:
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 # Variation in partial pressure
diff --git a/docs/source/users/pmodel/pmodel_details/pmodel_overview.md b/docs/source/users/pmodel/pmodel_details/pmodel_overview.md
index 9626ee98..475f9d6a 100644
--- a/docs/source/users/pmodel/pmodel_details/pmodel_overview.md
+++ b/docs/source/users/pmodel/pmodel_details/pmodel_overview.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 <!-- markdownlint-disable-next-line MD041 -->
diff --git a/docs/source/users/pmodel/pmodel_details/quantum_yield.md b/docs/source/users/pmodel/pmodel_details/quantum_yield.md
index 92a3ae8c..d5162a1f 100644
--- a/docs/source/users/pmodel/pmodel_details/quantum_yield.md
+++ b/docs/source/users/pmodel/pmodel_details/quantum_yield.md
@@ -5,15 +5,26 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # Quantum yield efficiency of photosynthesis
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 # This code loads required packages and then creates a representative range of
@@ -79,7 +90,7 @@ estimate of $\phi_0$, following {cite:t}`Bernacchi:2003dc` for C3 plants and
 $\phi_0 = 0.081785$, following the BRC parameterisation in Table 1. of
 {cite:t}`Stocker:2020dh`.
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 # Calculate temperature dependence of quantum yield efficiency
@@ -117,7 +128,7 @@ $\phi_0$ values to an otherwise constant environment. As you would expect given
 $\text{LUE} = \phi_0 \cdot M_C \cdot m_j$, light use efficiency changes linearly along
 this gradient of $\phi_0$ values.
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 # A constant environment to show a range of kphio values
@@ -157,7 +168,7 @@ $\phi_{0A} = \dfrac{\phi_{0R}}{(1 + \textrm{AI}^m) ^ n}$
 
 This captures a decrease in maximum $\phi_0$ in arid conditions, as shown below.
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 n_vals = 51
@@ -188,7 +199,7 @@ approach also alters the temperature at which $\phi_0$ is maximised as a functio
 mean growth temperature ($T_g$) in a location. The plot below shows how aridity and mean
 growth temperature interact to change the location and height of the peak $\phi_0$.
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 n_vals = 51
diff --git a/docs/source/users/pmodel/pmodel_details/rpmodel.md b/docs/source/users/pmodel/pmodel_details/rpmodel.md
index 7e691cc2..c46c12ba 100644
--- a/docs/source/users/pmodel/pmodel_details/rpmodel.md
+++ b/docs/source/users/pmodel/pmodel_details/rpmodel.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # The `rpmodel` implementation
diff --git a/docs/source/users/pmodel/pmodel_details/soil_moisture.md b/docs/source/users/pmodel/pmodel_details/soil_moisture.md
index 2b4b8a53..8d94a557 100644
--- a/docs/source/users/pmodel/pmodel_details/soil_moisture.md
+++ b/docs/source/users/pmodel/pmodel_details/soil_moisture.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # Soil moisture effects
@@ -85,7 +96,7 @@ varies with changing soil moisture for some different values of mean aridity. In
 the examples below, the default $\theta_0 = 0$ has been changed to $\theta_0 =
 0.1$ to make the lower bound more obvious.
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 from matplotlib import pyplot as plt
@@ -151,7 +162,7 @@ by the resulting factor. The example below shows how the predicted light use
 efficiency from the P Model changes across an aridity gradient both with and without the
 soil moisture factor.
 
-```{code-cell}
+```{code-cell} ipython3
 # Calculate the P Model in a constant environment
 tc = np.array([20] * 101)
 sm_gradient = np.linspace(0, 1.0, 101)
@@ -174,7 +185,7 @@ for mean_alpha in [0.9, 0.5, 0.3, 0.1, 0.0]:
     gpp_stressed[mean_alpha] = model.gpp * sm_stress
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 plt.plot(sm_gradient, model.gpp, label="No soil moisture penalty")
@@ -234,7 +245,7 @@ y &= \min( a  \textrm{AI} ^ {b}, 1)\\
 \end{align*}
 $$
 
-```{code-cell}
+```{code-cell} ipython3
 from pyrealm.constants import PModelConst
 
 const = PModelConst()
@@ -269,7 +280,7 @@ $$
     \end{cases}
 $$
 
-```{code-cell}
+```{code-cell} ipython3
 # Calculate the soil moisture stress factor across a soil moisture
 # gradient for different aridity index values
 beta = {}
@@ -298,7 +309,7 @@ calculated and then applied to the GPP calculated for a model
 ({attr}`~pyrealm.pmodel.pmodel.PModel.gpp`). In the example below, the result is
 obviously just $\beta(\theta)$ from above scaled to the constant GPP.
 
-```{code-cell}
+```{code-cell} ipython3
 for ai in ai_vals:
 
     plt.plot(sm_gradient, model.gpp * beta[ai], label=f"AI = {ai}")
@@ -309,6 +320,6 @@ plt.legend()
 plt.show()
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 
 ```
diff --git a/docs/source/users/pmodel/pmodel_details/worked_examples.md b/docs/source/users/pmodel/pmodel_details/worked_examples.md
index f3db6947..f69cc924 100644
--- a/docs/source/users/pmodel/pmodel_details/worked_examples.md
+++ b/docs/source/users/pmodel/pmodel_details/worked_examples.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # Worked examples
@@ -38,9 +49,7 @@ The example shows the steps required using a single site with:
 
 ### Estimate photosynthetic environment
 
-```{code-cell}
-:trusted: true
-
+```{code-cell} ipython3
 from importlib import resources
 
 from matplotlib import pyplot as plt
@@ -60,15 +69,11 @@ terse - just the shape of the data - but the
 {meth}`~pyrealm.pmodel.pmodel_environment.PModelEnvironment.summarize` method provides a
 more detailed summary of the attributes.
 
-```{code-cell}
-:trusted: true
-
+```{code-cell} ipython3
 env
 ```
 
-```{code-cell}
-:trusted: true
-
+```{code-cell} ipython3
 env.summarize()
 ```
 
@@ -77,18 +82,14 @@ env.summarize()
 Next, the P Model can be fitted to the photosynthetic environment using the
 ({class}`~pyrealm.pmodel.pmodel.PModel`) class:
 
-```{code-cell}
-:trusted: true
-
+```{code-cell} ipython3
 model = PModel(env)
 ```
 
 The returned model object holds a lot of information. The representation of the
 model object shows a terse display of the settings used to run the model:
 
-```{code-cell}
-:trusted: true
-
+```{code-cell} ipython3
 model
 ```
 
@@ -98,9 +99,7 @@ displays a summary of calculated predictions. Initially, this shows two measures
 photosynthetic efficiency: the intrinsic water use efficiency (``iwue``) and the light
 use efficiency (``lue``).
 
-```{code-cell}
-:trusted: true
-
+```{code-cell} ipython3
 model.summarize()
 ```
 
@@ -112,9 +111,7 @@ recording key parameters from the [calculation of $\chi$](./optimal_chi).
 This object also has a {meth}`~pyrealm.pmodel.optimal_chi.OptimalChiABC.summarize`
 method:
 
-```{code-cell}
-:trusted: true
-
+```{code-cell} ipython3
 model.optchi.summarize()
 ```
 
@@ -130,9 +127,7 @@ Here we are using:
 * An absorption fraction of 0.91 (-), and
 * a PPFD of 834 µmol m-2 s-1.
 
-```{code-cell}
-:trusted: true
-
+```{code-cell} ipython3
 model.estimate_productivity(fapar=0.91, ppfd=834)
 model.summarize()
 ```
@@ -164,9 +159,7 @@ to be the same size so some of the variables have repeated data across dimension
   cell.
 * Elevation is constant across months, so the data for each month is repeated.
 
-```{code-cell}
-:trusted: true
-
+```{code-cell} ipython3
 # Load an example dataset containing the forcing variables.
 data_path = resources.files("pyrealm_build_data.rpmodel") / "pmodel_global.nc"
 ds = xarray.load_dataset(data_path)
@@ -185,9 +178,7 @@ The model can now be run using that data. The first step is to convert the eleva
 data to atmospheric pressure, and then this is used to set the photosynthetic
 environment for the model:
 
-```{code-cell}
-:trusted: true
-
+```{code-cell} ipython3
 # Convert elevation to atmospheric pressure
 patm = calc_patm(elev)
 
@@ -205,9 +196,7 @@ env.summarize()
 That environment can then be run to calculate the P model predictions for light use
 efficiency:
 
-```{code-cell}
-:trusted: true
-
+```{code-cell} ipython3
 # Run the P model
 model = PModel(env)
 
@@ -220,9 +209,7 @@ plt.title("Light use efficiency")
 Finally, the light use efficiency can be used to calculate GPP given the
 photosynthetic photon flux density and fAPAR.
 
-```{code-cell}
-:trusted: true
-
+```{code-cell} ipython3
 # Scale the outputs from values per unit iabs to realised values
 model.estimate_productivity(fapar, ppfd)
 
diff --git a/docs/source/users/pmodel/subdaily_details/acclimation.md b/docs/source/users/pmodel/subdaily_details/acclimation.md
index 193b38f9..31805af1 100644
--- a/docs/source/users/pmodel/subdaily_details/acclimation.md
+++ b/docs/source/users/pmodel/subdaily_details/acclimation.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # Estimating acclimation
@@ -28,7 +39,7 @@ modelling approach to representing slow responses within the P Model, following
 
 * The interpolation of realised daily values back onto the subdaily timescale.
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 from copy import copy
@@ -58,7 +69,7 @@ interpolating data back to subdaily timescales. In practice {cite:t}`mengoli:202
 present results using one hour windows around noon or even the single value closest to
 noon.
 
-```{code-cell}
+```{code-cell} ipython3
 # Define a set of observations at a subdaily timescale
 fast_datetimes = np.arange(
     np.datetime64("1970-01-01"), np.datetime64("1970-01-08"), np.timedelta64(30, "m")
@@ -78,7 +89,7 @@ demo_scaler.set_window(window_center=np.timedelta64(12, "h"), half_width=half_wi
 The plot below shows the rapidly changing variable and the defined daily acclimation
 windows.
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 fig, ax = plt.subplots()
@@ -129,7 +140,7 @@ applies the memory effect with three different values of $\alpha$. When $\alpha
 the realised values are identical to the daily optimum value within the acclimation
 window.
 
-```{code-cell}
+```{code-cell} ipython3
 # Extract the optimal values within the daily acclimation windows
 daily_mean = demo_scaler.get_daily_means(fast_data)
 
@@ -139,7 +150,7 @@ real_3 = memory_effect(daily_mean, alpha=1 / 3)
 real_1 = memory_effect(daily_mean, alpha=1)
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 fig, ax = plt.subplots()
@@ -194,7 +205,7 @@ The code below shows how the
 used to interpolate realised values back to the subdaily scale, using different settings
 for the update point and interpolation method.
 
-```{code-cell}
+```{code-cell} ipython3
 # Fill to the subdaily scale using the default settings:
 # - update at the end of the acclimation window
 # - hold the value constant between update points
@@ -248,7 +259,7 @@ Plot D
 : The daily optimal realised value is again able to instantaneously adopt the daily
   optimal value, but the one day offset for linear interpolation is applied.
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 # Create the figure
diff --git a/docs/source/users/pmodel/subdaily_details/subdaily_calculations.md b/docs/source/users/pmodel/subdaily_details/subdaily_calculations.md
index 1a134c84..f32823a8 100644
--- a/docs/source/users/pmodel/subdaily_details/subdaily_calculations.md
+++ b/docs/source/users/pmodel/subdaily_details/subdaily_calculations.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # Subdaily P Model calculations
@@ -19,9 +30,8 @@ steps used in the estimation process in order to show intermediates results but
 practice, as shown in the [worked example](worked_example), most of these calculations
 are handled internally by the model fitting in `pyrealm`.
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
-:trusted: true
 
 from importlib import resources
 
@@ -47,9 +57,7 @@ The code below uses half hourly data from 2014 for the [BE-Vie FluxNET
 site](https://fluxnet.org/doi/FLUXNET2015/BE-Vie), which was also used as a
 demonstration in {cite:t}`mengoli:2022a`.
 
-```{code-cell}
-:trusted: true
-
+```{code-cell} ipython3
 data_path = resources.files("pyrealm_build_data.subdaily") / "subdaily_BE_Vie_2014.csv"
 
 data = pandas.read_csv(str(data_path))
@@ -70,9 +78,7 @@ This dataset can then be used to calculate the photosynthetic environment at the
 subdaily timescale. The code below also estimates GPP under the standard P Model with no
 slow responses for comparison.
 
-```{code-cell}
-:trusted: true
-
+```{code-cell} ipython3
 # Calculate the photosynthetic environment
 subdaily_env = PModelEnvironment(
     tc=temp_subdaily,
@@ -98,9 +104,7 @@ best to sample those conditions. Typically those might be the observed environme
 conditions at the observation closest to noon, or the mean environmental conditions in a
 window around noon.
 
-```{code-cell}
-:trusted: true
-
+```{code-cell} ipython3
 # Create the fast slow scaler
 fsscaler = SubdailyScaler(datetime_subdaily)
 
@@ -120,9 +124,8 @@ pmodel_subdaily = SubdailyPModel(
 )
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
-:trusted: true
 
 idx = np.arange(48 * 120, 48 * 130)
 plt.figure(figsize=(10, 4))
@@ -145,9 +148,7 @@ The daily average conditions during the acclimation window can be sampled and us
 inputs to the standard P Model to calculate the optimal behaviour of plants under those
 conditions.
 
-```{code-cell}
-:trusted: true
-
+```{code-cell} ipython3
 # Get the daily acclimation conditions for the forcing variables
 temp_acclim = fsscaler.get_daily_means(temp_subdaily)
 co2_acclim = fsscaler.get_daily_means(co2_subdaily)
@@ -178,9 +179,7 @@ temperatures so $J_{max}$ and $V_{cmax}$ must first be standardised to expected
 at 25°C. This is acheived by multiplying by the reciprocal of the exponential part of
 the Arrhenius equation ($h^{-1}$ in {cite}`mengoli:2022a`).
 
-```{code-cell}
-:trusted: true
-
+```{code-cell} ipython3
 # Are these any of the existing values in the constants?
 ha_vcmax25 = 65330
 ha_jmax25 = 43900
@@ -195,9 +194,7 @@ jmax25_acclim = pmodel_acclim.jmax * (1 / calc_ftemp_arrh(tk_acclim, ha_jmax25))
 The memory effect can now be applied to the three parameters with slow
 responses to calculate realised values, here using the default 15 day window.
 
-```{code-cell}
-:trusted: true
-
+```{code-cell} ipython3
 # Calculation of memory effect in xi, vcmax25 and jmax25
 xi_real = memory_effect(pmodel_acclim.optchi.xi, alpha=1 / 15)
 vcmax25_real = memory_effect(vcmax25_acclim, alpha=1 / 15, allow_holdover=True)
@@ -208,9 +205,8 @@ The plots below show the instantaneously acclimated values for  $J_{max25}$,
 $V_{cmax25}$ and $\xi$ in grey along with the realised slow reponses, after
 application of the memory effect.
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
-:trusted: true
 
 fig, axes = plt.subplots(1, 3, figsize=(16, 5))
 
@@ -243,9 +239,7 @@ temperature at fast scales:
 * These values are adjusted to the actual half hourly temperatures to give the fast
   responses of $J_{max}$ and $V_{cmax}$.
 
-```{code-cell}
-:trusted: true
-
+```{code-cell} ipython3
 tk_subdaily = subdaily_env.tc + pmodel_subdaily.env.core_const.k_CtoK
 
 # Fill the realised jmax and vcmax from subdaily to daily
@@ -265,9 +259,7 @@ passing the realised values of $\xi$ as a fixed constraint to the calculation of
 optimal $\chi$, rather than calculating the instantaneously optimal values of $\xi$
 as is the case in the standard P Model.
 
-```{code-cell}
-:trusted: true
-
+```{code-cell} ipython3
 # Interpolate xi to subdaily scale
 xi_subdaily = fsscaler.fill_daily_to_subdaily(xi_real)
 
@@ -286,9 +278,7 @@ Model, where $c_i$ includes the slow responses of $\xi$ and $V_{cmax}$ and $J_{m
 include the slow responses of $V_{cmax25}$ and $J_{max25}$ and fast responses to
 temperature.
 
-```{code-cell}
-:trusted: true
-
+```{code-cell} ipython3
 # Calculate Ac
 Ac_subdaily = (
     vcmax_subdaily
@@ -318,8 +308,6 @@ diff = GPP_subdaily - pmodel_subdaily.gpp
 print(np.nanmin(diff), np.nanmax(diff))
 ```
 
-```{code-cell}
-:trusted: true
-
+```{code-cell} ipython3
 
 ```
diff --git a/docs/source/users/pmodel/subdaily_details/subdaily_model_and_missing_data.md b/docs/source/users/pmodel/subdaily_details/subdaily_model_and_missing_data.md
index 881f71b8..61d730d8 100644
--- a/docs/source/users/pmodel/subdaily_details/subdaily_model_and_missing_data.md
+++ b/docs/source/users/pmodel/subdaily_details/subdaily_model_and_missing_data.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # Missing data in the subdaily model
@@ -46,7 +57,7 @@ The code below gives a concrete example - a time series that starts and ends dur
 the middle of a one hour acclimation window around noon. Only two of the three
 observations are provided for the first and last day
 
-```{code-cell}
+```{code-cell} ipython3
 import numpy as np
 
 from pyrealm.pmodel.scaler import SubdailyScaler
@@ -88,7 +99,7 @@ problem of the missing data clearly:
 * One day has a single missing 12:00 data point within the acclimation window.
 * One day has no data within the acclimation window.
 
-```{code-cell}
+```{code-cell} ipython3
 fsscaler.get_window_values(data)
 ```
 
@@ -97,7 +108,7 @@ The daily average conditions are calculated using the
 partial data are not allowed - which is the default - the daily average conditions for
 all days with missing data is also missing (`np.nan`).
 
-```{code-cell}
+```{code-cell} ipython3
 partial_not_allowed = fsscaler.get_daily_means(data)
 partial_not_allowed
 ```
@@ -106,7 +117,7 @@ Setting `allow_partial_data = True` allows the daily average conditions to be ca
 from the partial available information. This does not solve the problem for the day with
 no data in the acclimation window, which still results in a missing value.
 
-```{code-cell}
+```{code-cell} ipython3
 partial_allowed = fsscaler.get_daily_means(data, allow_partial_data=True)
 partial_allowed
 ```
@@ -115,7 +126,7 @@ The :func:`~pyrealm.pmodel.subdaily.memory_effect` function is used to calculate
 realised values of a variable from the optimal values. By default, this function *will
 raise an error* when missing data are present:
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [raises-exception]
 
 memory_effect(partial_not_allowed)
@@ -125,14 +136,14 @@ The `allow_holdover` option allows the function to be run - the value for the fi
 is still `np.nan` but the missing observations on day 3, 5 and 7 are filled by holding
 over the valid observations from the previous day.
 
-```{code-cell}
+```{code-cell} ipython3
 memory_effect(partial_not_allowed, allow_holdover=True)
 ```
 
 When the partial data is allowed, the `allow_holdover` is still required to fill the
 gap on day 5 by holding over the data from day 4.
 
-```{code-cell}
+```{code-cell} ipython3
 memory_effect(partial_allowed, allow_holdover=True)
 ```
 
diff --git a/docs/source/users/pmodel/subdaily_details/subdaily_overview.md b/docs/source/users/pmodel/subdaily_details/subdaily_overview.md
index 3278b056..e96a57c9 100644
--- a/docs/source/users/pmodel/subdaily_details/subdaily_overview.md
+++ b/docs/source/users/pmodel/subdaily_details/subdaily_overview.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # The P Model with acclimation
diff --git a/docs/source/users/pmodel/subdaily_details/worked_example.md b/docs/source/users/pmodel/subdaily_details/worked_example.md
index fac12b6a..f30d4155 100644
--- a/docs/source/users/pmodel/subdaily_details/worked_example.md
+++ b/docs/source/users/pmodel/subdaily_details/worked_example.md
@@ -5,17 +5,26 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # Worked example of the Subdaily P Model
 
-```{code-cell}
-:trusted: true
-
+```{code-cell} ipython3
 from importlib import resources
 
 import xarray
@@ -60,9 +69,7 @@ fitting basically takes all of the same arguments as the standard
 The test data use some UK WFDE data for three sites in order to compare predictions
 over a time series.
 
-```{code-cell}
-:trusted: true
-
+```{code-cell} ipython3
 # Loading the example dataset:
 dpath = (
     resources.files("pyrealm_build_data.uk_data") / "UK_WFDE5_FAPAR_2018_JuneJuly.nc"
@@ -83,9 +90,7 @@ sites = xarray.Dataset(
 The WFDE data need some conversion for use in the PModel, along with the definition of
 the atmospheric CO2 concentration.
 
-```{code-cell}
-:trusted: true
-
+```{code-cell} ipython3
 # Variable set up
 # Air temperature in °C from Tair in Kelvin
 tc = (ds["Tair"] - 273.15).to_numpy()
@@ -104,9 +109,7 @@ co2 = np.ones_like(tc) * 400
 
 The code below then calculates the photosynthetic environment.
 
-```{code-cell}
-:trusted: true
-
+```{code-cell} ipython3
 # Generate and check the PModelEnvironment
 pm_env = PModelEnvironment(tc=tc, patm=patm, vpd=vpd, co2=co2)
 pm_env.summarize()
@@ -117,9 +120,7 @@ pm_env.summarize()
 The standard implementation of the P Model used below assumes that plants can
 instantaneously adopt optimal behaviour.
 
-```{code-cell}
-:trusted: true
-
+```{code-cell} ipython3
 # Standard PModels
 pmodC3 = PModel(
     env=pm_env, method_kphio="fixed", reference_kphio=1 / 8, method_optchi="prentice14"
@@ -128,9 +129,7 @@ pmodC3.estimate_productivity(fapar=fapar, ppfd=ppfd)
 pmodC3.summarize()
 ```
 
-```{code-cell}
-:trusted: true
-
+```{code-cell} ipython3
 pmodC4 = PModel(
     env=pm_env, method_kphio="fixed", reference_kphio=1 / 8, method_optchi="c4_no_gamma"
 )
@@ -146,9 +145,7 @@ values to holdover previous realised values to cover missing data within the
 calculations: essentially the plant does not acclimate until the optimal values can be
 calculated again to update those realised estimates.
 
-```{code-cell}
-:trusted: true
-
+```{code-cell} ipython3
 # Set the acclimation window to an hour either side of noon
 fsscaler = SubdailyScaler(datetimes)
 fsscaler.set_window(
@@ -187,9 +184,7 @@ The code below then extracts the time series for the two months from the three s
 shown above and plots the instantaneous predictions against predictions including slow
 photosynthetic responses.
 
-```{code-cell}
-:trusted: true
-
+```{code-cell} ipython3
 # Store the predictions in the xarray Dataset to use indexing
 ds["GPP_pmodC3"] = (ds["Tair"].dims, pmodC3.gpp)
 ds["GPP_subdailyC3"] = (ds["Tair"].dims, subdailyC3.gpp)
@@ -246,9 +241,7 @@ plt.tight_layout()
 The subdaily models can also be obtained directly from the standard models, using the
 `convert_pmodel_to_subdaily` method:
 
-```{code-cell}
-:trusted: true
-
+```{code-cell} ipython3
 # Convert standard C3 model
 converted_C3 = convert_pmodel_to_subdaily(
     pmodel=pmodC3,
@@ -269,9 +262,7 @@ converted_C4 = convert_pmodel_to_subdaily(
 This produces the same outputs as the `SubdailyPModel` class, but is convenient and more
 compact when the two models are going to be compared.
 
-```{code-cell}
-:trusted: true
-
+```{code-cell} ipython3
 # Models have identical GPP - maximum absolute difference is zero.
 print(np.nanmax(abs(subdailyC3.gpp.flatten() - converted_C3.gpp.flatten())))
 print(np.nanmax(abs(subdailyC4.gpp.flatten() - converted_C4.gpp.flatten())))
diff --git a/docs/source/users/splash.md b/docs/source/users/splash.md
index e5c4d283..20236db4 100644
--- a/docs/source/users/splash.md
+++ b/docs/source/users/splash.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # The `splash` submodule
@@ -46,7 +57,7 @@ The data below provides a 2 year daily time series of precipitation, temperature
 solar fraction (1 - cloud cover) for 0.5° resolution grid cells in a 10° by 10° block
 of the North Western USA. It also provides the mean elevation of those cells.
 
-```{code-cell}
+```{code-cell} ipython3
 from importlib import resources
 import numpy as np
 import xarray
@@ -74,7 +85,7 @@ data
 The plot below shows the elevation for the example data area, along with the locations
 of three sites that will be used to compare SPLASH outputs.
 
-```{code-cell}
+```{code-cell} ipython3
 # Get the latitude and longitude extents
 extent = (
     data["lon"].min(),
@@ -95,7 +106,7 @@ The three sites capture wetter coastal conditions with milder temperatures (San
 Francisco), intermediate rainfall with colder temperatures (Yosemite) and arid
 conditions with extreme temperatures (Death Valley).
 
-```{code-cell}
+```{code-cell} ipython3
 # Get three sites to show time series for locations
 site_data = data.sel(sites, method="nearest")
 
@@ -149,7 +160,7 @@ may well be constant across the longitude dimension for gridded data - but, at t
 moment, you need to broadcast these variables to match.
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 splash = SplashModel(
     lat=np.broadcast_to(data.lat.data[None, :, None], data.sf.data.shape),
     elv=np.broadcast_to(data.elev.data[None, :, :], data.sf.data.shape),
@@ -176,7 +187,7 @@ give the expected soil moisture at the end of the year. If this is sufficiently
 to the start values, the estimate is returned, otherwise the end of year expectations
 are used as a starting point to recalculate the annual water balances.
 
-```{code-cell}
+```{code-cell} ipython3
 init_soil_moisture = splash.estimate_initial_soil_moisture(verbose=False)
 ```
 
@@ -200,7 +211,7 @@ The plots show the soil moisture for the first day, along with the changes in so
 moisture from the initial estimates (the 'previous day'). Note the saturated soil
 moisture of 150mm near the coast and in the mountains.
 
-```{code-cell}
+```{code-cell} ipython3
 # Calculate the water balance equation for the first day from the initial soil
 #  moisture estimates.
 aet, wn, ro = splash.estimate_daily_water_balance(init_soil_moisture, day_idx=0)
@@ -223,13 +234,13 @@ the daily estimation across all of the dates in the input data from initial soil
 moisture estimates. It returns a set of time series of soil moisture, runoff and AET for
 all sites.
 
-```{code-cell}
+```{code-cell} ipython3
 aet_out, wn_out, ro_out = splash.calculate_soil_moisture(init_soil_moisture)
 ```
 
 The plots below show the resulting soil moisture and a time series for the three
 
-```{code-cell}
+```{code-cell} ipython3
 # Add the outputs to the xarray to select the three sites easily.
 data["aet"] = xarray.DataArray(aet_out, dims=("time", "lat", "lon"))
 data["wn"] = xarray.DataArray(wn_out, dims=("time", "lat", "lon"))
diff --git a/docs/source/users/tmodel/canopy.md b/docs/source/users/tmodel/canopy.md
index 3fab2df8..d760ca77 100644
--- a/docs/source/users/tmodel/canopy.md
+++ b/docs/source/users/tmodel/canopy.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # Canopy model
@@ -36,7 +47,7 @@ This sketch:
   object.
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 import numpy as np
 import matplotlib.pyplot as plt
 from scipy.optimize import root_scalar
@@ -54,11 +65,11 @@ The scaling of a set of trees is automatically calculated using the initial diam
 the `TTree` instance. This automatically calculates the other dimensions, such as
 height, using the underlying scaling equations of the T Model.
 
-```{code-cell}
+```{code-cell} ipython3
 pft.height
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 :lines_to_next_cell: 2
 
 pft.crown_area
@@ -97,7 +108,7 @@ r_0 &= \frac{1}{q_m}\sqrt{\frac{A_c}{\pi}}
 \end{align}
 $$
 
-```{code-cell}
+```{code-cell} ipython3
 :lines_to_next_cell: 2
 
 def calculate_qm(m, n):
@@ -146,7 +157,7 @@ r(z) &= r_0 \; q(z)
 \end{align}
 $$
 
-```{code-cell}
+```{code-cell} ipython3
 def calculate_relative_canopy_radius_at_z(z, H, m, n):
     """Calculate q(z)"""
 
@@ -155,20 +166,20 @@ def calculate_relative_canopy_radius_at_z(z, H, m, n):
     return m * n * z_over_H ** (n - 1) * (1 - z_over_H**n) ** (m - 1)
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 # Validate that zm and r0 generate the predicted maximum crown area
 q_zm = calculate_relative_canopy_radius_at_z(zm, pft.height, m, n)
 rm = r0 * q_zm
 print("rm = ", rm)
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 np.allclose(rm**2 * np.pi, pft.crown_area)
 ```
 
 Vertical crown radius profiles can now be calculated for each stem:
 
-```{code-cell}
+```{code-cell} ipython3
 # Create an interpolation from ground to maximum stem height, with 5 cm resolution.
 # Also append a set of values _fractionally_ less than the exact height  of stems
 # so that the height at the top of each stem is included but to avoid floating
@@ -190,7 +201,7 @@ np.cumsum(np.convolve(rm, np.ones(2), "valid") + 0.1)
 
 Those can be plotted out to show the vertical crown radius profiles
 
-```{code-cell}
+```{code-cell} ipython3
 # Separate the stems along the x axis for plotting
 stem_x = np.concatenate(
     [np.array([0]), np.cumsum(np.convolve(rm, np.ones(2), "valid") + 0.4)]
@@ -237,7 +248,7 @@ A_c \left(\dfrac{q(z)}{q_m}\right)^2, & H > z > z_m \\
 \end{cases}
 $$
 
-```{code-cell}
+```{code-cell} ipython3
 Stems = float | np.ndarray
 
 
@@ -281,7 +292,7 @@ def calculate_projected_area(
 The code below calculates the projected crown area for each stem and then plots the
 vertical profile for individual stems and across the community.
 
-```{code-cell}
+```{code-cell} ipython3
 :lines_to_next_cell: 2
 
 # Calculate the projected area for each stem
@@ -331,7 +342,7 @@ $$
 l_m = \left\lceil \frac{\sum_1^{N_s}{ A_c}}{ A(1 - f_G)}\right\rceil
 $$
 
-```{code-cell}
+```{code-cell} ipython3
 def solve_canopy_closure_height(
     z: float,
     l: int,
@@ -395,7 +406,7 @@ def calculate_canopy_heights(
 The example below calculates the projected crown area above ground level for the example
 stems. These should be identical to the crown area of the stems.
 
-```{code-cell}
+```{code-cell} ipython3
 # Set the total available canopy space and community gap fraction
 canopy_area = 32
 community_gap_fraction = 2 / 32
@@ -412,7 +423,7 @@ superimpose the calculated $z^*_l$ values and the cumulative canopy area for eac
 to confirm that the calculated values coincide with the profile. Note here that the
 total area at each closed layer height is omitting the community gap fraction.
 
-```{code-cell}
+```{code-cell} ipython3
 community_Ap_z = np.nansum(Ap_z, axis=1)
 
 fig, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(10, 5))
@@ -481,14 +492,14 @@ plt.tight_layout()
 The projected area from individual stems to each canopy layer can then be calculated at
 $z^*_l$ and hence the projected area of canopy **within each layer**.
 
-```{code-cell}
+```{code-cell} ipython3
 # Calculate the canopy area above z_star for each stem
 Ap_z_star = calculate_projected_area(z=z_star[:, None], pft=pft, m=m, n=n, qm=qm, zm=zm)
 
 print(Ap_z_star)
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 :lines_to_next_cell: 2
 
 # Calculate the contribution _within_ each layer per stem
@@ -527,7 +538,7 @@ $$
 
 The function below calculates $\tilde{A}_{cp}(z)$.
 
-```{code-cell}
+```{code-cell} ipython3
 def calculate_leaf_area(
     z: float,
     fg: float,
@@ -572,7 +583,7 @@ there are no crown gaps and hence all of the leaf area is within the crown surfa
 $f_g \to 1$, more of the leaf area is displaced deeper into the canopy, leaves in the
 lower crown intercepting light coming through holes in the upper canopy.
 
-```{code-cell}
+```{code-cell} ipython3
 fig, ax1 = plt.subplots(1, 1, figsize=(6, 5))
 
 for fg in np.arange(0, 1.01, 0.05):
@@ -600,7 +611,7 @@ ax1.legend(frameon=False)
 We can now calculate the crown area occupied by leaves above the height of each closed
 layer $z^*_l$:
 
-```{code-cell}
+```{code-cell} ipython3
 # Calculate the leaf area above z_star for each stem
 crown_gap_fraction = 0.05
 Acp_z_star = calculate_leaf_area(
@@ -614,7 +625,7 @@ And from that, the area occupied by leaves **within each layer**. These values a
 similar to the projected crown area within layers (`Ap_within_layer`, above) but
 leaf area is displaced into lower layers because $f_g > 0$.
 
-```{code-cell}
+```{code-cell} ipython3
 # Calculate the contribution _within_ each layer per stem
 Acp_within_layer = np.diff(Acp_z_star, axis=0, prepend=0)
 
@@ -631,7 +642,7 @@ $f_{abs} = 1 - e ^ {-kL}$,
 where $k$ is the light extinction coefficient ($k$) and $L$ is the leaf area index
 (LAI). The LAI can be calculated for each stem and layer:
 
-```{code-cell}
+```{code-cell} ipython3
 LAI = Acp_within_layer / canopy_area
 print(LAI)
 ```
@@ -639,7 +650,7 @@ print(LAI)
 This can be used to calculate the LAI of individual stems but also the LAI of each layer
 in the canopy:
 
-```{code-cell}
+```{code-cell} ipython3
 LAI_stem = LAI.sum(axis=0)
 LAI_layer = LAI.sum(axis=1)
 
@@ -650,7 +661,7 @@ print("LAI layer = ", LAI_layer)
 The layer LAI values can now be used to calculate the light transmission of each layer and
 hence the cumulative light extinction profile through the canopy.
 
-```{code-cell}
+```{code-cell} ipython3
 f_abs = 1 - np.exp(-pft.traits.par_ext * LAI_layer)
 ext = np.cumprod(f_abs)
 
@@ -662,7 +673,7 @@ One issue that needs to be resolved is that the T Model implementation in `pyrea
 follows the original implementation of the T Model in having LAI as a fixed trait of
 a given plant functional type, so is constant for all stems of that PFT.
 
-```{code-cell}
+```{code-cell} ipython3
 print("f_abs = ", (1 - np.exp(-pft.traits.par_ext * pft.traits.lai)))
 ```
 
diff --git a/docs/source/users/tmodel/tmodel.md b/docs/source/users/tmodel/tmodel.md
index fda39de0..d7e75b31 100644
--- a/docs/source/users/tmodel/tmodel.md
+++ b/docs/source/users/tmodel/tmodel.md
@@ -5,10 +5,21 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
+    jupytext_version: 1.16.4
 kernelspec:
   display_name: Python 3
   language: python
   name: python3
+language_info:
+  codemirror_mode:
+    name: ipython
+    version: 3
+  file_extension: .py
+  mimetype: text/x-python
+  name: python
+  nbconvert_exporter: python
+  pygments_lexer: ipython3
+  version: 3.11.9
 ---
 
 # The T Model
@@ -31,7 +42,7 @@ class description.
 
 The class can be used to create a default T Model trait set:
 
-```{code-cell}
+```{code-cell} ipython3
 import numpy as np
 from pyrealm import tmodel
 
@@ -42,7 +53,7 @@ print(traits1)
 
 It can also be edited to generate different growth patterns:
 
-```{code-cell}
+```{code-cell} ipython3
 # A slower growing tree with a higher maximum height
 traits2 = tmodel.TModelTraits(a_hd=50, h_max=40)
 print(traits2)
@@ -71,7 +82,7 @@ diameters and an optional set of traits as a
 {class}`~pyrealm.constants.tmodel_const.TModelTraits` object. If no traits are provided,
 the default {class}`~pyrealm.constants.tmodel_const.TModelTraits` settings are used.
 
-```{code-cell}
+```{code-cell} ipython3
 # Use a sequence of diameters from sapling to large tree
 diameters = np.linspace(0.02, 2, 100)
 tree1 = tmodel.TTree(diameters=diameters)  # Using default traits
@@ -92,7 +103,7 @@ These inputs are then immediately used to calculate the following properties of
 Using an array of diameter values provides an immediate way to visualise the geometric
 scaling resulting from a particular set of plant traits:
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 from matplotlib import pyplot
@@ -144,12 +155,12 @@ provide estimates of the following growth parameters:
 The code below calculates growth estimates at each diameter under a constant GPP of 7
 TODO - UNITS!.
 
-```{code-cell}
+```{code-cell} ipython3
 tree1.calculate_growth(np.array([7]))
 tree2.calculate_growth(np.array([7]))
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [hide-input]
 
 fig, (ax1, ax2, ax3) = pyplot.subplots(1, 3, figsize=(12, 4))
@@ -183,7 +194,7 @@ The {meth}`~pyrealm.tmodel.TTree.reset_diameters` can be used to update an exist
 {meth}`~pyrealm.tmodel.TTree.reset_diameters` automatically resets any calculated growth
 parameters: they will need to be recalculated for the new diameters.
 
-```{code-cell}
+```{code-cell} ipython3
 tree1.reset_diameters(np.array([0.0001]))
 print(tree1.height)
 ```
diff --git a/pyproject.toml b/pyproject.toml
index 2823f7cc..e436ab8a 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -31,15 +31,15 @@ repository = "https://github.com/ImperialCollegeLondon/pyrealm"
 version = "1.0.0"
 
 [tool.poetry.dependencies]
-dacite = "^1.6.0" 
-numpy = "^2.0.0" 
+dacite = "^1.6.0"
+numpy = "^2.0.0"
 python = ">=3.10"
-scipy = "^1.7.3" 
-tabulate = "^0.8.10" 
+scipy = "^1.7.3"
+tabulate = "^0.8.10"
 
 marshmallow = "^3.22.0"
-pandas = "^2.2.2"
 marshmallow-dataclass = "^8.7.0"
+pandas = "^2.2.2"
 pandas-stubs = "^2.2.2.240909"
 [tool.poetry.group.types.dependencies]
 pandas-stubs = "^2.2.0.240218"
@@ -132,7 +132,7 @@ select = [
   "I", # isort
   "UP", # pyupgrade
   "RUF", # RUF specific checks
-  "NPY201"
+  "NPY201",
 ]
 
 # On top of the Google convention, disable:
@@ -147,6 +147,6 @@ convention = "google"
 
 [tool.jupytext]
 # Stop jupytext from removing mystnb and other settings in MyST Notebook YAML headers
-notebook_metadata_filter = "-jupytext.text_representation.jupytext_version,settings,mystnb"
-# Also stop it from stripping cell metadata.
-cell_metadata_filter =  "all"
\ No newline at end of file
+notebook_metadata_filter = "settings,mystnb,language_info"
+# Also stop it from stripping cell metadata, except for specific ones to lose.
+cell_metadata_filter = "all,-trusted"