diff --git a/docs/source/conf.py b/docs/source/conf.py
index 2de60f86..2fe69fc6 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -99,6 +99,7 @@ class MyReferenceStyle(AuthorYearReferenceStyle):
     ("py:class", "numpy._typing._array_like._ScalarType_co"),
     ("py:class", "numpy._typing._generic_alias.ScalarType"),
     ("py:class", "numpy.float32"),
+    ("py:class", "numpy.float64"),
     ("py:class", "numpy.int64"),
     ("py:class", "numpy.timedelta64"),
     ("py:class", "numpy.bool_"),
diff --git a/pyrealm/constants/core_const.py b/pyrealm/constants/core_const.py
index 97f7408c..935d5cfb 100644
--- a/pyrealm/constants/core_const.py
+++ b/pyrealm/constants/core_const.py
@@ -114,7 +114,7 @@ class CoreConst(ConstantsClass):
     :cite:t:`berger:1978a`."""
 
     # Hygro constants
-    magnus_coef: NDArray[np.float32] = field(
+    magnus_coef: NDArray[np.float64] = field(
         default_factory=lambda: np.array((611.2, 17.62, 243.12))
     )
     """Three coefficients of the Magnus equation for saturated vapour pressure,
@@ -132,7 +132,7 @@ class CoreConst(ConstantsClass):
     """Set the method used for calculating water density ('fisher' or 'chen')."""
 
     # Fisher Dial
