Comprehensive modelling of the Yellowstone caldera: insights into thermal evolution, magmatic behavior, and lithospheric strength

The Yellowstone Volcanic Complex (YVC) in Yellowstone National Park (Wyoming, USA) attracts significant geological interest as one of the largest active continental silicic volcanic fields in the world. Despite extensive research on its high heat flow and abundant geothermal features, a detailed quantitative analysis of the brittle-ductile transition remains absent. This study aims to deepen the understanding of the subsurface geological and geophysical properties of the Yellowstone area, with a particular focus on developing an optimized lithospheric thermal profile, essential for reliable rheological and lithospheric strength analyses. Initially, the Curie isothermal surface depth was extensively mapped using high-resolution aeromagnetic data and innovative spectral analysis techniques. This mapping revealed a shallow Curie isothermal surface, ranging from 2 km to 4 km beneath the YVC. The retrieved iso-Curie depth was subsequently used as a key constraint to validate a 3D stationary Finite Element (FE) thermal model. Specifically, this isothermal surface served as experimental data for optimizing the thermal state of the crust through the Optimization Module in COMSOL Multiphysics® 6.2. The 3D model of the Yellowstone lithosphere covers approximately 40 km², with a lithospheric thickness of about 35.000  km. The domain is subdivided into three litho-thermal units: the upper crust, the lower crust, and the magmatic body. The geometry of the magmatic heat source was derived from tomographic data and incorporated into COMSOL Multiphysics® to create a consistent subsurface image of the magmatic heat source. The thermal state of the crust was simulated using the Heat Transfer in Solids Module under a purely conductive regime. To further validate the thermal model, the DBSCAN clustering algorithm was applied to analyze seismic data. A comprehensive rheological model was also developed to delineate the brittle-ductile transition within the lithospheric volume. The results revealed a brittle region well-aligned with the earthquake distribution and a complex, layered ductile zone structure, reflecting the stratified nature of the local lithospheric architecture. This study contributes to a deeper understanding of the YVC’s subsurface dynamics, offering insights into its complex geodynamic processes and providing methodologies applicable to similar studies in other volcanic and geothermally active regions with large calderas.