diff --git a/ctapipe/containers.py b/ctapipe/containers.py
index bcccaef7f40..ba8bddf8fe2 100644
--- a/ctapipe/containers.py
+++ b/ctapipe/containers.py
@@ -586,8 +586,11 @@ class SimulatedShowerContainer(Container):
     core_x = Field(nan * u.m, "Simulated core position (x)", unit=u.m)
     core_y = Field(nan * u.m, "Simulated core position (y)", unit=u.m)
     h_first_int = Field(nan * u.m, "Height of first interaction", unit=u.m)
-    x_max = Field(
-        nan * u.g / (u.cm**2), "Simulated Xmax value", unit=u.g / (u.cm**2)
+    x_max = Field(nan * u.g / u.cm**2, "Simulated Xmax value", unit=u.g / u.cm**2)
+    starting_grammage = Field(
+        nan * u.g / u.cm**2,
+        "Grammage (mass overburden) where the particle was injected into the atmosphere",
+        unit=u.g / u.cm**2,
     )
     shower_primary_id = Field(
         np.int16(np.iinfo(np.int16).max),
@@ -692,9 +695,6 @@ class SimulationConfigContainer(Container):
     core_pos_mode = Field(
         nan, description="Core Position Mode (0=Circular, 1=Rectangular)"
     )
-    injection_height = Field(
-        nan * u.m, description="Height of particle injection", unit=u.m
-    )
     atmosphere = Field(nan * u.m, description="Atmospheric model number")
     corsika_iact_options = Field(
         nan, description="CORSIKA simulation options for IACTs"
diff --git a/ctapipe/io/simteleventsource.py b/ctapipe/io/simteleventsource.py
index 5dd00a226e7..235f8edd737 100644
--- a/ctapipe/io/simteleventsource.py
+++ b/ctapipe/io/simteleventsource.py
@@ -88,6 +88,7 @@
 
 _half_pi = 0.5 * np.pi
 _half_pi_maxval = (1 + 1e-6) * _half_pi
+_float32_nan = np.float32(np.nan)
 
 
 def _clip_altitude_if_close(altitude):
@@ -117,7 +118,7 @@ class MirrorClass(enum.Enum):
     DUAL_MIRROR = 2
 
 
-X_MAX_UNIT = u.g / (u.cm**2)
+GRAMMAGE_UNIT = u.g / (u.cm**2)
 NANOSECONDS_PER_DAY = (1 * u.day).to_value(u.ns)
 
 
@@ -990,7 +991,6 @@ def _parse_simulation_header(self):
             max_scatter_range=mc_run_head["core_range"][1] * u.m,
             min_scatter_range=mc_run_head["core_range"][0] * u.m,
             core_pos_mode=mc_run_head["core_pos_mode"],
-            injection_height=mc_run_head["injection_height"] * u.m,
             atmosphere=mc_run_head["atmosphere"],
             corsika_iact_options=mc_run_head["corsika_iact_options"],
             corsika_low_E_model=mc_run_head["corsika_low_E_model"],
@@ -1053,6 +1053,9 @@ def _fill_simulated_event_information(array_event):
             core_x=u.Quantity(mc_event["xcore"], u.m),
             core_y=u.Quantity(mc_event["ycore"], u.m),
             h_first_int=u.Quantity(mc_shower["h_first_int"], u.m),
-            x_max=u.Quantity(mc_shower["xmax"], X_MAX_UNIT),
+            x_max=u.Quantity(mc_shower["xmax"], GRAMMAGE_UNIT),
             shower_primary_id=mc_shower["primary_id"],
+            starting_grammage=u.Quantity(
+                mc_shower.get("depth_start", _float32_nan), GRAMMAGE_UNIT
+            ),
         )
diff --git a/ctapipe/io/tests/test_simteleventsource.py b/ctapipe/io/tests/test_simteleventsource.py
index 04cf133d414..e9e7ea31f9a 100644
--- a/ctapipe/io/tests/test_simteleventsource.py
+++ b/ctapipe/io/tests/test_simteleventsource.py
@@ -591,3 +591,11 @@ def test_float32_pihalf(sign):
     assert shower.alt.value == sign * np.pi / 2
     # check we cana create a Latitude:
     Latitude(shower.alt.value, u.rad)
+
+
+def test_starting_grammage():
+    path = "dataset://lst_muons.simtel.zst"
+
+    with SimTelEventSource(path, focal_length_choice="EQUIVALENT") as source:
+        e = next(iter(source))
+        assert e.simulation.shower.starting_grammage == 580 * u.g / u.cm**2