diff --git a/pyrealm/constants/pmodel_const.py b/pyrealm/constants/pmodel_const.py
index 6e418569..31bcf5a9 100644
--- a/pyrealm/constants/pmodel_const.py
+++ b/pyrealm/constants/pmodel_const.py
@@ -170,34 +170,42 @@ class PModelConst(ConstantsClass):
     (:math:`V_{\infty}`)."""
 
     # Chen water density
-    chen_po: NDArray[np.float32] = np.array(
-        [
-            0.99983952,
-            6.788260e-5,
-            -9.08659e-6,
-            1.022130e-7,
-            -1.35439e-9,
-            1.471150e-11,
-            -1.11663e-13,
-            5.044070e-16,
-            -1.00659e-18,
-        ]
+    chen_po: NDArray[np.float32] = field(
+        default_factory=lambda: np.array(
+            [
+                0.99983952,
+                6.788260e-5,
+                -9.08659e-6,
+                1.022130e-7,
+                -1.35439e-9,
+                1.471150e-11,
+                -1.11663e-13,
+                5.044070e-16,
+                -1.00659e-18,
+            ]
+        )
     )
     r"""Polynomial relationship of water density at 1 atm (kg/m^3) with temperature."""
 
-    chen_ko: NDArray[np.float32] = np.array(
-        [19652.17, 148.1830, -2.29995, 0.01281, -4.91564e-5, 1.035530e-7]
+    chen_ko: NDArray[np.float32] = field(
+        default_factory=lambda: np.array(
+            [19652.17, 148.1830, -2.29995, 0.01281, -4.91564e-5, 1.035530e-7]
+        )
     )
     r"""Polynomial relationship of bulk modulus of water at 1 atm (kg/m^3) with
      temperature."""
 
-    chen_ca: NDArray[np.float32] = np.array(
-        [3.26138, 5.223e-4, 1.324e-4, -7.655e-7, 8.584e-10]
+    chen_ca: NDArray[np.float32] = field(
+        default_factory=lambda: np.array(
+            [3.26138, 5.223e-4, 1.324e-4, -7.655e-7, 8.584e-10]
+        )
     )
     r"""Polynomial temperature dependent coefficient :math:`c_{a}`."""
 
-    chen_cb: NDArray[np.float32] = np.array(
-        [7.2061e-5, -5.8948e-6, 8.69900e-8, -1.0100e-9, 4.3220e-12]
+    chen_cb: NDArray[np.float32] = field(
+        default_factory=lambda: np.array(
+            [7.2061e-5, -5.8948e-6, 8.69900e-8, -1.0100e-9, 4.3220e-12]
+        )
     )
     r"""Polynomial temperature dependent coefficient :math:`c_{b}`."""
 
diff --git a/tests/pmodel/test_functions.py b/tests/pmodel/test_functions.py
index 2e2f23cd..9e7af3c6 100644
--- a/tests/pmodel/test_functions.py
+++ b/tests/pmodel/test_functions.py
@@ -1,4 +1,4 @@
-"""Test the pmodel functions
+"""Test the pmodel functions.
 
 TODO - note that there are parallel tests in test_pmodel that benchmark against the
 rpmodel outputs and test a wider range of inputs. Those could be moved here. These tests
@@ -10,7 +10,7 @@
 
 @pytest.mark.parametrize(argnames="shape", argvalues=[(1,), (6, 9), (4, 7, 3)])
 def test_calc_density_h20_fisher(shape):
-    """Test the fisher method"""
+    """Test the fisher method."""
     from pyrealm.pmodel.functions import calc_density_h2o_fisher
 
     rho = calc_density_h2o_fisher(
@@ -22,7 +22,7 @@ def test_calc_density_h20_fisher(shape):
 
 @pytest.mark.parametrize(argnames="shape", argvalues=[(1,), (6, 9), (4, 7, 3)])
 def test_calc_density_h20_chen(shape):
-    """Test the chen method"""
+    """Test the chen method."""
     from pyrealm.pmodel.functions import calc_density_h2o_chen
 
     rho = calc_density_h2o_chen(
@@ -41,7 +41,7 @@ def test_calc_density_h20_chen(shape):
     ],
 )
 def test_calc_density_h20(const_args, exp):
-    """Test the wrapper method dispatches as expected"""
+    """Test the wrapper method dispatches as expected."""
     from pyrealm.constants import PModelConst
     from pyrealm.pmodel.functions import calc_density_h2o
 
@@ -57,7 +57,7 @@ def test_calc_density_h20(const_args, exp):
 
 @pytest.mark.parametrize(argnames="shape", argvalues=[(1,), (6, 9), (4, 7, 3)])
 def test_calc_viscosity_h20(shape):
-    """Test the viscosity calculation"""
+    """Test the viscosity calculation."""
     from pyrealm.pmodel.functions import calc_viscosity_h2o
 
     eta = calc_viscosity_h2o(
@@ -69,7 +69,7 @@ def test_calc_viscosity_h20(shape):
 
 @pytest.mark.parametrize(argnames="shape", argvalues=[(1,), (6, 9), (4, 7, 3)])
 def test_calc_viscosity_h20_matrix(shape):
-    """Test the viscosity calculation"""
+    """Test the viscosity calculation."""
     from pyrealm.pmodel.functions import calc_viscosity_h2o_matrix
 
     eta = calc_viscosity_h2o_matrix(