arXiv:2610.23177v1 Announce Type: new Abstract: Mathematical modelling of burn injury healing involves coupling biological processes with mechanical deformation, resulting in computationally stiff, nonlinear partial differential equations that severely limit the scalability of existing numerical solutions. We present a computational framework addressing these challenges. By nondimensionalizing the model and performing scaling analysis, we show that chemotaxis and inertia are negligible, permitting a quasi-static momentum approach. To enforce non-negativity constraints, we implement a Variational Inequality formulation using the Finite Element Method in FEniCSx, and verify it via numerical convergence studies and mass balances. We identify high-order parameter interactions, specifically equilibrium collagen density and signalling dynamics, as the primary drivers of solver stability. A sensitivity analysis of nearly four million simulations reveals that isolated parameter variations predominantly govern the maximum wound contraction. In contrast, higher-order interactions determine the timing of maximum contraction and the permanent contracture. However, no significant pairwise parameter interactions were found, suggesting that the model's behaviour is governed by higher-order interactions, highlighting the inherent complexity of the wound-healing process.