diff --git a/docs/source/_static/ve_full_model_configuration.toml b/docs/source/_static/ve_full_model_configuration.toml
index c110584f0..27e5027f9 100644
--- a/docs/source/_static/ve_full_model_configuration.toml
+++ b/docs/source/_static/ve_full_model_configuration.toml
@@ -49,10 +49,6 @@ var_name = "surface_runoff"
 file_path = "/private/tmp/ve_example/data/example_soil_data.nc"
 var_name = "pH"
 
-[[core.data.variable]]
-file_path = "/private/tmp/ve_example/data/example_soil_data.nc"
-var_name = "bulk_density"
-
 [[core.data.variable]]
 file_path = "/private/tmp/ve_example/data/example_soil_data.nc"
 var_name = "clay_fraction"
diff --git a/docs/source/using_the_ve/example_data.md b/docs/source/using_the_ve/example_data.md
index 04428a39c..28c5e77f2 100644
--- a/docs/source/using_the_ve/example_data.md
+++ b/docs/source/using_the_ve/example_data.md
@@ -241,10 +241,6 @@ The `example_soil_data.nc` file provides:
   - `pH`
   - unitless
   - XY
-* - Bulk soil density
-  - `bulk_density`
-  - kg $\textrm{m}^{-3}$
-  - XY
 * - Soil clay fraction
   - `clay_fraction`
   - unitless
diff --git a/tests/models/soil/conftest.py b/tests/models/soil/conftest.py
index 934dd546a..7a4bfe2a7 100644
--- a/tests/models/soil/conftest.py
+++ b/tests/models/soil/conftest.py
@@ -46,7 +46,6 @@ def dummy_carbon_data(fixture_core_components):
         "soil_p_pool_secondary": [0.00705668, 0.03816896, 0.01152589, 0.00733107],
         "soil_p_pool_labile": [1.0582393e-5, 3.252961e-5, 6.806745e-5, 1.945635e-4],
         "pH": [3.0, 7.5, 9.0, 5.7],
-        "bulk_density": [1350.0, 1800.0, 1000.0, 1500.0],
         "clay_fraction": [0.8, 0.3, 0.1, 0.9],
         "litter_C_mineralisation_rate": [0.00212106, 0.00106053, 0.00049000, 0.0055],
         "litter_N_mineralisation_rate": [3.5351e-5, 7.0702e-5, 0.000183, 1.63333e-5],
diff --git a/tests/models/soil/test_soil_model.py b/tests/models/soil/test_soil_model.py
index 0ebc79d98..2ad30ef4a 100644
--- a/tests/models/soil/test_soil_model.py
+++ b/tests/models/soil/test_soil_model.py
@@ -556,7 +556,6 @@ def test_order_independance(
     # Add all the non-pool data into the new data object
     not_pools = [
         "pH",
-        "bulk_density",
         "soil_moisture",
         "matric_potential",
         "vertical_flow",
diff --git a/virtual_ecosystem/data_variables.toml b/virtual_ecosystem/data_variables.toml
index 560abd19e..d115ac2ae 100644
--- a/virtual_ecosystem/data_variables.toml
+++ b/virtual_ecosystem/data_variables.toml
@@ -712,13 +712,6 @@ name = "soil_p_pool_labile"
 unit = "kg P m^-3"
 variable_type = "float"
 
-[[variable]]
-axis = ["spatial"]
-description = "Bulk density of soil"
-name = "bulk_density"
-unit = "kg m^-3"
-variable_type = "float"
-
 [[variable]]
 axis = ["spatial"]
 description = "Soil pH values for each grid cell"
diff --git a/virtual_ecosystem/example_data/config/data_config.toml b/virtual_ecosystem/example_data/config/data_config.toml
index fabde284a..56fc1c76f 100644
--- a/virtual_ecosystem/example_data/config/data_config.toml
+++ b/virtual_ecosystem/example_data/config/data_config.toml
@@ -40,9 +40,6 @@ file_path = "../data/example_soil_data.nc"
 var_name = "pH"
 [[core.data.variable]]
 file_path = "../data/example_soil_data.nc"
-var_name = "bulk_density"
-[[core.data.variable]]
-file_path = "../data/example_soil_data.nc"
 var_name = "clay_fraction"
 [[core.data.variable]]
 file_path = "../data/example_soil_data.nc"
diff --git a/virtual_ecosystem/example_data/generation_scripts/soil_example_data.py b/virtual_ecosystem/example_data/generation_scripts/soil_example_data.py
index 3e14c96c9..400ddcc2f 100644
--- a/virtual_ecosystem/example_data/generation_scripts/soil_example_data.py
+++ b/virtual_ecosystem/example_data/generation_scripts/soil_example_data.py
@@ -16,9 +16,6 @@
 # Generate a range of plausible values (3.5-4.5) for the soil pH [unitless].
 pH_values = 3.5 + 1.00 * gradient / 64.0
 
-# Generate a range of plausible values (1200-1800) for the bulk density [kg m^-3].
-bulk_density_values = 1200.0 + 600.0 * gradient / 64.0
-
 # Generate a range of plausible values (0.27-0.40) for the clay fraction [fraction].
 clay_fraction_values = 0.27 + 0.13 * gradient / 64.0
 
@@ -104,7 +101,6 @@
 example_soil_data = Dataset(
     data_vars=dict(
         pH=(["x", "y"], pH_values),
-        bulk_density=(["x", "y"], bulk_density_values),
         clay_fraction=(["x", "y"], clay_fraction_values),
         soil_c_pool_lmwc=(["x", "y"], lmwc_values),
         soil_c_pool_maom=(["x", "y"], maom_values),