-    fisher_dial_lambda: NDArray[np.float32] = field(
+    fisher_dial_lambda: NDArray[np.float64] = field(
         default_factory=lambda: np.array(
             [1788.316, 21.55053, -0.4695911, 0.003096363, -7.341182e-06]
         )
@@ -140,7 +140,7 @@ class CoreConst(ConstantsClass):
     r"""Coefficients of the temperature dependent polynomial for :math:`\lambda`
      in the Tumlirz equation."""
 
-    fisher_dial_Po: NDArray[np.float32] = field(
+    fisher_dial_Po: NDArray[np.float64] = field(
         default_factory=lambda: np.array(
             [5918.499, 58.05267, -1.1253317, 0.0066123869, -1.4661625e-05]
         )
@@ -148,7 +148,7 @@ class CoreConst(ConstantsClass):
     """Coefficients of the temperature dependent polynomial for :math:`P_0` in the
     Tumlirz equation."""
 
-    fisher_dial_Vinf: NDArray[np.float32] = field(
+    fisher_dial_Vinf: NDArray[np.float64] = field(
         default_factory=lambda: np.array(
             [
                 0.6980547,
@@ -168,7 +168,7 @@ class CoreConst(ConstantsClass):
     in the Tumlirz equation."""
 
     # Chen water density
-    chen_po: NDArray[np.float32] = field(
+    chen_po: NDArray[np.float64] = field(
         default_factory=lambda: np.array(
             [
                 0.99983952,
@@ -186,7 +186,7 @@ class CoreConst(ConstantsClass):
     r"""Coefficients of the polynomial relationship of water density with temperature at
     1 atm (:math:`P^0`, kg/m^3) from :cite:t:`chen:2008a`."""
 
-    chen_ko: NDArray[np.float32] = field(
+    chen_ko: NDArray[np.float64] = field(
         default_factory=lambda: np.array(
             [19652.17, 148.1830, -2.29995, 0.01281, -4.91564e-5, 1.035530e-7]
         )
@@ -194,7 +194,7 @@ class CoreConst(ConstantsClass):
     r"""Polynomial relationship of bulk modulus of water with temperature at 1 atm
      (:math:`K^0`, kg/m^3) from :cite:t:`chen:2008a`."""
 
-    chen_ca: NDArray[np.float32] = field(
+    chen_ca: NDArray[np.float64] = field(
         default_factory=lambda: np.array(
             [3.26138, 5.223e-4, 1.324e-4, -7.655e-7, 8.584e-10]
         )
@@ -202,7 +202,7 @@ class CoreConst(ConstantsClass):
     r"""Coefficients of the polynomial temperature dependent coefficient :math:`A` from
      :cite:t:`chen:2008a`."""
 
-    chen_cb: NDArray[np.float32] = field(
+    chen_cb: NDArray[np.float64] = field(
         default_factory=lambda: np.array(
             [7.2061e-5, -5.8948e-6, 8.69900e-8, -1.0100e-9, 4.3220e-12]
         )
@@ -220,11 +220,11 @@ class CoreConst(ConstantsClass):
     huber_mu_ast: float = 1e-06
     r"""Huber reference pressure (:math:`\mu_{ast}` 1.0e-6, Pa s)"""
 
-    huber_H_i: NDArray[np.float32] = field(
+    huber_H_i: NDArray[np.float64] = field(
         default_factory=lambda: np.array([1.67752, 2.20462, 0.6366564, -0.241605])
     )
     """Temperature dependent parameterisation of Hi in Huber."""
-    huber_H_ij: NDArray[np.float32] = field(
+    huber_H_ij: NDArray[np.float64] = field(
         default_factory=lambda: np.array(
             [
                 [0.520094, 0.0850895, -1.08374, -0.289555, 0.0, 0.0],
diff --git a/pyrealm/constants/pmodel_const.py b/pyrealm/constants/pmodel_const.py
index 5e3746a6..fb5c7ebf 100644
--- a/pyrealm/constants/pmodel_const.py
+++ b/pyrealm/constants/pmodel_const.py
@@ -139,11 +139,11 @@ class PModelConst(ConstantsClass):
     # - note that kphio_C4 has been updated to account for an unintended double
     #   8 fold downscaling to account for the fraction of light reaching PS2.
     #   from original values of [-0.008, 0.00375, -0.58e-4]
-    kphio_C4: NDArray[np.float32] = field(
+    kphio_C4: NDArray[np.float64] = field(
         default_factory=lambda: np.array((-0.064, 0.03, -0.000464))
     )
     """Quadratic scaling of Kphio with temperature for C4 plants"""
-    kphio_C3: NDArray[np.float32] = field(
+    kphio_C3: NDArray[np.float64] = field(
         default_factory=lambda: np.array((0.352, 0.022, -0.00034))
     )
     """Quadratic scaling of Kphio with temperature for C3 plants"""
@@ -193,11 +193,11 @@ class PModelConst(ConstantsClass):
     """Exponent of the threshold function for Mengoli soil moisture"""
 
     # Unit cost ratio (beta) values for different CalcOptimalChi methods
-    beta_cost_ratio_prentice14: NDArray[np.float32] = field(
+    beta_cost_ratio_prentice14: NDArray[np.float64] = field(
         default_factory=lambda: np.array([146.0])
     )
     r"""Unit cost ratio for C3 plants (:math:`\beta`, 146.0)."""
-    beta_cost_ratio_c4: NDArray[np.float32] = field(
+    beta_cost_ratio_c4: NDArray[np.float64] = field(
         default_factory=lambda: np.array([146.0 / 9])
     )
     r"""Unit cost ratio for C4 plants (:math:`\beta`, 16.222)."""
diff --git a/pyrealm/core/hygro.py b/pyrealm/core/hygro.py
index e83880e6..fa48482d 100644
--- a/pyrealm/core/hygro.py
+++ b/pyrealm/core/hygro.py
@@ -12,7 +12,9 @@
 from pyrealm.core.utilities import bounds_checker, evaluate_horner_polynomial
 
 
-def calc_vp_sat(ta: NDArray, core_const: CoreConst = CoreConst()) -> NDArray:
+def calc_vp_sat(
+    ta: NDArray[np.float64], core_const: CoreConst = CoreConst()
+) -> NDArray[np.float64]:
     r"""Calculate vapour pressure of saturated air.
 
     This function calculates the vapour pressure of saturated air in kPa at a given
@@ -56,8 +58,10 @@ def calc_vp_sat(ta: NDArray, core_const: CoreConst = CoreConst()) -> NDArray:
 
 
 def convert_vp_to_vpd(
-    vp: NDArray, ta: NDArray, core_const: CoreConst = CoreConst()
-) -> NDArray:
+    vp: NDArray[np.float64],
+    ta: NDArray[np.float64],
+    core_const: CoreConst = CoreConst(),
+) -> NDArray[np.float64]:
     """Convert vapour pressure to vapour pressure deficit.
 
     Args:
@@ -86,8 +90,10 @@ def convert_vp_to_vpd(
 
 
 def convert_rh_to_vpd(
-    rh: NDArray, ta: NDArray, core_const: CoreConst = CoreConst()
-) -> NDArray:
+    rh: NDArray[np.float64],
+    ta: NDArray[np.float64],
+    core_const: CoreConst = CoreConst(),
+) -> NDArray[np.float64]:
     """Convert relative humidity to vapour pressure deficit.
 
     Args:
@@ -124,8 +130,10 @@ def convert_rh_to_vpd(
 
 
 def convert_sh_to_vp(
-    sh: NDArray, patm: NDArray, core_const: CoreConst = CoreConst()
-) -> NDArray:
+    sh: NDArray[np.float64],
+    patm: NDArray[np.float64],
+    core_const: CoreConst = CoreConst(),
+) -> NDArray[np.float64]:
     """Convert specific humidity to vapour pressure.
 
     Args:
@@ -149,8 +157,11 @@ def convert_sh_to_vp(
 
 
 def convert_sh_to_vpd(
-    sh: NDArray, ta: NDArray, patm: NDArray, core_const: CoreConst = CoreConst()
-) -> NDArray:
+    sh: NDArray[np.float64],
+    ta: NDArray[np.float64],
+    patm: NDArray[np.float64],
+    core_const: CoreConst = CoreConst(),
+) -> NDArray[np.float64]:
     """Convert specific humidity to vapour pressure deficit.
 
     Args:
@@ -185,7 +196,9 @@ def convert_sh_to_vpd(
 # The following functions are integrated from the evap.py implementation of SPLASH v1.
 
 
-def calc_saturation_vapour_pressure_slope(tc: NDArray) -> NDArray:
+def calc_saturation_vapour_pressure_slope(
+    tc: NDArray[np.float64],
+) -> NDArray[np.float64]:
     """Calculate the slope of the saturation vapour pressure curve.
 
     Calculates the slope of the saturation pressure temperature curve, following
@@ -207,7 +220,7 @@ def calc_saturation_vapour_pressure_slope(tc: NDArray) -> NDArray:
     )
 
 
-def calc_enthalpy_vaporisation(tc: NDArray) -> NDArray:
+def calc_enthalpy_vaporisation(tc: NDArray[np.float64]) -> NDArray[np.float64]:
     """Calculate the enthalpy of vaporization.
 
     Calculates the latent heat of vaporization of water as a function of
@@ -224,7 +237,7 @@ def calc_enthalpy_vaporisation(tc: NDArray) -> NDArray:
     return 1.91846e6 * ((tc + 273.15) / (tc + 273.15 - 33.91)) ** 2
 
 
-def calc_specific_heat(tc: NDArray) -> NDArray:
+def calc_specific_heat(tc: NDArray[np.float64]) -> NDArray[np.float64]:
     """Calculate the specific heat of air.
 
     Calculates the specific heat of air at a constant pressure (:math:`c_{pm}`, J/kg/K)
@@ -257,8 +270,8 @@ def calc_specific_heat(tc: NDArray) -> NDArray:
 
 
 def calc_psychrometric_constant(
-    tc: NDArray, p: NDArray, core_const: CoreConst = CoreConst()
-) -> NDArray:
+    tc: NDArray[np.float64], p: NDArray[np.float64], core_const: CoreConst = CoreConst()
+) -> NDArray[np.float64]:
     r"""Calculate the psychrometric constant.
 
     Calculates the psychrometric constant (:math:`\lambda`, Pa/K) given the temperature
diff --git a/pyrealm/core/pressure.py b/pyrealm/core/pressure.py
index 2a2e44d5..ed8334c5 100644
--- a/pyrealm/core/pressure.py
+++ b/pyrealm/core/pressure.py
@@ -2,12 +2,15 @@
 atmospheric pressure.
 """  # noqa D210, D415
 
+import numpy as np
 from numpy.typing import NDArray
 
 from pyrealm.constants import CoreConst
 
 
-def calc_patm(elv: NDArray, core_const: CoreConst = CoreConst()) -> NDArray:
+def calc_patm(
+    elv: NDArray[np.float64], core_const: CoreConst = CoreConst()
+) -> NDArray[np.float64]:
     r"""Calculate atmospheric pressure from elevation.
 
     Calculates atmospheric pressure as a function of elevation with reference to the
diff --git a/pyrealm/core/solar.py b/pyrealm/core/solar.py
index 6fc0db38..0f7ea2ec 100644
--- a/pyrealm/core/solar.py
+++ b/pyrealm/core/solar.py
@@ -24,7 +24,9 @@
 
 
 def calc_heliocentric_longitudes(
-    julian_day: NDArray, n_days: NDArray, core_const: CoreConst = CoreConst()
+    julian_day: NDArray[np.float64],
+    n_days: NDArray[np.float64],
+    core_const: CoreConst = CoreConst(),
 ) -> tuple[NDArray, NDArray]:
     """Calculate heliocentric longitude and anomaly.
 
diff --git a/pyrealm/core/utilities.py b/pyrealm/core/utilities.py
index 2e86d367..15ad9498 100644
--- a/pyrealm/core/utilities.py
+++ b/pyrealm/core/utilities.py
@@ -211,13 +211,13 @@ def _get_interval_functions(interval_type: str = "[]") -> tuple[np.ufunc, np.ufu
 
 
 def bounds_checker(
-    values: NDArray,
+    values: NDArray[np.float64],
     lower: float = -np.inf,
     upper: float = np.inf,
     interval_type: str = "[]",
     label: str = "",
     unit: str = "",
-) -> NDArray:
+) -> NDArray[np.float64]:
     r"""Check inputs fall within bounds.
 
     This is a simple pass through function that tests whether the values fall within
@@ -255,7 +255,7 @@ def bounds_checker(
 
 
 def bounds_mask(
-    inputs: NDArray,
+    inputs: NDArray[np.float64],
     lower: float = -np.inf,
     upper: float = np.inf,
     interval_type: str = "[]",
@@ -312,7 +312,7 @@ def bounds_mask(
     # modifying the original input. Using type
     if not np.issubdtype(inputs.dtype, np.floating):
         # Copies implicitly
-        outputs = inputs.astype(np.float32)
+        outputs = inputs.astype(np.float64)
     else:
         outputs = inputs.copy()
 
@@ -336,7 +336,9 @@ def bounds_mask(
     return outputs
 
 
-def evaluate_horner_polynomial(x: NDArray, cf: list | NDArray) -> NDArray:
+def evaluate_horner_polynomial(
+    x: NDArray[np.float64], cf: list | NDArray
+) -> NDArray[np.float64]:
     r"""Evaluates a polynomial with coefficients `cf` at `x` using Horner's method.
 
     Horner's method is a fast way to evaluate polynomials, especially for large degrees,
diff --git a/pyrealm/core/water.py b/pyrealm/core/water.py
index 72c7568b..953027c8 100644
--- a/pyrealm/core/water.py
+++ b/pyrealm/core/water.py
@@ -10,10 +10,10 @@
 
 
 def calc_density_h2o_chen(
-    tc: NDArray,
-    p: NDArray,
+    tc: NDArray[np.float64],
+    p: NDArray[np.float64],
     core_const: CoreConst = CoreConst(),
-) -> NDArray:
+) -> NDArray[np.float64]:
     """Calculate the density of water using Chen et al 2008.
 
     This function calculates the density of water at a given temperature and pressure
@@ -65,10 +65,10 @@ def calc_density_h2o_chen(
 
 
 def calc_density_h2o_fisher(
-    tc: NDArray,
-    patm: NDArray,
+    tc: NDArray[np.float64],
+    patm: NDArray[np.float64],
     core_const: CoreConst = CoreConst(),
-) -> NDArray:
+) -> NDArray[np.float64]:
     """Calculate water density.
 
     Calculates the density of water as a function of temperature and atmospheric
@@ -124,11 +124,11 @@ def calc_density_h2o_fisher(
 
 
 def calc_density_h2o(
-    tc: NDArray,
-    patm: NDArray,
+    tc: NDArray[np.float64],
+    patm: NDArray[np.float64],
     core_const: CoreConst = CoreConst(),
     safe: bool = True,
-) -> NDArray:
+) -> NDArray[np.float64]:
     """Calculate water density.
 
     Calculates the density of water as a function of temperature and atmospheric
@@ -179,11 +179,11 @@ def calc_density_h2o(
 
 
 def calc_viscosity_h2o(
-    tc: NDArray,
-    patm: NDArray,
+    tc: NDArray[np.float64],
+    patm: NDArray[np.float64],
     core_const: CoreConst = CoreConst(),
     simple: bool = False,
-) -> NDArray:
+) -> NDArray[np.float64]:
     r"""Calculate the viscosity of water.
 
     Calculates the viscosity of water (:math:`\eta`) as a function of temperature and
@@ -247,11 +247,11 @@ def calc_viscosity_h2o(
 
 
 def calc_viscosity_h2o_matrix(
-    tc: NDArray,
-    patm: NDArray,
+    tc: NDArray[np.float64],
+    patm: NDArray[np.float64],
     core_const: CoreConst = CoreConst(),
     simple: bool = False,
-) -> NDArray:
+) -> NDArray[np.float64]:
     r"""Calculate the viscosity of water.
 
     Calculates the viscosity of water (:math:`\eta`) as a function of temperature and
diff --git a/pyrealm/demography/canopy.py b/pyrealm/demography/canopy.py
index 49f47176..282b20d4 100644
--- a/pyrealm/demography/canopy.py
+++ b/pyrealm/demography/canopy.py
@@ -15,16 +15,16 @@
 
 def solve_canopy_area_filling_height(
     z: float,
-    stem_height: NDArray[np.float32],
-    crown_area: NDArray[np.float32],
-    m: NDArray[np.float32],
-    n: NDArray[np.float32],
-    q_m: NDArray[np.float32],
-    z_max: NDArray[np.float32],
-    n_individuals: NDArray[np.float32],
+    stem_height: NDArray[np.float64],
+    crown_area: NDArray[np.float64],
+    m: NDArray[np.float64],
+    n: NDArray[np.float64],
+    q_m: NDArray[np.float64],
+    z_max: NDArray[np.float64],
+    n_individuals: NDArray[np.float64],
     target_area: float = 0,
     validate: bool = True,
-) -> NDArray[np.float32]:
+) -> NDArray[np.float64]:
     """Solver function for finding the height where a canopy occupies a given area.
 
     This function takes the number of individuals in each cohort along with the stem
@@ -87,7 +87,7 @@ def fit_perfect_plasticity_approximation(
     canopy_gap_fraction: float,
     max_stem_height: float,
     solver_tolerance: float,
-) -> NDArray[np.float32]:
+) -> NDArray[np.float64]:
     r"""Find canopy layer heights under the PPA model.
 
     Finds the closure heights of the canopy layers under the perfect plasticity
@@ -126,7 +126,7 @@ def fit_perfect_plasticity_approximation(
     # Initialise the layer heights array and then loop over the layers indices,
     # except for the final layer, which will be the partial remaining vegetation below
     # the last closed layer.
-    layer_heights = np.zeros(n_layers, dtype=np.float32)
+    layer_heights = np.zeros(n_layers, dtype=np.float64)
     upper_bound = max_stem_height
 
     for layer in np.arange(n_layers - 1):
@@ -191,7 +191,7 @@ class Canopy:
     def __init__(
         self,
         community: Community,
-        layer_heights: NDArray[np.float32] | None = None,
+        layer_heights: NDArray[np.float64] | None = None,
         fit_ppa: bool = False,
         canopy_gap_fraction: float = 0,
         solver_tolerance: float = 0.001,
@@ -209,32 +209,32 @@ def __init__(
         """Total number of canopy layers."""
         self.n_cohorts: int
         """Total number of cohorts in the canopy."""
-        self.heights: NDArray[np.float32]
+        self.heights: NDArray[np.float64]
         """The vertical heights at which the canopy structure is calculated."""
 
         self.crown_profile: CrownProfile
         """The crown profiles of the community stems at the provided layer heights."""
-        self.stem_leaf_area: NDArray[np.float32]
+        self.stem_leaf_area: NDArray[np.float64]
         """The leaf area of the crown model for each cohort by layer."""
-        self.cohort_lai: NDArray[np.float32]
+        self.cohort_lai: NDArray[np.float64]
         """The leaf area index for each cohort by layer."""
-        self.cohort_f_trans: NDArray[np.float32]
+        self.cohort_f_trans: NDArray[np.float64]
         """The fraction of light transmitted by each cohort by layer."""
-        self.cohort_f_abs: NDArray[np.float32]
+        self.cohort_f_abs: NDArray[np.float64]
         """The fraction of light absorbed by each cohort by layer."""
-        self.f_trans: NDArray[np.float32]
+        self.f_trans: NDArray[np.float64]
         """The fraction of light transmitted by the whole community by layer."""
-        self.f_abs: NDArray[np.float32]
+        self.f_abs: NDArray[np.float64]
         """The fraction of light absorbed by the whole community by layer."""
-        self.transmission_profile: NDArray[np.float32]
+        self.transmission_profile: NDArray[np.float64]
         """The light transmission profile for the whole community by layer."""
-        self.extinction_profile: NDArray[np.float32]
+        self.extinction_profile: NDArray[np.float64]
         """The light extinction profile for the whole community by layer."""
-        self.fapar: NDArray[np.float32]
+        self.fapar: NDArray[np.float64]
         """The fraction of absorbed radiation for the whole community by layer."""
-        self.cohort_fapar: NDArray[np.float32]
+        self.cohort_fapar: NDArray[np.float64]
         """The fraction of absorbed radiation for each cohort by layer."""
-        self.stem_fapar: NDArray[np.float32]
+        self.stem_fapar: NDArray[np.float64]
         """The fraction of absorbed radiation for each stem by layer."""
         self.filled_community_area: float
         """The area filled by crown after accounting for the crown gap fraction."""
diff --git a/pyrealm/demography/community.py b/pyrealm/demography/community.py
index 11502f35..2f593471 100644
--- a/pyrealm/demography/community.py
+++ b/pyrealm/demography/community.py
@@ -353,7 +353,7 @@ class Community:
     flora: Flora
 
     # - arrays representing properties of cohorts
-    cohort_dbh_values: InitVar[NDArray[np.float32]]
+    cohort_dbh_values: InitVar[NDArray[np.float64]]
     cohort_n_individuals: InitVar[NDArray[np.int_]]
     cohort_pft_names: InitVar[NDArray[np.str_]]
 
@@ -365,7 +365,7 @@ class Community:
 
     def __post_init__(
         self,
-        cohort_dbh_values: NDArray[np.float32],
+        cohort_dbh_values: NDArray[np.float64],
         cohort_n_individuals: NDArray[np.int_],
         cohort_pft_names: NDArray[np.str_],
     ) -> None:
diff --git a/pyrealm/demography/crown.py b/pyrealm/demography/crown.py
index 90481507..5f0ac6ee 100644
--- a/pyrealm/demography/crown.py
+++ b/pyrealm/demography/crown.py
@@ -14,9 +14,9 @@
 
 
 def _validate_z_qz_args(
-    z: NDArray[np.float32],
-    stem_properties: list[NDArray[np.float32]],
-    q_z: NDArray[np.float32] | None = None,
+    z: NDArray[np.float64],
+    stem_properties: list[NDArray[np.float64]],
+    q_z: NDArray[np.float64] | None = None,
 ) -> None:
     """Shared validation of for crown function arguments.
 
@@ -95,12 +95,12 @@ def _validate_z_qz_args(
 
 
 def calculate_relative_crown_radius_at_z(
-    z: NDArray[np.float32],
-    stem_height: NDArray[np.float32],
-    m: NDArray[np.float32],
-    n: NDArray[np.float32],
+    z: NDArray[np.float64],
+    stem_height: NDArray[np.float64],
+    m: NDArray[np.float64],
+    n: NDArray[np.float64],
     validate: bool = True,
-) -> NDArray[np.float32]:
+) -> NDArray[np.float64]:
     r"""Calculate relative crown radius at a given height.
 
     The crown shape parameters ``m`` and ``n`` define the vertical distribution of
@@ -143,10 +143,10 @@ def calculate_relative_crown_radius_at_z(
 
 
 def calculate_crown_radius(
-    q_z: NDArray[np.float32],
-    r0: NDArray[np.float32],
+    q_z: NDArray[np.float64],
+    r0: NDArray[np.float64],
     validate: bool = True,
-) -> NDArray[np.float32]:
+) -> NDArray[np.float64]:
     r"""Calculate crown radius from relative crown radius and crown r0.
 
     The relative crown radius (:math:`q(z)`) at a given height :math:`z` describes the
@@ -175,14 +175,14 @@ def calculate_crown_radius(
 
 
 def calculate_stem_projected_crown_area_at_z(
-    z: NDArray[np.float32],
-    q_z: NDArray[np.float32],
-    stem_height: NDArray[np.float32],
-    crown_area: NDArray[np.float32],
-    q_m: NDArray[np.float32],
-    z_max: NDArray[np.float32],
+    z: NDArray[np.float64],
+    q_z: NDArray[np.float64],
+    stem_height: NDArray[np.float64],
+    crown_area: NDArray[np.float64],
+    q_m: NDArray[np.float64],
+    z_max: NDArray[np.float64],
     validate: bool = True,
-) -> NDArray[np.float32]:
+) -> NDArray[np.float64]:
     """Calculate stem projected crown area above a given height.
 
     This function calculates the projected crown area of a set of stems with given
@@ -221,15 +221,15 @@ def calculate_stem_projected_crown_area_at_z(
 
 
 def calculate_stem_projected_leaf_area_at_z(
-    z: NDArray[np.float32],
-    q_z: NDArray[np.float32],
-    stem_height: NDArray[np.float32],
-    crown_area: NDArray[np.float32],
-    f_g: NDArray[np.float32],
-    q_m: NDArray[np.float32],
-    z_max: NDArray[np.float32],
+    z: NDArray[np.float64],
+    q_z: NDArray[np.float64],
+    stem_height: NDArray[np.float64],
+    crown_area: NDArray[np.float64],
+    f_g: NDArray[np.float64],
+    q_m: NDArray[np.float64],
+    z_max: NDArray[np.float64],
     validate: bool = True,
-) -> NDArray[np.float32]:
+) -> NDArray[np.float64]:
     """Calculate projected leaf area above a given height.
 
     This function calculates the projected leaf area of a set of stems with given
@@ -328,16 +328,16 @@ class CrownProfile:
     """A Flora or StemTraits instance providing plant functional trait data."""
     stem_allometry: InitVar[StemAllometry]
     """A StemAllometry instance setting the stem allometries for the crown profile."""
-    z: NDArray[np.float32]
+    z: NDArray[np.float64]
     """An array of vertical height values at which to calculate crown profiles."""
 
-    relative_crown_radius: NDArray[np.float32] = field(init=False)
+    relative_crown_radius: NDArray[np.float64] = field(init=False)
     """An array of the relative crown radius of stems at z heights"""
-    crown_radius: NDArray[np.float32] = field(init=False)
+    crown_radius: NDArray[np.float64] = field(init=False)
     """An array of the actual crown radius of stems at z heights"""
-    projected_crown_area: NDArray[np.float32] = field(init=False)
+    projected_crown_area: NDArray[np.float64] = field(init=False)
     """An array of the projected crown area of stems at z heights"""
-    projected_leaf_area: NDArray[np.float32] = field(init=False)
+    projected_leaf_area: NDArray[np.float64] = field(init=False)
     """An array of the projected leaf area of stems at z heights"""
 
     # Information attributes
diff --git a/pyrealm/demography/flora.py b/pyrealm/demography/flora.py
index 41767d5f..a12463f9 100644
--- a/pyrealm/demography/flora.py
+++ b/pyrealm/demography/flora.py
@@ -46,8 +46,8 @@
 
 
 def calculate_crown_q_m(
-    m: float | NDArray[np.float32], n: float | NDArray[np.float32]
-) -> float | NDArray[np.float32]:
+    m: float | NDArray[np.float64], n: float | NDArray[np.float64]
+) -> float | NDArray[np.float64]:
     """Calculate the crown scaling trait ``q_m``.
 
     The value of q_m is a constant crown scaling parameter derived from the ``m`` and
@@ -66,8 +66,8 @@ def calculate_crown_q_m(
 
 
 def calculate_crown_z_max_proportion(
-    m: float | NDArray[np.float32], n: float | NDArray[np.float32]
-) -> float | NDArray[np.float32]:
+    m: float | NDArray[np.float64], n: float | NDArray[np.float64]
+) -> float | NDArray[np.float64]:
     r"""Calculate the z_m trait.
 
     The z_m proportion (:math:`p_{zm}`) is the constant proportion of stem height at
@@ -267,45 +267,45 @@ class Flora:
     # - trait arrays
     name: NDArray[np.str_] = field(init=False)
     r"""The name of the plant functional type."""
-    a_hd: NDArray[np.float32] = field(init=False)
+    a_hd: NDArray[np.float64] = field(init=False)
     r"""Initial slope of height-diameter relationship (:math:`a`, -)"""
-    ca_ratio: NDArray[np.float32] = field(init=False)
+    ca_ratio: NDArray[np.float64] = field(init=False)
     r"""Initial ratio of crown area to stem cross-sectional area
     (:math:`c`, -)"""
-    h_max: NDArray[np.float32] = field(init=False)
+    h_max: NDArray[np.float64] = field(init=False)
     r"""Maximum tree height (:math:`H_m`, m)"""
-    rho_s: NDArray[np.float32] = field(init=False)
+    rho_s: NDArray[np.float64] = field(init=False)
     r"""Sapwood density (:math:`\rho_s`, kg Cm-3)"""
-    lai: NDArray[np.float32] = field(init=False)
+    lai: NDArray[np.float64] = field(init=False)
     """Leaf area index within the crown (:math:`L`,  -)"""
-    sla: NDArray[np.float32] = field(init=False)
+    sla: NDArray[np.float64] = field(init=False)
     r"""Specific leaf area (:math:`\sigma`,  m2 kg-1 C)"""
-    tau_f: NDArray[np.float32] = field(init=False)
+    tau_f: NDArray[np.float64] = field(init=False)
     r"""Foliage turnover time (:math:`\tau_f`,years)"""
-    tau_r: NDArray[np.float32] = field(init=False)
+    tau_r: NDArray[np.float64] = field(init=False)
     r"""Fine-root turnover time (:math:`\tau_r`,  years)"""
-    par_ext: NDArray[np.float32] = field(init=False)
+    par_ext: NDArray[np.float64] = field(init=False)
     r"""Extinction coefficient of photosynthetically active radiation (PAR) (:math:`k`,
      -)"""
-    yld: NDArray[np.float32] = field(init=False)
+    yld: NDArray[np.float64] = field(init=False)
     r"""Yield factor (:math:`y`,  -)"""
-    zeta: NDArray[np.float32] = field(init=False)
+    zeta: NDArray[np.float64] = field(init=False)
     r"""Ratio of fine-root mass to foliage area (:math:`\zeta`, kg C m-2)"""
-    resp_r: NDArray[np.float32] = field(init=False)
+    resp_r: NDArray[np.float64] = field(init=False)
     r"""Fine-root specific respiration rate (:math:`r_r`, year-1)"""
-    resp_s: NDArray[np.float32] = field(init=False)
+    resp_s: NDArray[np.float64] = field(init=False)
     r"""Sapwood-specific respiration rate (:math:`r_s`,  year-1)"""
-    resp_f: NDArray[np.float32] = field(init=False)
+    resp_f: NDArray[np.float64] = field(init=False)
     r"""Foliage maintenance respiration fraction (:math:`r_f`,  -)"""
-    m: NDArray[np.float32] = field(init=False)
+    m: NDArray[np.float64] = field(init=False)
     r"""Crown shape parameter (:math:`m`, -)"""
-    n: NDArray[np.float32] = field(init=False)
+    n: NDArray[np.float64] = field(init=False)
     r"""Crown shape parameter (:math:`n`, -)"""
-    f_g: NDArray[np.float32] = field(init=False)
+    f_g: NDArray[np.float64] = field(init=False)
     r"""Crown gap fraction (:math:`f_g`, -)"""
-    q_m: NDArray[np.float32] = field(init=False)
+    q_m: NDArray[np.float64] = field(init=False)
     """Scaling factor to derive maximum crown radius from crown area."""
-    z_max_prop: NDArray[np.float32] = field(init=False)
+    z_max_prop: NDArray[np.float64] = field(init=False)
     """Proportion of stem height at which maximum crown radius is found."""
 
     # - other instance attributes
@@ -464,45 +464,45 @@ class StemTraits:
     # Instance trait attributes
     name: NDArray[np.str_]
     r"""The name of the plant functional type."""
-    a_hd: NDArray[np.float32]
+    a_hd: NDArray[np.float64]
     r"""Initial slope of height-diameter relationship (:math:`a`, -)"""
-    ca_ratio: NDArray[np.float32]
+    ca_ratio: NDArray[np.float64]
     r"""Initial ratio of crown area to stem cross-sectional area
     (:math:`c`, -)"""
-    h_max: NDArray[np.float32]
+    h_max: NDArray[np.float64]
     r"""Maximum tree height (:math:`H_m`, m)"""
-    rho_s: NDArray[np.float32]
+    rho_s: NDArray[np.float64]
     r"""Sapwood density (:math:`\rho_s`, kg Cm-3)"""
-    lai: NDArray[np.float32]
+    lai: NDArray[np.float64]
     """Leaf area index within the crown (:math:`L`,  -)"""
-    sla: NDArray[np.float32]
+    sla: NDArray[np.float64]
     r"""Specific leaf area (:math:`\sigma`,  m2 kg-1 C)"""
-    tau_f: NDArray[np.float32]
+    tau_f: NDArray[np.float64]
     r"""Foliage turnover time (:math:`\tau_f`,years)"""
-    tau_r: NDArray[np.float32]
+    tau_r: NDArray[np.float64]
     r"""Fine-root turnover time (:math:`\tau_r`,  years)"""
-    par_ext: NDArray[np.float32]
+    par_ext: NDArray[np.float64]
     r"""Extinction coefficient of photosynthetically active radiation (PAR) (:math:`k`,
      -)"""
-    yld: NDArray[np.float32]
+    yld: NDArray[np.float64]
     r"""Yield factor (:math:`y`,  -)"""
-    zeta: NDArray[np.float32]
+    zeta: NDArray[np.float64]
     r"""Ratio of fine-root mass to foliage area (:math:`\zeta`, kg C m-2)"""
-    resp_r: NDArray[np.float32]
+    resp_r: NDArray[np.float64]
     r"""Fine-root specific respiration rate (:math:`r_r`, year-1)"""
-    resp_s: NDArray[np.float32]
+    resp_s: NDArray[np.float64]
     r"""Sapwood-specific respiration rate (:math:`r_s`,  year-1)"""
-    resp_f: NDArray[np.float32]
+    resp_f: NDArray[np.float64]
     r"""Foliage maintenance respiration fraction (:math:`r_f`,  -)"""
-    m: NDArray[np.float32]
+    m: NDArray[np.float64]
     r"""Crown shape parameter (:math:`m`, -)"""
-    n: NDArray[np.float32]
+    n: NDArray[np.float64]
     r"""Crown shape parameter (:math:`n`, -)"""
-    f_g: NDArray[np.float32]
+    f_g: NDArray[np.float64]
     r"""Crown gap fraction (:math:`f_g`, -)"""
-    q_m: NDArray[np.float32]
+    q_m: NDArray[np.float64]
     """Scaling factor to derive maximum crown radius from crown area."""
-    z_max_prop: NDArray[np.float32]
+    z_max_prop: NDArray[np.float64]
     """Proportion of stem height at which maximum crown radius is found."""
 
     # Post init attributes
diff --git a/pyrealm/demography/t_model_functions.py b/pyrealm/demography/t_model_functions.py
index dbedb70f..8cabeaa4 100644
--- a/pyrealm/demography/t_model_functions.py
+++ b/pyrealm/demography/t_model_functions.py
@@ -53,11 +53,11 @@ def _validate_t_model_args(pft_args: list[NDArray], size_args: list[NDArray]) ->
 
 
 def calculate_heights(
-    h_max: NDArray[np.float32],
-    a_hd: NDArray[np.float32],
-    dbh: NDArray[np.float32],
+    h_max: NDArray[np.float64],
+    a_hd: NDArray[np.float64],
+    dbh: NDArray[np.float64],
     validate: bool = True,
-) -> NDArray[np.float32]:
+) -> NDArray[np.float64]:
     r"""Calculate tree height under the T Model.
 
     The height of trees (:math:`H`) are calculated from individual diameters at breast
@@ -83,11 +83,11 @@ def calculate_heights(
 
 
 def calculate_dbh_from_height(
-    h_max: NDArray[np.float32],
-    a_hd: NDArray[np.float32],
-    stem_height: NDArray[np.float32],
+    h_max: NDArray[np.float64],
+    a_hd: NDArray[np.float64],
+    stem_height: NDArray[np.float64],
     validate: bool = True,
-) -> NDArray[np.float32]:
+) -> NDArray[np.float64]:
     r"""Calculate diameter at breast height from stem height under the T Model.
 
     This function inverts the normal calculation of stem height (:math:`H`) from
@@ -133,12 +133,12 @@ def calculate_dbh_from_height(
 
 
 def calculate_crown_areas(
-    ca_ratio: NDArray[np.float32],
-    a_hd: NDArray[np.float32],
-    dbh: NDArray[np.float32],
-    stem_height: NDArray[np.float32],
+    ca_ratio: NDArray[np.float64],
+    a_hd: NDArray[np.float64],
+    dbh: NDArray[np.float64],
+    stem_height: NDArray[np.float64],
     validate: bool = True,
-) -> NDArray[np.float32]:
+) -> NDArray[np.float64]:
     r"""Calculate tree crown area under the T Model.
 
     The tree crown area (:math:`A_{c}`)is calculated from individual diameters at breast
@@ -166,11 +166,11 @@ def calculate_crown_areas(
 
 
 def calculate_crown_fractions(
-    a_hd: NDArray[np.float32],
-    stem_height: NDArray[np.float32],
-    dbh: NDArray[np.float32],
+    a_hd: NDArray[np.float64],
+    stem_height: NDArray[np.float64],
+    dbh: NDArray[np.float64],
     validate: bool = True,
-) -> NDArray[np.float32]:
+) -> NDArray[np.float64]:
     r"""Calculate tree crown fraction under the T Model.
 
     The crown fraction (:math:`f_{c}`)is calculated from individual diameters at breast
@@ -195,11 +195,11 @@ def calculate_crown_fractions(
 
 
 def calculate_stem_masses(
-    rho_s: NDArray[np.float32],
-    stem_height: NDArray[np.float32],
-    dbh: NDArray[np.float32],
+    rho_s: NDArray[np.float64],
+    stem_height: NDArray[np.float64],
+    dbh: NDArray[np.float64],
     validate: bool = True,
-) -> NDArray[np.float32]:
+) -> NDArray[np.float64]:
     r"""Calculate stem mass under the T Model.
 
     The stem mass (:math:`W_{s}`) is calculated from individual diameters at breast
@@ -223,11 +223,11 @@ def calculate_stem_masses(
 
 
 def calculate_foliage_masses(
-    sla: NDArray[np.float32],
-    lai: NDArray[np.float32],
-    crown_area: NDArray[np.float32],
+    sla: NDArray[np.float64],
+    lai: NDArray[np.float64],
+    crown_area: NDArray[np.float64],
     validate: bool = True,
-) -> NDArray[np.float32]:
+) -> NDArray[np.float64]:
     r"""Calculate foliage mass under the T Model.
 
     The foliage mass (:math:`W_{f}`) is calculated from the crown area (:math:`A_{c}`),
@@ -251,13 +251,13 @@ def calculate_foliage_masses(
 
 
 def calculate_sapwood_masses(
-    rho_s: NDArray[np.float32],
-    ca_ratio: NDArray[np.float32],
-    stem_height: NDArray[np.float32],
-    crown_area: NDArray[np.float32],
-    crown_fraction: NDArray[np.float32],
+    rho_s: NDArray[np.float64],
+    ca_ratio: NDArray[np.float64],
+    stem_height: NDArray[np.float64],
+    crown_area: NDArray[np.float64],
+    crown_fraction: NDArray[np.float64],
     validate: bool = True,
-) -> NDArray[np.float32]:
+) -> NDArray[np.float64]:
     r"""Calculate sapwood mass under the T Model.
 
     The sapwood mass (:math:`W_{\cdot s}`) is calculated from the individual crown area
@@ -287,8 +287,8 @@ def calculate_sapwood_masses(
 
 
 def calculate_crown_z_max(
-    z_max_prop: NDArray[np.float32], stem_height: NDArray[np.float32]
-) -> NDArray[np.float32]:
+    z_max_prop: NDArray[np.float64], stem_height: NDArray[np.float64]
+) -> NDArray[np.float64]:
     r"""Calculate height of maximum crown radius.
 
     The height of the maximum crown radius (:math:`z_m`) is derived from the crown
@@ -313,8 +313,8 @@ def calculate_crown_z_max(
 
 
 def calculate_crown_r0(
-    q_m: NDArray[np.float32], crown_area: NDArray[np.float32]
-) -> NDArray[np.float32]:
+    q_m: NDArray[np.float64], crown_area: NDArray[np.float64]
+) -> NDArray[np.float64]:
     r"""Calculate scaling factor for width of maximum crown radius.
 
     This scaling factor (:math:`r_0`) is derived from the crown shape parameters
@@ -339,12 +339,12 @@ def calculate_crown_r0(
 
 
 def calculate_whole_crown_gpp(
-    potential_gpp: NDArray[np.float32],
-    crown_area: NDArray[np.float32],
-    par_ext: NDArray[np.float32],
-    lai: NDArray[np.float32],
+    potential_gpp: NDArray[np.float64],
+    crown_area: NDArray[np.float64],
+    par_ext: NDArray[np.float64],
+    lai: NDArray[np.float64],
     validate: bool = True,
-) -> NDArray[np.float32]:
+) -> NDArray[np.float64]:
     r"""Calculate whole crown gross primary productivity.
 
     This function calculates individual GPP across the whole crown, given the individual
@@ -373,10 +373,10 @@ def calculate_whole_crown_gpp(
 
 
 def calculate_sapwood_respiration(
-    resp_s: NDArray[np.float32],
-    sapwood_mass: NDArray[np.float32],
+    resp_s: NDArray[np.float64],
+    sapwood_mass: NDArray[np.float64],
     validate: bool = True,
-) -> NDArray[np.float32]:
+) -> NDArray[np.float64]:
     r"""Calculate sapwood respiration.
 
     Calculates the total sapwood respiration (:math:`R_{\cdot s}`) given the individual
@@ -398,10 +398,10 @@ def calculate_sapwood_respiration(
 
 
 def calculate_foliar_respiration(
-    resp_f: NDArray[np.float32],
-    whole_crown_gpp: NDArray[np.float32],
+    resp_f: NDArray[np.float64],
+    whole_crown_gpp: NDArray[np.float64],
     validate: bool = True,
-) -> NDArray[np.float32]:
+) -> NDArray[np.float64]:
     r"""Calculate foliar respiration.
 
     Calculates the total foliar respiration (:math:`R_{f}`) given the individual crown
@@ -425,12 +425,12 @@ def calculate_foliar_respiration(
 
 
 def calculate_fine_root_respiration(
-    zeta: NDArray[np.float32],
-    sla: NDArray[np.float32],
-    resp_r: NDArray[np.float32],
-    foliage_mass: NDArray[np.float32],
+    zeta: NDArray[np.float64],
+    sla: NDArray[np.float64],
+    resp_r: NDArray[np.float64],
+    foliage_mass: NDArray[np.float64],
     validate: bool = True,
-) -> NDArray[np.float32]:
+) -> NDArray[np.float64]:
     r"""Calculate foliar respiration.
 
     Calculates the total fine root respiration (:math:`R_{r}`) given the individual
@@ -455,13 +455,13 @@ def calculate_fine_root_respiration(
 
 
 def calculate_net_primary_productivity(
-    yld: NDArray[np.float32],
-    whole_crown_gpp: NDArray[np.float32],
-    foliar_respiration: NDArray[np.float32],
-    fine_root_respiration: NDArray[np.float32],
-    sapwood_respiration: NDArray[np.float32],
+    yld: NDArray[np.float64],
+    whole_crown_gpp: NDArray[np.float64],
+    foliar_respiration: NDArray[np.float64],
+    fine_root_respiration: NDArray[np.float64],
+    sapwood_respiration: NDArray[np.float64],
     validate: bool = True,
-) -> NDArray[np.float32]:
+) -> NDArray[np.float64]:
     r"""Calculate net primary productivity.
 
     The net primary productivity (NPP, :math:`P_{net}`) is calculated as a plant
@@ -507,13 +507,13 @@ def calculate_net_primary_productivity(
 
 
 def calculate_foliage_and_fine_root_turnover(
-    sla: NDArray[np.float32],
-    zeta: NDArray[np.float32],
-    tau_f: NDArray[np.float32],
-    tau_r: NDArray[np.float32],
-    foliage_mass: NDArray[np.float32],
+    sla: NDArray[np.float64],
+    zeta: NDArray[np.float64],
+    tau_f: NDArray[np.float64],
+    tau_r: NDArray[np.float64],
+    foliage_mass: NDArray[np.float64],
     validate: bool = True,
-) -> NDArray[np.float32]:
+) -> NDArray[np.float64]:
     r"""Calculate turnover costs.
 
     This function calculates the costs associated with the turnover of fine roots and
@@ -544,19 +544,19 @@ def calculate_foliage_and_fine_root_turnover(
 
 
 def calculate_growth_increments(
-    rho_s: NDArray[np.float32],
-    a_hd: NDArray[np.float32],
-    h_max: NDArray[np.float32],
-    lai: NDArray[np.float32],
-    ca_ratio: NDArray[np.float32],
-    sla: NDArray[np.float32],
-    zeta: NDArray[np.float32],
-    npp: NDArray[np.float32],
-    turnover: NDArray[np.float32],
-    dbh: NDArray[np.float32],
-    stem_height: NDArray[np.float32],
+    rho_s: NDArray[np.float64],
+    a_hd: NDArray[np.float64],
+    h_max: NDArray[np.float64],
+    lai: NDArray[np.float64],
+    ca_ratio: NDArray[np.float64],
+    sla: NDArray[np.float64],
+    zeta: NDArray[np.float64],
+    npp: NDArray[np.float64],
+    turnover: NDArray[np.float64],
+    dbh: NDArray[np.float64],
+    stem_height: NDArray[np.float64],
     validate: bool = True,
-) -> tuple[NDArray[np.float32], NDArray[np.float32], NDArray[np.float32]]:
+) -> tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.float64]]:
     r"""Calculate growth increments.
 
     Given an estimate of net primary productivity (:math:`P_{net}`), less associated  
@@ -692,28 +692,28 @@ class StemAllometry:
     """ An instance of :class:`~pyrealm.demography.flora.Flora` or 
     :class:`~pyrealm.demography.flora.StemTraits`, providing plant functional trait data
     for a set of stems."""
-    at_dbh: InitVar[NDArray[np.float32]]
+    at_dbh: InitVar[NDArray[np.float64]]
     """An array of diameter at breast height values at which to predict stem allometry 
     values."""
 
     # Post init allometry attributes
-    dbh: NDArray[np.float32] = field(init=False)
+    dbh: NDArray[np.float64] = field(init=False)
     """Diameter at breast height (metres)"""
-    stem_height: NDArray[np.float32] = field(init=False)
+    stem_height: NDArray[np.float64] = field(init=False)
     """Stem height (metres)"""
-    crown_area: NDArray[np.float32] = field(init=False)
+    crown_area: NDArray[np.float64] = field(init=False)
     """Crown area (square metres)"""
-    crown_fraction: NDArray[np.float32] = field(init=False)
+    crown_fraction: NDArray[np.float64] = field(init=False)
     """Vertical fraction of the stem covered by the crown (-)"""
-    stem_mass: NDArray[np.float32] = field(init=False)
+    stem_mass: NDArray[np.float64] = field(init=False)
     """Stem mass (kg)"""
-    foliage_mass: NDArray[np.float32] = field(init=False)
+    foliage_mass: NDArray[np.float64] = field(init=False)
     """Foliage mass (kg)"""
-    sapwood_mass: NDArray[np.float32] = field(init=False)
+    sapwood_mass: NDArray[np.float64] = field(init=False)
     """Sapwood mass (kg)"""
-    crown_r0: NDArray[np.float32] = field(init=False)
+    crown_r0: NDArray[np.float64] = field(init=False)
     """Crown radius scaling factor (-)"""
-    crown_z_max: NDArray[np.float32] = field(init=False)
+    crown_z_max: NDArray[np.float64] = field(init=False)
     """Height of maximum crown radius (metres)"""
 
     # Information attributes
@@ -723,7 +723,7 @@ class StemAllometry:
     """The number of stems."""
 
     def __post_init__(
-        self, stem_traits: Flora | StemTraits, at_dbh: NDArray[np.float32]
+        self, stem_traits: Flora | StemTraits, at_dbh: NDArray[np.float64]
     ) -> None:
         """Populate the stem allometry attributes from the traits and size data."""
 
@@ -835,30 +835,30 @@ class StemAllocation:
     stem_allometry: InitVar[StemAllometry]
     """An instance of :class:`~pyrealm.demography.t_model_functions.StemAllometry`
     providing the stem size data for which to calculate allocation."""
-    at_potential_gpp: InitVar[NDArray[np.float32]]
+    at_potential_gpp: InitVar[NDArray[np.float64]]
     """An array of potential gross primary productivity for each stem that should be
     allocated to respiration, turnover and growth."""
 
     # Post init allometry attributes
-    potential_gpp: NDArray[np.float32] = field(init=False)
+    potential_gpp: NDArray[np.float64] = field(init=False)
     """Potential GPP per unit area (g C m2)"""
-    whole_crown_gpp: NDArray[np.float32] = field(init=False)
+    whole_crown_gpp: NDArray[np.float64] = field(init=False)
     """Estimated GPP across the whole crown (g C)"""
-    sapwood_respiration: NDArray[np.float32] = field(init=False)
+    sapwood_respiration: NDArray[np.float64] = field(init=False)
     """Allocation to sapwood respiration (g C)"""
-    foliar_respiration: NDArray[np.float32] = field(init=False)
+    foliar_respiration: NDArray[np.float64] = field(init=False)
     """Allocation to foliar respiration (g C)"""
-    fine_root_respiration: NDArray[np.float32] = field(init=False)
+    fine_root_respiration: NDArray[np.float64] = field(init=False)
     """Allocation to fine root respiration (g C)"""
-    npp: NDArray[np.float32] = field(init=False)
+    npp: NDArray[np.float64] = field(init=False)
     """Net primary productivity (g C)"""
-    turnover: NDArray[np.float32] = field(init=False)
+    turnover: NDArray[np.float64] = field(init=False)
     """Allocation to leaf and fine root turnover (g C)"""
-    delta_dbh: NDArray[np.float32] = field(init=False)
+    delta_dbh: NDArray[np.float64] = field(init=False)
     """Predicted increase in stem diameter from growth allocation (g C)"""
-    delta_stem_mass: NDArray[np.float32] = field(init=False)
+    delta_stem_mass: NDArray[np.float64] = field(init=False)
     """Predicted increase in stem mass from growth allocation (g C)"""
-    delta_foliage_mass: NDArray[np.float32] = field(init=False)
+    delta_foliage_mass: NDArray[np.float64] = field(init=False)
     """Predicted increase in foliar mass from growth allocation (g C)"""
 
     # Information attributes
@@ -871,7 +871,7 @@ def __post_init__(
         self,
         stem_traits: Flora | StemTraits,
         stem_allometry: StemAllometry,
-        at_potential_gpp: NDArray[np.float32],
+        at_potential_gpp: NDArray[np.float64],
     ) -> None:
         """Populate stem allocation attributes from the traits, allometry and GPP."""
 
diff --git a/pyrealm/pmodel/competition.py b/pyrealm/pmodel/competition.py
index 3be943b5..5c5755fb 100644
--- a/pyrealm/pmodel/competition.py
+++ b/pyrealm/pmodel/competition.py
@@ -12,8 +12,10 @@
 
 
 def convert_gpp_advantage_to_c4_fraction(
-    gpp_adv_c4: NDArray, treecover: NDArray, c3c4_const: C3C4Const = C3C4Const()
-) -> NDArray:
+    gpp_adv_c4: NDArray[np.float64],
+    treecover: NDArray[np.float64],
+    c3c4_const: C3C4Const = C3C4Const(),
+) -> NDArray[np.float64]:
     r"""Convert C4 GPP advantage to C4 fraction.
 
     This function calculates an initial estimate of the fraction of C4 plants based on
@@ -60,8 +62,8 @@ def convert_gpp_advantage_to_c4_fraction(
 
 
 def calculate_tree_proportion(
-    gppc3: NDArray, c3c4_const: C3C4Const = C3C4Const()
-) -> NDArray:
+    gppc3: NDArray[np.float64], c3c4_const: C3C4Const = C3C4Const()
+) -> NDArray[np.float64]:
     r"""Calculate the proportion of GPP from C3 trees.
 
     This function estimates the proportion of C3 trees in the community, which can then
@@ -181,11 +183,11 @@ class C3C4Competition:
 
     def __init__(
         self,
-        gpp_c3: NDArray,
-        gpp_c4: NDArray,
-        treecover: NDArray,
-        below_t_min: NDArray,
-        cropland: NDArray,
+        gpp_c3: NDArray[np.float64],
+        gpp_c4: NDArray[np.float64],
+        treecover: NDArray[np.float64],
+        below_t_min: NDArray[np.float64],
+        cropland: NDArray[np.float64],
         c3c4_const: C3C4Const = C3C4Const(),
     ):
         # Check inputs are congruent
@@ -198,7 +200,7 @@ def __init__(
         # annual total GPP estimates for C3 and C4 plants. This uses use
         # np.full to handle division by zero without raising warnings
         gpp_adv_c4 = np.full(self.shape, np.nan)
-        self.gpp_adv_c4: NDArray = np.divide(
+        self.gpp_adv_c4: NDArray[np.float64] = np.divide(
             gpp_c4 - gpp_c3, gpp_c3, out=gpp_adv_c4, where=gpp_c3 > 0
         )
         """The proportional advantage in GPP of C4 over C3 plants"""
@@ -224,22 +226,22 @@ def __init__(
         # Step 5: remove cropland areas
         frac_c4[cropland] = np.nan  # type: ignore
 
-        self.frac_c4: NDArray = frac_c4
+        self.frac_c4: NDArray[np.float64] = frac_c4
         """The estimated fraction of C4 plants."""
 
-        self.gpp_c3_contrib: NDArray = gpp_c3 * (1 - self.frac_c4)
+        self.gpp_c3_contrib: NDArray[np.float64] = gpp_c3 * (1 - self.frac_c4)
         """The estimated contribution of C3 plants to GPP (gC m-2 yr-1)"""
         self.gpp_c4_contrib = gpp_c4 * self.frac_c4
         """The estimated contribution of C4 plants to GPP (gC m-2 yr-1)"""
 
         # Define attributes used elsewhere
-        self.Delta13C_C3: NDArray
+        self.Delta13C_C3: NDArray[np.float64]
         r"""Contribution from C3 plants to (:math:`\Delta\ce{^13C}`, permil)."""
-        self.Delta13C_C4: NDArray
+        self.Delta13C_C4: NDArray[np.float64]
         r"""Contribution from C4 plants to (:math:`\Delta\ce{^13C}`, permil)."""
-        self.d13C_C3: NDArray
+        self.d13C_C3: NDArray[np.float64]
         r"""Contribution from C3 plants to (:math:`d\ce{^13C}`, permil)."""
-        self.d13C_C4: NDArray
+        self.d13C_C4: NDArray[np.float64]
         r"""Contribution from C3 plants to (:math:`d\ce{^13C}`, permil)."""
 
     def __repr__(self) -> str:
@@ -247,7 +249,10 @@ def __repr__(self) -> str:
         return f"C3C4Competition(shape={self.shape})"
 
     def estimate_isotopic_discrimination(
-        self, d13CO2: NDArray, Delta13C_C3_alone: NDArray, Delta13C_C4_alone: NDArray
+        self,
+        d13CO2: NDArray[np.float64],
+        Delta13C_C3_alone: NDArray[np.float64],
+        Delta13C_C4_alone: NDArray[np.float64],
     ) -> None:
         r"""Estimate CO2 isotopic discrimination values.
 
diff --git a/pyrealm/pmodel/functions.py b/pyrealm/pmodel/functions.py
index 69f4c98e..a9f3938a 100644
--- a/pyrealm/pmodel/functions.py
+++ b/pyrealm/pmodel/functions.py
@@ -12,12 +12,12 @@
 
 
 def calc_ftemp_arrh(
-    tk: NDArray,
+    tk: NDArray[np.float64],
     ha: float | NDArray,
     tk_ref: float | NDArray | None = None,
     pmodel_const: PModelConst = PModelConst(),
     core_const: CoreConst = CoreConst(),
-) -> NDArray:
+) -> NDArray[np.float64]:
     r"""Calculate enzyme kinetics scaling factor.
 
     Calculates the temperature-scaling factor :math:`f` for enzyme kinetics following an
@@ -81,9 +81,9 @@ def calc_ftemp_arrh(
 
 
 def calc_ftemp_inst_rd(
-    tc: NDArray,
+    tc: NDArray[np.float64],
     pmodel_const: PModelConst = PModelConst(),
-) -> NDArray:
+) -> NDArray[np.float64]:
     r"""Calculate temperature scaling of dark respiration.
 
     Calculates the temperature-scaling factor for dark respiration at a given
@@ -123,14 +123,14 @@ def calc_ftemp_inst_rd(
 
 
 def calc_modified_arrhenius_factor(
-    tk: NDArray,
+    tk: NDArray[np.float64],
     Ha: float | NDArray,
     Hd: float | NDArray,
     deltaS: float | NDArray,
     tk_ref: float | NDArray,
     mode: str = "M2002",
     core_const: CoreConst = CoreConst(),
-) -> NDArray:
+) -> NDArray[np.float64]:
     r"""Calculate the modified Arrhenius factor with temperature for an enzyme.
 
     This function returns a temperature-determined factor expressing the rate of an
@@ -198,8 +198,8 @@ def calc_modified_arrhenius_factor(
 
 
 def calc_ftemp_kphio(
-    tc: NDArray, c4: bool = False, pmodel_const: PModelConst = PModelConst()
-) -> NDArray:
+    tc: NDArray[np.float64], c4: bool = False, pmodel_const: PModelConst = PModelConst()
+) -> NDArray[np.float64]:
     r"""Calculate temperature dependence of quantum yield efficiency.
 
     Calculates the temperature dependence of the quantum yield efficiency, as a
@@ -257,11 +257,11 @@ def calc_ftemp_kphio(
 
 
 def calc_gammastar(
-    tc: NDArray,
-    patm: NDArray,
+    tc: NDArray[np.float64],
+    patm: NDArray[np.float64],
     pmodel_const: PModelConst = PModelConst(),
     core_const: CoreConst = CoreConst(),
-) -> NDArray:
+) -> NDArray[np.float64]:
     r"""Calculate the photorespiratory CO2 compensation point.
 
     Calculates the photorespiratory **CO2 compensation point** in absence of dark
@@ -311,10 +311,10 @@ def calc_gammastar(
 
 
 def calc_ns_star(
-    tc: NDArray,
-    patm: NDArray,
+    tc: NDArray[np.float64],
+    patm: NDArray[np.float64],
     core_const: CoreConst = CoreConst(),
-) -> NDArray:
+) -> NDArray[np.float64]:
     r"""Calculate the relative viscosity of water.
 
     Calculates the relative viscosity of water (:math:`\eta^*`), given the standard
@@ -355,11 +355,11 @@ def calc_ns_star(
 
 
 def calc_kmm(
-    tc: NDArray,
-    patm: NDArray,
+    tc: NDArray[np.float64],
+    patm: NDArray[np.float64],
     pmodel_const: PModelConst = PModelConst(),
     core_const: CoreConst = CoreConst(),
-) -> NDArray:
+) -> NDArray[np.float64]:
     r"""Calculate the Michaelis Menten coefficient of Rubisco-limited assimilation.
 
     Calculates the Michaelis Menten coefficient of Rubisco-limited assimilation
@@ -431,11 +431,11 @@ def calc_kmm(
 
 
 def calc_kp_c4(
-    tc: NDArray,
-    patm: NDArray,
+    tc: NDArray[np.float64],
+    patm: NDArray[np.float64],
     pmodel_const: PModelConst = PModelConst(),
     core_const: CoreConst = CoreConst(),
-) -> NDArray:
+) -> NDArray[np.float64]:
     r"""Calculate the Michaelis Menten coefficient of PEPc.
 
     Calculates the Michaelis Menten coefficient of phosphoenolpyruvate carboxylase
@@ -474,10 +474,10 @@ def calc_kp_c4(
 
 
 def calc_soilmstress_stocker(
-    soilm: NDArray,
-    meanalpha: NDArray = np.array(1.0),
+    soilm: NDArray[np.float64],
+    meanalpha: NDArray[np.float64] = np.array(1.0),
     pmodel_const: PModelConst = PModelConst(),
-) -> NDArray:
+) -> NDArray[np.float64]:
     r"""Calculate Stocker's empirical soil moisture stress factor.
 
     This function calculates a penalty factor :math:`\beta(\theta)` for well-watered GPP
@@ -562,10 +562,10 @@ def calc_soilmstress_stocker(
 
 
 def calc_soilmstress_mengoli(
-    soilm: NDArray = np.array(1.0),
-    aridity_index: NDArray = np.array(1.0),
+    soilm: NDArray[np.float64] = np.array(1.0),
+    aridity_index: NDArray[np.float64] = np.array(1.0),
     pmodel_const: PModelConst = PModelConst(),
-) -> NDArray:
+) -> NDArray[np.float64]:
     r"""Calculate the Mengoli et al. empirical soil moisture stress factor.
 
     This function calculates a penalty factor :math:`\beta(\theta)` for well-watered GPP
@@ -644,7 +644,9 @@ def calc_soilmstress_mengoli(
     return np.where(soilm >= psi, y, (y / psi) * soilm)
 
 
-def calc_co2_to_ca(co2: NDArray, patm: NDArray) -> NDArray:
+def calc_co2_to_ca(
+    co2: NDArray[np.float64], patm: NDArray[np.float64]
+) -> NDArray[np.float64]:
     r"""Convert :math:`\ce{CO2}` ppm to Pa.
 
     Converts ambient :math:`\ce{CO2}` (:math:`c_a`) in part per million to Pascals,
diff --git a/pyrealm/pmodel/isotopes.py b/pyrealm/pmodel/isotopes.py
index 50caeaab..ca7da2f4 100644
--- a/pyrealm/pmodel/isotopes.py
+++ b/pyrealm/pmodel/isotopes.py
@@ -5,6 +5,7 @@
 
 from warnings import warn
 
+import numpy as np
 from numpy.typing import NDArray
 
 from pyrealm.constants import IsotopesConst
@@ -42,8 +43,8 @@ class CalcCarbonIsotopes:
     def __init__(
         self,
         pmodel: PModel,
-        D14CO2: NDArray,
-        d13CO2: NDArray,
+        D14CO2: NDArray[np.float64],
+        d13CO2: NDArray[np.float64],
         isotopes_const: IsotopesConst = IsotopesConst(),
     ):
         # Check inputs are congruent
@@ -57,22 +58,22 @@ def __init__(
         """Indicates if estimates calculated for C3 or C4 photosynthesis."""
 
         # Attributes defined by methods below
-        self.Delta13C_simple: NDArray
+        self.Delta13C_simple: NDArray[np.float64]
         r"""Discrimination against carbon 13 (:math:`\Delta\ce{^{13}C}`, permil)
         excluding photorespiration."""
-        self.Delta14C: NDArray
+        self.Delta14C: NDArray[np.float64]
         r"""Discrimination against carbon 13 (:math:`\Delta\ce{^{13}C}`, permil)
         including photorespiration."""
-        self.Delta13C: NDArray
+        self.Delta13C: NDArray[np.float64]
         r"""Discrimination against carbon 14 (:math:`\Delta\ce{^{14}C}`, permil)
         including photorespiration."""
-        self.d13C_leaf: NDArray
+        self.d13C_leaf: NDArray[np.float64]
         r"""Isotopic ratio of carbon 13 in leaves
         (:math:`\delta\ce{^{13}C}`, permil)."""
-        self.d14C_leaf: NDArray
+        self.d14C_leaf: NDArray[np.float64]
         r"""Isotopic ratio of carbon 14 in leaves
         (:math:`\delta\ce{^{14}C}`, permil)."""
-        self.d13C_wood: NDArray
+        self.d13C_wood: NDArray[np.float64]
         r"""Isotopic ratio of carbon 13 in wood (:math:`\delta\ce{^{13}C}`, permil),
         given a parameterized post-photosynthetic fractionation."""
 
diff --git a/pyrealm/pmodel/jmax_limitation.py b/pyrealm/pmodel/jmax_limitation.py
index 27628f3f..cedc6793 100644
--- a/pyrealm/pmodel/jmax_limitation.py
+++ b/pyrealm/pmodel/jmax_limitation.py
@@ -89,14 +89,14 @@ def __init__(
 
         # Attributes populated by alternative method - two should always be populated by
         # the methods used below, but omega and omega_star only apply to smith19
-        self.f_j: NDArray
+        self.f_j: NDArray[np.float64]
         """:math:`J_{max}` limitation factor, calculated using the method."""
-        self.f_v: NDArray
+        self.f_v: NDArray[np.float64]
         """:math:`V_{cmax}` limitation factor, calculated using the method."""
-        self.omega: NDArray | None = None
+        self.omega: NDArray[np.float64] | None = None
         """Component of :math:`J_{max}` calculation for method ``smith19``
         (:cite:`Smith:2019dv`)."""
-        self.omega_star: NDArray | None = None
+        self.omega_star: NDArray[np.float64] | None = None
         """Component of :math:`J_{max}` calculation for method ``smith19``
         (:cite:`Smith:2019dv`)."""
 
diff --git a/pyrealm/pmodel/optimal_chi.py b/pyrealm/pmodel/optimal_chi.py
index 3b5b5ba7..86513b09 100644
--- a/pyrealm/pmodel/optimal_chi.py
+++ b/pyrealm/pmodel/optimal_chi.py
@@ -110,23 +110,23 @@ def __init__(
         # default value as they must be populated by the set_beta and estimate_chi
         # methods, which are called below, and so will be populated before __init__
         # returns.
-        self.beta: NDArray
+        self.beta: NDArray[np.float64]
         """The ratio of carboxylation to transpiration cost factors."""
-        self.xi: NDArray
+        self.xi: NDArray[np.float64]
         r"""Defines the sensitivity of :math:`\chi` to the vapour pressure deficit,
         related to the carbon cost of water (Medlyn et al. 2011; Prentice et 2014)."""
-        self.chi: NDArray
+        self.chi: NDArray[np.float64]
         r"""The ratio of leaf internal to ambient :math:`\ce{CO2}` partial pressure
         (:math:`\chi`)."""
-        self.mc: NDArray
+        self.mc: NDArray[np.float64]
         r""":math:`\ce{CO2}` limitation factor for RuBisCO-limited assimilation
         (:math:`m_c`)."""
-        self.mj: NDArray
+        self.mj: NDArray[np.float64]
         r""":math:`\ce{CO2}` limitation factor for light-limited assimilation
         (:math:`m_j`)."""
-        self.ci: NDArray
+        self.ci: NDArray[np.float64]
         r"""The leaf internal :math:`\ce{CO2}` partial pressure (:math:`c_i`)."""
-        self.mjoc: NDArray
+        self.mjoc: NDArray[np.float64]
         r"""Ratio of :math:`m_j/m_c`."""
 
         # Run the calculation methods after checking for any required variables
@@ -144,7 +144,7 @@ def set_beta(self) -> None:
         """Set the beta values."""
 
     @abstractmethod
-    def estimate_chi(self, xi_values: NDArray | None = None) -> None:
+    def estimate_chi(self, xi_values: NDArray[np.float64] | None = None) -> None:
         """Estimate xi, chi and other variables."""
 
     def _check_requires(self) -> None:
@@ -250,7 +250,7 @@ def set_beta(self) -> None:
         # leaf-internal-to-ambient CO2 partial pressure (ci/ca) ratio
         self.beta = self.pmodel_const.beta_cost_ratio_prentice14
 
-    def estimate_chi(self, xi_values: NDArray | None = None) -> None:
+    def estimate_chi(self, xi_values: NDArray[np.float64] | None = None) -> None:
         """Estimate ``chi`` for C3 plants."""
 
         if xi_values is not None:
@@ -319,7 +319,7 @@ def set_beta(self) -> None:
         # leaf-internal-to-ambient CO2 partial pressure (ci/ca) ratio
         self.beta = self.pmodel_const.beta_cost_ratio_prentice14
 
-    def estimate_chi(self, xi_values: NDArray | None = None) -> None:
+    def estimate_chi(self, xi_values: NDArray[np.float64] | None = None) -> None:
         """Estimate ``chi`` for C3 plants."""
 
         if xi_values is not None:
@@ -382,7 +382,7 @@ def set_beta(self) -> None:
         # leaf-internal-to-ambient CO2 partial pressure (ci/ca) ratio
         self.beta = self.pmodel_const.beta_cost_ratio_c4
 
-    def estimate_chi(self, xi_values: NDArray | None = None) -> None:
+    def estimate_chi(self, xi_values: NDArray[np.float64] | None = None) -> None:
         """Estimate ``chi`` for C4 plants, setting ``mj`` and ``mc`` to 1."""
         if xi_values is not None:
             _ = check_input_shapes(self.env.ca, xi_values)
@@ -449,7 +449,7 @@ def set_beta(self) -> None:
         # leaf-internal-to-ambient CO2 partial pressure (ci/ca) ratio
         self.beta = self.pmodel_const.beta_cost_ratio_c4
 
-    def estimate_chi(self, xi_values: NDArray | None = None) -> None:
+    def estimate_chi(self, xi_values: NDArray[np.float64] | None = None) -> None:
         """Estimate ``chi`` for C4 plants, setting ``mj`` and ``mc`` to 1."""
         if xi_values is not None:
             _ = check_input_shapes(self.env.ca, xi_values)
@@ -534,7 +534,7 @@ def set_beta(self) -> None:
             + self.pmodel_const.lavergne_2020_a_c3
         )
 
-    def estimate_chi(self, xi_values: NDArray | None = None) -> None:
+    def estimate_chi(self, xi_values: NDArray[np.float64] | None = None) -> None:
         """Estimate ``chi`` for C3 plants."""
 
         if xi_values is not None:
@@ -625,7 +625,7 @@ def set_beta(self) -> None:
             + self.pmodel_const.lavergne_2020_a_c4
         )
 
-    def estimate_chi(self, xi_values: NDArray | None = None) -> None:
+    def estimate_chi(self, xi_values: NDArray[np.float64] | None = None) -> None:
         """Estimate ``chi`` for C4 plants excluding photorespiration."""
 
         # Calculate chi and xi as in Prentice 14 but removing gamma terms.
@@ -702,7 +702,7 @@ def set_beta(self) -> None:
         # Calculate chi and xi as in Prentice 14 but removing gamma terms.
         self.beta = self.pmodel_const.beta_cost_ratio_c4
 
-    def estimate_chi(self, xi_values: NDArray | None = None) -> None:
+    def estimate_chi(self, xi_values: NDArray[np.float64] | None = None) -> None:
         """Estimate ``chi`` for C4 plants excluding photorespiration."""
 
         # Calculate chi and xi as in Prentice 14 but removing gamma terms.
@@ -771,7 +771,7 @@ def set_beta(self) -> None:
         # Calculate chi and xi as in Prentice 14 but removing gamma terms.
         self.beta = self.pmodel_const.beta_cost_ratio_c4
 
-    def estimate_chi(self, xi_values: NDArray | None = None) -> None:
+    def estimate_chi(self, xi_values: NDArray[np.float64] | None = None) -> None:
         """Estimate ``chi`` for C4 plants excluding photorespiration."""
 
         # Calculate chi and xi as in Prentice 14 but removing gamma terms.
diff --git a/pyrealm/pmodel/pmodel.py b/pyrealm/pmodel/pmodel.py
index c3ce20ca..abad0e80 100644
--- a/pyrealm/pmodel/pmodel.py
+++ b/pyrealm/pmodel/pmodel.py
@@ -238,14 +238,14 @@ def __init__(
         # in Pascals, but more commonly reported in µmol mol-1. The standard equation
         # (ca - ci) / 1.6 expects inputs in ppm, so the pascal versions are back
         # converted here.
-        self.iwue: NDArray = (5 / 8 * (env.ca - self.optchi.ci)) / (
+        self.iwue: NDArray[np.float64] = (5 / 8 * (env.ca - self.optchi.ci)) / (
             1e-6 * self.env.patm
         )
         """Intrinsic water use efficiency (iWUE, µmol mol-1)"""
 
         # The basic calculation of LUE = phi0 * M_c * m with an added penalty term
         # for jmax limitation
-        self.lue: NDArray = (
+        self.lue: NDArray[np.float64] = (
             self.kphio.kphio
             * self.optchi.mj
             * self.jmaxlim.f_v
@@ -258,14 +258,14 @@ def __init__(
         # no defaults and are only populated by estimate_productivity. Their getter
         # methods have a check to raise an informative error
         # -----------------------------------------------------------------------
-        self._vcmax: NDArray
-        self._vcmax25: NDArray
-        self._rd: NDArray
-        self._jmax: NDArray
-        self._gpp: NDArray
-        self._gs: NDArray
-        self._ppfd: NDArray
-        self._fapar: NDArray
+        self._vcmax: NDArray[np.float64]
+        self._vcmax25: NDArray[np.float64]
+        self._rd: NDArray[np.float64]
+        self._jmax: NDArray[np.float64]
+        self._gpp: NDArray[np.float64]
+        self._gs: NDArray[np.float64]
+        self._ppfd: NDArray[np.float64]
+        self._fapar: NDArray[np.float64]
 
     def _check_estimated(self, varname: str) -> None:
         """Raise error when accessing unpopulated parameters.
@@ -277,50 +277,50 @@ def _check_estimated(self, varname: str) -> None:
             raise RuntimeError(f"{varname} not calculated: use estimate_productivity")
 
     @property
-    def gpp(self) -> NDArray:
+    def gpp(self) -> NDArray[np.float64]:
         """Gross primary productivity (µg C m-2 s-1)."""
         self._check_estimated("gpp")
 
         return self._gpp
 
     @property
-    def vcmax(self) -> NDArray:
+    def vcmax(self) -> NDArray[np.float64]:
         """Maximum rate of carboxylation (µmol m-2 s-1)."""
         self._check_estimated("vcmax")
         return self._vcmax
 
     @property
-    def vcmax25(self) -> NDArray:
+    def vcmax25(self) -> NDArray[np.float64]:
         """Maximum rate of carboxylation at standard temperature (µmol m-2 s-1)."""
         self._check_estimated("vcmax25")
         return self._vcmax25
 
     @property
-    def rd(self) -> NDArray:
+    def rd(self) -> NDArray[np.float64]:
         """Dark respiration (µmol m-2 s-1)."""
         self._check_estimated("rd")
         return self._rd
 
     @property
-    def jmax(self) -> NDArray:
+    def jmax(self) -> NDArray[np.float64]:
         """Maximum rate of electron transport (µmol m-2 s-1)."""
         self._check_estimated("jmax")
         return self._jmax
 
     @property
-    def gs(self) -> NDArray:
+    def gs(self) -> NDArray[np.float64]:
         """Stomatal conductance (µmol m-2 s-1)."""
         self._check_estimated("gs")
         return self._gs
 
     @property
-    def ppfd(self) -> NDArray:
+    def ppfd(self) -> NDArray[np.float64]:
         """Photosynthetic photon flux density (PPFD, µmol m-2 s-1)."""
         self._check_estimated("gs")
         return self._ppfd
 
     @property
-    def fapar(self) -> NDArray:
+    def fapar(self) -> NDArray[np.float64]:
         """Fraction of absorbed photosynthetically active radiation
         (:math:`f_{APAR}` unitless).
         """  # noqa: D205
diff --git a/pyrealm/pmodel/pmodel_environment.py b/pyrealm/pmodel/pmodel_environment.py
index 34722c4c..a63083a7 100644
--- a/pyrealm/pmodel/pmodel_environment.py
+++ b/pyrealm/pmodel/pmodel_environment.py
@@ -75,27 +75,29 @@ class PModelEnvironment:
 
     def __init__(
         self,
-        tc: NDArray,
-        vpd: NDArray,
-        co2: NDArray,
-        patm: NDArray,
-        theta: NDArray | None = None,
-        rootzonestress: NDArray | None = None,
-        aridity_index: NDArray | None = None,
-        mean_growth_temperature: NDArray | None = None,
+        tc: NDArray[np.float64],
+        vpd: NDArray[np.float64],
+        co2: NDArray[np.float64],
+        patm: NDArray[np.float64],
+        theta: NDArray[np.float64] | None = None,
+        rootzonestress: NDArray[np.float64] | None = None,
+        aridity_index: NDArray[np.float64] | None = None,
+        mean_growth_temperature: NDArray[np.float64] | None = None,
         pmodel_const: PModelConst = PModelConst(),
         core_const: CoreConst = CoreConst(),
     ):
         self.shape: tuple = check_input_shapes(tc, vpd, co2, patm)
 
         # Validate and store the forcing variables
-        self.tc: NDArray = bounds_checker(tc, -25, 80, "[]", "tc", "°C")
+        self.tc: NDArray[np.float64] = bounds_checker(tc, -25, 80, "[]", "tc", "°C")
         """The temperature at which to estimate photosynthesis, °C"""
-        self.vpd: NDArray = bounds_checker(vpd, 0, 10000, "[]", "vpd", "Pa")
+        self.vpd: NDArray[np.float64] = bounds_checker(vpd, 0, 10000, "[]", "vpd", "Pa")
         """Vapour pressure deficit, Pa"""
-        self.co2: NDArray = bounds_checker(co2, 0, 1000, "[]", "co2", "ppm")
+        self.co2: NDArray[np.float64] = bounds_checker(co2, 0, 1000, "[]", "co2", "ppm")
         """CO2 concentration, ppm"""
-        self.patm: NDArray = bounds_checker(patm, 30000, 110000, "[]", "patm", "Pa")
+        self.patm: NDArray[np.float64] = bounds_checker(
+            patm, 30000, 110000, "[]", "patm", "Pa"
+        )
         """Atmospheric pressure, Pa"""
 
         # Guard against calc_density issues
@@ -112,7 +114,7 @@ def __init__(
                 "zero or explicitly set to np.nan"
             )
 
-        self.ca: NDArray = calc_co2_to_ca(self.co2, self.patm)
+        self.ca: NDArray[np.float64] = calc_co2_to_ca(self.co2, self.patm)
         """Ambient CO2 partial pressure, Pa"""
 
         self.gammastar = calc_gammastar(
@@ -140,13 +142,13 @@ def __init__(
         #        Easy to add the attributes dynamically, but bounds checking less
         #        obvious.
 
-        self.theta: NDArray
+        self.theta: NDArray[np.float64]
         """Volumetric soil moisture (m3/m3)"""
-        self.rootzonestress: NDArray
+        self.rootzonestress: NDArray[np.float64]
         """Rootzone stress factor (experimental) (-)"""
-        self.aridity_index: NDArray
+        self.aridity_index: NDArray[np.float64]
         """Climatological aridity index as PET/P (-)"""
-        self.mean_growth_temperature: NDArray
+        self.mean_growth_temperature: NDArray[np.float64]
         """Mean temperature > 0°C during growing degree days (°C)"""
 
         if theta is not None:
diff --git a/pyrealm/pmodel/quantum_yield.py b/pyrealm/pmodel/quantum_yield.py
index 63ece4be..9a7b0c8d 100644
--- a/pyrealm/pmodel/quantum_yield.py
+++ b/pyrealm/pmodel/quantum_yield.py
@@ -146,7 +146,7 @@ def __init__(
                     "of reference kphio values"
                 )
 
-        self.reference_kphio: NDArray = reference_kphio
+        self.reference_kphio: NDArray[np.float64] = reference_kphio
         """The kphio reference value for the method."""
         self.use_c4: bool = use_c4
         """Use a C4 parameterisation if available."""
@@ -154,7 +154,7 @@ def __init__(
         # Declare attributes populated by methods. These are typed but not assigned a
         # default value as they must are populated by the subclass specific
         # calculate_kphio method, which is called below to populate the values.
-        self.kphio: NDArray
+        self.kphio: NDArray[np.float64]
         """The calculated intrinsic quantum yield of photosynthesis."""
 
         # Run the calculation methods after checking for any required variables
@@ -300,7 +300,7 @@ class QuantumYieldSandoval(
     defaulting to the ratio of 1/9 in the absence of a Q cycle :cite:`long:1993a`.
     """
 
-    def peak_quantum_yield(self, aridity: NDArray) -> NDArray:
+    def peak_quantum_yield(self, aridity: NDArray[np.float64]) -> NDArray[np.float64]:
         """Calculate the peak quantum yield as a function of the aridity index.
 
         Args:
diff --git a/pyrealm/pmodel/scaler.py b/pyrealm/pmodel/scaler.py
index c7e9e7e2..a781a0e7 100644
--- a/pyrealm/pmodel/scaler.py
+++ b/pyrealm/pmodel/scaler.py
@@ -69,7 +69,7 @@ class SubdailyScaler:
 
     def __init__(
         self,
-        datetimes: NDArray,
+        datetimes: NDArray[np.datetime64],
     ) -> None:
         # Datetime validation. The inputs must be:
         # - one dimensional datetime64
@@ -306,7 +306,7 @@ def set_nearest(self, time: np.timedelta64) -> None:
 
     #     self.method = "Around max"
 
-    def pad_values(self, values: NDArray) -> NDArray:
+    def pad_values(self, values: NDArray[np.float64]) -> NDArray[np.float64]:
         """Pad values array to full days.
 
         This method takes an array representing daily values and pads the first and
@@ -329,7 +329,7 @@ def pad_values(self, values: NDArray) -> NDArray:
 
         return np.pad(values, padding_dims, constant_values=(np.nan, np.nan))
 
-    def get_window_values(self, values: NDArray) -> NDArray:
+    def get_window_values(self, values: NDArray[np.float64]) -> NDArray[np.float64]:
         """Extract acclimation window values for a variable.
 
         This method takes an array of values which has the same shape along the first
@@ -368,8 +368,8 @@ def get_window_values(self, values: NDArray) -> NDArray:
         return values_by_day[:, self.include, ...]
 
     def get_daily_means(
-        self, values: NDArray, allow_partial_data: bool = False
-    ) -> NDArray:
+        self, values: NDArray[np.float64], allow_partial_data: bool = False
+    ) -> NDArray[np.float64]:
         """Get the daily means of a variable during the acclimation window.
 
         This method extracts values from a given variable during a defined acclimation
@@ -406,12 +406,12 @@ def get_daily_means(
 
     def fill_daily_to_subdaily(
         self,
-        values: NDArray,
-        previous_value: NDArray | None = None,
+        values: NDArray[np.float64],
+        previous_value: NDArray[np.float64] | None = None,
         update_point: str = "max",
         kind: str = "previous",
         fill_from: np.timedelta64 | None = None,
-    ) -> NDArray:
+    ) -> NDArray[np.float64]:
         """Resample daily variables onto the subdaily time scale.
 
         This method takes an array representing daily values and interpolates those
diff --git a/pyrealm/pmodel/subdaily.py b/pyrealm/pmodel/subdaily.py
index 24d871d9..cd8277a7 100644
--- a/pyrealm/pmodel/subdaily.py
+++ b/pyrealm/pmodel/subdaily.py
@@ -41,11 +41,11 @@
 
 
 def memory_effect(
-    values: NDArray,
-    previous_values: NDArray | None = None,
+    values: NDArray[np.float64],
+    previous_values: NDArray[np.float64] | None = None,
     alpha: float = 0.067,
     allow_holdover: bool = False,
-) -> NDArray:
+) -> NDArray[np.float64]:
     r"""Apply a memory effect to a variable.
 
     Three key photosynthetic parameters (:math:`\xi`, :math:`V_{cmax25}` and
@@ -109,7 +109,7 @@ def memory_effect(
 
     # Initialise the output storage and set the first values to be a slice along the
     # first axis of the input values
-    memory_values = np.empty_like(values, dtype=np.float32)
+    memory_values = np.empty_like(values, dtype=np.float64)
     if previous_values is None:
         memory_values[0] = values[0]
     else:
@@ -245,8 +245,8 @@ def __init__(
         self,
         env: PModelEnvironment,
         fs_scaler: SubdailyScaler,
-        fapar: NDArray,
-        ppfd: NDArray,
+        fapar: NDArray[np.float64],
+        ppfd: NDArray[np.float64],
         method_optchi: str = "prentice14",
         method_jmaxlim: str = "wang17",
         method_kphio: str = "temperature",
@@ -427,21 +427,21 @@ def __init__(
             ]
 
         # 5) Calculate the realised daily values from the instantaneous optimal values
-        self.xi_real: NDArray = memory_effect(
+        self.xi_real: NDArray[np.float64] = memory_effect(
             self.pmodel_acclim.optchi.xi,
             previous_values=previous_xi_real,
             alpha=alpha,
             allow_holdover=allow_holdover,
         )
         r"""Realised daily slow responses in :math:`\xi`"""
-        self.vcmax25_real: NDArray = memory_effect(
+        self.vcmax25_real: NDArray[np.float64] = memory_effect(
             self.vcmax25_opt,
             previous_values=previous_vcmax25_real,
             alpha=alpha,
             allow_holdover=allow_holdover,
         )
         r"""Realised daily slow responses in :math:`V_{cmax25}`"""
-        self.jmax25_real: NDArray = memory_effect(
+        self.jmax25_real: NDArray[np.float64] = memory_effect(
             self.jmax25_opt,
             previous_values=previous_jmax25_real,
             alpha=alpha,
@@ -465,13 +465,19 @@ def __init__(
         # 7) Adjust subdaily jmax25 and vcmax25 back to jmax and vcmax given the
         #    actual subdaily temperatures.
         subdaily_tk = self.env.tc + self.env.core_const.k_CtoK
-        self.subdaily_vcmax: NDArray = self.subdaily_vcmax25 * calc_ftemp_arrh(
-            tk=subdaily_tk, ha=self.env.pmodel_const.subdaily_vcmax25_ha
+        self.subdaily_vcmax: NDArray[np.float64] = (
+            self.subdaily_vcmax25
+            * calc_ftemp_arrh(
+                tk=subdaily_tk, ha=self.env.pmodel_const.subdaily_vcmax25_ha
+            )
         )
         """Estimated subdaily :math:`V_{cmax}`."""
 
-        self.subdaily_jmax: NDArray = self.subdaily_jmax25 * calc_ftemp_arrh(
-            tk=subdaily_tk, ha=self.env.pmodel_const.subdaily_jmax25_ha
+        self.subdaily_jmax: NDArray[np.float64] = (
+            self.subdaily_jmax25
+            * calc_ftemp_arrh(
+                tk=subdaily_tk, ha=self.env.pmodel_const.subdaily_jmax25_ha
+            )
         )
         """Estimated subdaily :math:`J_{max}`."""
 
@@ -484,7 +490,9 @@ def __init__(
         """Estimated subdaily :math:`c_i`."""
 
         # Calculate Ac, J and Aj at subdaily scale to calculate assimilation
-        self.subdaily_Ac: NDArray = self.subdaily_vcmax * self.optimal_chi.mc
+        self.subdaily_Ac: NDArray[np.float64] = (
+            self.subdaily_vcmax * self.optimal_chi.mc
+        )
         """Estimated subdaily :math:`A_c`."""
 
         iabs = fapar * ppfd
@@ -493,11 +501,11 @@ def __init__(
             1 + ((4 * self.kphio.kphio * iabs) / self.subdaily_jmax) ** 2
         )
 
-        self.subdaily_Aj: NDArray = (subdaily_J / 4) * self.optimal_chi.mj
+        self.subdaily_Aj: NDArray[np.float64] = (subdaily_J / 4) * self.optimal_chi.mj
         """Estimated subdaily :math:`A_j`."""
 
         # Calculate GPP and convert from mol to gC
-        self.gpp: NDArray = (
+        self.gpp: NDArray[np.float64] = (
             np.minimum(self.subdaily_Aj, self.subdaily_Ac)
             * self.env.core_const.k_c_molmass
         )
@@ -598,8 +606,8 @@ def __init__(
         self,
         env: PModelEnvironment,
         fs_scaler: SubdailyScaler,
-        ppfd: NDArray,
-        fapar: NDArray,
+        ppfd: NDArray[np.float64],
+        fapar: NDArray[np.float64],
         alpha: float = 1 / 15,
         allow_holdover: bool = False,
         kphio: float = 1 / 8,
@@ -681,11 +689,11 @@ def __init__(
         )
 
         # Calculate the realised values from the instantaneous optimal values
-        self.vcmax25_real: NDArray = memory_effect(
+        self.vcmax25_real: NDArray[np.float64] = memory_effect(
             self.vcmax25_opt, alpha=alpha, allow_holdover=allow_holdover
         )
         r"""Realised daily slow responses in :math:`V_{cmax25}`"""
-        self.jmax25_real: NDArray = memory_effect(
+        self.jmax25_real: NDArray[np.float64] = memory_effect(
             self.jmax25_opt, alpha=alpha, allow_holdover=allow_holdover
         )
         r"""Realised daily slow responses in :math:`J_{max25}`"""
@@ -724,18 +732,24 @@ def __init__(
             self.jmax25_real, fill_from=fill_from
         )
 
-        self.subdaily_vcmax: NDArray = self.subdaily_vcmax25 * calc_ftemp_arrh(
-            tk=subdaily_tk, ha=self.env.pmodel_const.subdaily_vcmax25_ha
+        self.subdaily_vcmax: NDArray[np.float64] = (
+            self.subdaily_vcmax25
+            * calc_ftemp_arrh(
+                tk=subdaily_tk, ha=self.env.pmodel_const.subdaily_vcmax25_ha
+            )
         )
         """Estimated subdaily :math:`V_{cmax}`."""
 
-        self.subdaily_jmax: NDArray = self.subdaily_jmax25 * calc_ftemp_arrh(
-            tk=subdaily_tk, ha=self.env.pmodel_const.subdaily_jmax25_ha
+        self.subdaily_jmax: NDArray[np.float64] = (
+            self.subdaily_jmax25
+            * calc_ftemp_arrh(
+                tk=subdaily_tk, ha=self.env.pmodel_const.subdaily_jmax25_ha
+            )
         )
         """Estimated subdaily :math:`J_{max}`."""
 
         # Calculate Ac, J and Aj at subdaily scale to calculate assimilation
-        self.subdaily_Ac: NDArray = (
+        self.subdaily_Ac: NDArray[np.float64] = (
             self.subdaily_vcmax
             * (self.subdaily_ci - self.env.gammastar)
             / (self.subdaily_ci + self.env.kmm)
@@ -751,7 +765,7 @@ def __init__(
             1 + ((4 * self.kphio.kphio * iabs) / self.subdaily_jmax) ** 2
         )
 
-        self.subdaily_Aj: NDArray = (
+        self.subdaily_Aj: NDArray[np.float64] = (
             (subdaily_J / 4)
             * (self.subdaily_ci - self.env.gammastar)
             / (self.subdaily_ci + 2 * self.env.gammastar)
@@ -759,7 +773,7 @@ def __init__(
         """Estimated subdaily :math:`A_j`."""
 
         # Calculate GPP, converting from mol m2 s1 to grams carbon m2 s1
-        self.gpp: NDArray = (
+        self.gpp: NDArray[np.float64] = (
             np.minimum(self.subdaily_Aj, self.subdaily_Ac)
             * self.env.core_const.k_c_molmass
         )
diff --git a/pyrealm/splash/evap.py b/pyrealm/splash/evap.py
index 501f9eab..281720ff 100644
--- a/pyrealm/splash/evap.py
+++ b/pyrealm/splash/evap.py
@@ -49,29 +49,29 @@ class DailyEvapFluxes:
     solar: DailySolarFluxes
     pa: InitVar[NDArray]
     tc: InitVar[NDArray]
-    kWm: NDArray = field(default_factory=lambda: np.array([150.0]))
+    kWm: NDArray[np.float64] = field(default_factory=lambda: np.array([150.0]))
     core_const: CoreConst = field(default_factory=lambda: CoreConst())
 
-    sat: NDArray = field(init=False)
+    sat: NDArray[np.float64] = field(init=False)
     """Slope of saturation vapour pressure temperature curve, Pa/K"""
-    lv: NDArray = field(init=False)
+    lv: NDArray[np.float64] = field(init=False)
     """Enthalpy of vaporization, J/kg"""
-    pw: NDArray = field(init=False)
+    pw: NDArray[np.float64] = field(init=False)
     """Density of water, kg/m^3"""
-    psy: NDArray = field(init=False)
+    psy: NDArray[np.float64] = field(init=False)
     """Psychrometric constant, Pa/K"""
-    econ: NDArray = field(init=False)
+    econ: NDArray[np.float64] = field(init=False)
     """Water-to-energy conversion factor"""
-    cond: NDArray = field(init=False)
+    cond: NDArray[np.float64] = field(init=False)
     """Daily condensation, mm"""
-    eet_d: NDArray = field(init=False)
+    eet_d: NDArray[np.float64] = field(init=False)
     """Daily equilibrium evapotranspiration (EET), mm"""
-    pet_d: NDArray = field(init=False)
+    pet_d: NDArray[np.float64] = field(init=False)
     """Daily potential evapotranspiration (PET), mm"""
-    rx: NDArray = field(init=False)
+    rx: NDArray[np.float64] = field(init=False)
     """Variable substitute, (mm/hr)/(W/m^2)"""
 
-    def __post_init__(self, pa: NDArray, tc: NDArray) -> None:
+    def __post_init__(self, pa: NDArray[np.float64], tc: NDArray[np.float64]) -> None:
         """Calculate invariant components of evapotranspiration.
 
         The post_init method calculates the components of the evaporative fluxes that
@@ -107,8 +107,8 @@ def __post_init__(self, pa: NDArray, tc: NDArray) -> None:
         self.rx = (3.6e6) * (1.0 + self.core_const.k_w) * self.econ
 
     def estimate_aet(
-        self, wn: NDArray, day_idx: int | None = None, only_aet: bool = True
-    ) -> NDArray | tuple[NDArray, NDArray, NDArray]:
+        self, wn: NDArray[np.float64], day_idx: int | None = None, only_aet: bool = True
+    ) -> NDArray[np.float64] | tuple[NDArray, NDArray, NDArray]:
         """Estimate actual evapotranspiration.
 
         This method estimates the daily actual evapotranspiration (AET, mm/day), given
diff --git a/pyrealm/splash/solar.py b/pyrealm/splash/solar.py
index 3613fe0a..fc1e1831 100644
--- a/pyrealm/splash/solar.py
+++ b/pyrealm/splash/solar.py
@@ -37,39 +37,43 @@ class DailySolarFluxes:
     tc: InitVar[NDArray]
     core_const: CoreConst = field(default_factory=lambda: CoreConst())
 
-    nu: NDArray = field(init=False)
+    nu: NDArray[np.float64] = field(init=False)
     """True heliocentric anomaly, degrees"""
-    lambda_: NDArray = field(init=False)
+    lambda_: NDArray[np.float64] = field(init=False)
     """True heliocentric longitude, degrees"""
-    dr: NDArray = field(init=False)
+    dr: NDArray[np.float64] = field(init=False)
     """Distance factor, -"""
-    delta: NDArray = field(init=False)
+    delta: NDArray[np.float64] = field(init=False)
     """Declination angle, degrees"""
-    ru: NDArray = field(init=False)
+    ru: NDArray[np.float64] = field(init=False)
     """Intermediate variable, unitless"""
-    rv: NDArray = field(init=False)
+    rv: NDArray[np.float64] = field(init=False)
     """Intermediate variable, unitless"""
-    hs: NDArray = field(init=False)
+    hs: NDArray[np.float64] = field(init=False)
     """Sunset hour angle, degrees"""
-    ra_d: NDArray = field(init=False)
+    ra_d: NDArray[np.float64] = field(init=False)
     """Daily extraterrestrial solar radiation, J/m^2"""
-    tau: NDArray = field(init=False)
+    tau: NDArray[np.float64] = field(init=False)
     """Transmittivity, unitless"""
-    ppfd_d: NDArray = field(init=False)
+    ppfd_d: NDArray[np.float64] = field(init=False)
     """Daily PPFD, mol/m^2"""
-    rnl: NDArray = field(init=False)
+    rnl: NDArray[np.float64] = field(init=False)
     """Net longwave radiation, W/m^2"""
-    rw: NDArray = field(init=False)
+    rw: NDArray[np.float64] = field(init=False)
     """Intermediate variable,  W/m^2"""
-    hn: NDArray = field(init=False)
+    hn: NDArray[np.float64] = field(init=False)
     """Net radiation cross-over hour angle, degrees"""
-    rn_d: NDArray = field(init=False)
+    rn_d: NDArray[np.float64] = field(init=False)
     """Daytime net radiation, J/m^2"""
-    rnn_d: NDArray = field(init=False)
+    rnn_d: NDArray[np.float64] = field(init=False)
     """Nighttime net radiation (rnn_d), J/m^2"""
 
     def __post_init__(
-        self, lat: NDArray, elv: NDArray, sf: NDArray, tc: NDArray
+        self,
+        lat: NDArray[np.float64],
+        elv: NDArray[np.float64],
+        sf: NDArray[np.float64],
+        tc: NDArray[np.float64],
     ) -> None:
         """Populates key fluxes from input variables."""
 
diff --git a/pyrealm/splash/splash.py b/pyrealm/splash/splash.py
index d7116609..79c0230f 100644
--- a/pyrealm/splash/splash.py
+++ b/pyrealm/splash/splash.py
@@ -57,13 +57,13 @@ class SplashModel:
 
     def __init__(
         self,
-        lat: NDArray,
-        elv: NDArray,
-        sf: NDArray,
-        tc: NDArray,
-        pn: NDArray,
+        lat: NDArray[np.float64],
+        elv: NDArray[np.float64],
+        sf: NDArray[np.float64],
+        tc: NDArray[np.float64],
+        pn: NDArray[np.float64],
         dates: Calendar,
-        kWm: NDArray = np.array([150.0]),
+        kWm: NDArray[np.float64] = np.array([150.0]),
         core_const: CoreConst = CoreConst(),
     ):
         # Check input sizes are congurent
@@ -76,23 +76,31 @@ def __init__(
         if len(dates) != self.shape[0]:
             raise ValueError("Number of dates must match the first dimension of inputs")
 
-        self.elv: NDArray = elv
+        self.elv: NDArray[np.float64] = elv
         """The elevation of sites."""
-        self.lat: NDArray = bounds_checker(lat, -90, 90, label="lat", unit="°")
+        self.lat: NDArray[np.float64] = bounds_checker(
+            lat, -90, 90, label="lat", unit="°"
+        )
         """The latitude of sites."""
-        self.sf: NDArray = bounds_checker(sf, 0, 1, label="sf")
+        self.sf: NDArray[np.float64] = bounds_checker(sf, 0, 1, label="sf")
         """The sunshine fraction (0-1) of daily observations."""
-        self.tc: NDArray = bounds_checker(tc, -25, 80, label="tc", unit="°C")
+        self.tc: NDArray[np.float64] = bounds_checker(
+            tc, -25, 80, label="tc", unit="°C"
+        )
         """The air temperature in °C of daily observations."""
-        self.pn: NDArray = bounds_checker(pn, 0, 1e3, label="pn", unit="mm/day")
+        self.pn: NDArray[np.float64] = bounds_checker(
+            pn, 0, 1e3, label="pn", unit="mm/day"
+        )
         """The precipitation in mm of daily observations."""
         self.dates: Calendar = dates
         """The dates of observations along the first array axis."""
-        self.kWm: NDArray = bounds_checker(kWm, 0, 1e4, label="kWm", unit="mm")
+        self.kWm: NDArray[np.float64] = bounds_checker(
+            kWm, 0, 1e4, label="kWm", unit="mm"
+        )
         """The maximum soil water capacity for sites."""
 
         # TODO - potentially allow _actual_ climatic pressure data as an input
-        self.pa: NDArray = calc_patm(elv, core_const=core_const)
+        self.pa: NDArray[np.float64] = calc_patm(elv, core_const=core_const)
         """The atmospheric pressure at sites, derived from elevation"""
 
         # Calculate the daily solar fluxes - these are invariant across the simulation
@@ -109,12 +117,12 @@ def __init__(
 
     def estimate_initial_soil_moisture(  # noqa: max-complexity=12
         self,
-        wn_init: NDArray | None = None,
+        wn_init: NDArray[np.float64] | None = None,
         max_iter: int = 10,
         max_diff: float = 1.0,
         return_convergence: bool = False,
         verbose: bool = False,
-    ) -> NDArray:
+    ) -> NDArray[np.float64]:
         """Estimate initial soil moisture.
 
         This method uses the first year of data provided to a SplashModel instance to
@@ -222,7 +230,7 @@ def estimate_initial_soil_moisture(  # noqa: max-complexity=12
             return wn_start
 
     def estimate_daily_water_balance(
-        self, previous_wn: NDArray, day_idx: int | None = None
+        self, previous_wn: NDArray[np.float64], day_idx: int | None = None
     ) -> tuple[NDArray, NDArray, NDArray]:
         r"""Estimate the daily water balance.
 
@@ -286,7 +294,7 @@ def estimate_daily_water_balance(
 
     def calculate_soil_moisture(
         self,
-        wn_init: NDArray,
+        wn_init: NDArray[np.float64],
     ) -> tuple[NDArray, NDArray, NDArray]:
         """Calculate the soil moisture, AET and runoff from a SplashModel.
 
diff --git a/pyrealm/tmodel.py b/pyrealm/tmodel.py
index 35739d42..416c53c5 100644
--- a/pyrealm/tmodel.py
+++ b/pyrealm/tmodel.py
@@ -47,7 +47,7 @@ class TTree:
 
     def __init__(
         self,
-        diameters: NDArray,
+        diameters: NDArray[np.float64],
         traits: TModelTraits = TModelTraits(),
     ) -> None:
         self.traits: TModelTraits = traits
@@ -55,33 +55,33 @@ def __init__(
         # The diameter is used to define all of the geometric scaling
         # based on the trait parameters.
 
-        self._diameter: NDArray
-        self._height: NDArray
-        self._crown_fraction: NDArray
-        self._crown_area: NDArray
-        self._mass_stm: NDArray
-        self._mass_fol: NDArray
-        self._mass_swd: NDArray
+        self._diameter: NDArray[np.float64]
+        self._height: NDArray[np.float64]
+        self._crown_fraction: NDArray[np.float64]
+        self._crown_area: NDArray[np.float64]
+        self._mass_stm: NDArray[np.float64]
+        self._mass_fol: NDArray[np.float64]
+        self._mass_swd: NDArray[np.float64]
 
         self.reset_diameters(diameters)
 
         # Growth is then applied by providing estimated gpp using the
         # calculate_growth() method, which populates the following:
         self.growth_calculated: bool = False
-        self._gpp_raw: NDArray
-        self._gpp_actual: NDArray
-        self._npp: NDArray
-        self._resp_swd: NDArray
-        self._resp_frt: NDArray
-        self._resp_fol: NDArray
-        self._turnover: NDArray
-        self._d_mass_s: NDArray
-        self._d_mass_fr: NDArray
-        self._delta_d: NDArray
-        self._delta_mass_stm: NDArray
-        self._delta_mass_frt: NDArray
-
-    def _check_growth_calculated(self, property: str) -> NDArray:
+        self._gpp_raw: NDArray[np.float64]
+        self._gpp_actual: NDArray[np.float64]
+        self._npp: NDArray[np.float64]
+        self._resp_swd: NDArray[np.float64]
+        self._resp_frt: NDArray[np.float64]
+        self._resp_fol: NDArray[np.float64]
+        self._turnover: NDArray[np.float64]
+        self._d_mass_s: NDArray[np.float64]
+        self._d_mass_fr: NDArray[np.float64]
+        self._delta_d: NDArray[np.float64]
+        self._delta_mass_stm: NDArray[np.float64]
+        self._delta_mass_frt: NDArray[np.float64]
+
+    def _check_growth_calculated(self, property: str) -> NDArray[np.float64]:
         """Helper function to return growth values if calculated.
 
         This acts as a gatekeeper to make sure that a growth property is not returned
@@ -96,101 +96,101 @@ def _check_growth_calculated(self, property: str) -> NDArray:
         return getattr(self, property)
 
     @property
-    def diameter(self) -> NDArray:
+    def diameter(self) -> NDArray[np.float64]:
         """Individual diameter (m)."""
         return self._diameter
 
     @property
-    def height(self) -> NDArray:
+    def height(self) -> NDArray[np.float64]:
         """Individual height (m)."""
         return self._height
 
     @property
-    def crown_fraction(self) -> NDArray:
+    def crown_fraction(self) -> NDArray[np.float64]:
         """Individual crown fraction (unitless)."""
         return self._crown_fraction
 
     @property
-    def crown_area(self) -> NDArray:
+    def crown_area(self) -> NDArray[np.float64]:
         """Individual crown area (m2)."""
         return self._crown_area
 
     @property
-    def mass_swd(self) -> NDArray:
+    def mass_swd(self) -> NDArray[np.float64]:
         """Individual softwood mass (kg)."""
         return self._mass_swd
 
     @property
-    def mass_stm(self) -> NDArray:
+    def mass_stm(self) -> NDArray[np.float64]:
         """Individual stem mass (kg)."""
         return self._mass_stm
 
     @property
-    def mass_fol(self) -> NDArray:
+    def mass_fol(self) -> NDArray[np.float64]:
         """Individual foliage mass (kg)."""
         return self._mass_fol
 
     @property
-    def gpp_raw(self) -> NDArray:
+    def gpp_raw(self) -> NDArray[np.float64]:
         """Raw gross primary productivity."""
         return self._check_growth_calculated("_gpp_raw")
 
     @property
-    def gpp_actual(self) -> NDArray:
+    def gpp_actual(self) -> NDArray[np.float64]:
         """Actual gross primary productivity."""
         return self._check_growth_calculated("_gpp_actual")
 
     @property
-    def resp_swd(self) -> NDArray:
+    def resp_swd(self) -> NDArray[np.float64]:
         """Individual softwood respiration costs."""
         return self._check_growth_calculated("_resp_swd")
 
     @property
-    def resp_frt(self) -> NDArray:
+    def resp_frt(self) -> NDArray[np.float64]:
         """Individual fine root respiration costs."""
         return self._check_growth_calculated("_resp_frt")
 
     @property
-    def resp_fol(self) -> NDArray:
+    def resp_fol(self) -> NDArray[np.float64]:
         """Individual foliar respiration costs."""
         return self._check_growth_calculated("_resp_fol")
 
     @property
-    def npp(self) -> NDArray:
+    def npp(self) -> NDArray[np.float64]:
         """Net primary productivity."""
         return self._check_growth_calculated("_npp")
 
     @property
-    def turnover(self) -> NDArray:
+    def turnover(self) -> NDArray[np.float64]:
         """Plant turnover."""
         return self._check_growth_calculated("_turnover")
 
     @property
-    def d_mass_s(self) -> NDArray:
+    def d_mass_s(self) -> NDArray[np.float64]:
         """Individual relative change in mass."""
         return self._check_growth_calculated("_d_mass_s")
 
     @property
-    def d_mass_fr(self) -> NDArray:
+    def d_mass_fr(self) -> NDArray[np.float64]:
         """Individual relative change in fine root mass."""
         return self._check_growth_calculated("_d_mass_fr")
 
     @property
-    def delta_d(self) -> NDArray:
+    def delta_d(self) -> NDArray[np.float64]:
         """Individual change in diameter."""
         return self._check_growth_calculated("_delta_d")
 
     @property
-    def delta_mass_stm(self) -> NDArray:
+    def delta_mass_stm(self) -> NDArray[np.float64]:
         """Individual total change in stem mass."""
         return self._check_growth_calculated("_delta_mass_stm")
 
     @property
-    def delta_mass_frt(self) -> NDArray:
+    def delta_mass_frt(self) -> NDArray[np.float64]:
         """Individual total change in fine root mass."""
         return self._check_growth_calculated("_delta_mass_frt")
 
-    def reset_diameters(self, values: NDArray) -> None:
+    def reset_diameters(self, values: NDArray[np.float64]) -> None:
         """Reset the stem diameters for the T model.
 
         The set_diameter method can be used to reset the diameter values and then uses
@@ -239,7 +239,7 @@ def reset_diameters(self, values: NDArray) -> None:
         # Flag any calculated growth values as outdated
         self.growth_calculated = False
 
-    def calculate_growth(self, gpp: NDArray) -> None:
+    def calculate_growth(self, gpp: NDArray[np.float64]) -> None:
         """Calculate growth predictions given a GPP estimate.
 
         This method updates the instance with predicted changes in tree geometry, mass
@@ -318,8 +318,8 @@ def calculate_growth(self, gpp: NDArray) -> None:
 
 
 def grow_ttree(
-    gpp: NDArray,
-    d_init: NDArray,
+    gpp: NDArray[np.float64],
+    d_init: NDArray[np.float64],
     time_axis: int,
     traits: TModelTraits = TModelTraits(),
     outvars: tuple[str, ...] = ("diameter", "height", "crown_area", "delta_d"),