Magma is commonly thought to be stored as a crystal-rich mush within vertically extensive, crustal storage regions. A key unknown is how to remobilise and erupt such crystal-rich material, and whether the growth of gas bubbles within the mush could promote remobilisation. In order to investigate this, we need improved constraints on the timing of volatile saturation in magmas. The mineral apatite represents a potentially useful record of pre-eruptive magmatic volatiles, but data interpretation is complex because exchange reactions control the volatile partitioning. Model solutions are therefore non-unique. Here, we present a numerical forward modelling program with a sensitivity analysis function, which addresses non-uniqueness by identifying alternative sets of starting parameters that match a target compositional trend through a population of apatite crystals. The model is applied to a new dataset of volatiles in apatite from the 12.9 ka Laacher See eruption, Eifel volcanic region, Germany. The results indicate that the magma was initially strongly volatile-undersaturated and became saturated through progressive crystal fractionation. Apatite crystals are not in volatile or trace element equilibrium with their carrier melts, indicating dispersal of crystals into different chemical environments. Consideration of apatite diffusivities suggests that this reorganisation occurred shortly before eruption. Our modelling results also allow us to constrain directly the amount of pre-eruptive magmatic vapour emitted during the explosive eruption, highlighting the importance of considering the behaviour of halogens during magma storage. Overall, our approach confirms the value of measuring apatite volatile contents and highlights the potential of this method to provide quantitative constraints on magmatic evolution and storage conditions.