Total Petroleum Hydrocarbon (TPH) contamination poses a long-term threat to soil quality and ecosystem health in oil-producing regions. Accurate characterization of TPH drivers is challenging because environmental responses are frequently nonlinear, while available datasets are often sparse and uncertain. In this study, we evaluated whether incorporating process-informed shape constraints can improve statistical modelling of TPH concentrations. A common Gamma–log modelling framework was implemented and used to compare a Generalized Linear Model (GLM), a Generalized Additive Model (GAM), and a Shape-Constrained Additive Model (SCAM). The evaluation was performed using a site-informed simulation scenario that reflects established environmental expectations, in which TPH increases with Total Organic Carbon (TOC) and decreases with both distance from oil wells and soil depth. A synthetic dataset comprising 600 observations was generated through constrained random sampling over empirically informed covariate ranges. Model performance was examined using residual diagnostics, calibration analysis, and stratified 10-fold cross-validation. In contrast, uncertainty was assessed through Bias-Corrected and Accelerated (BCa) bootstrap intervals (B = 5,000) and local one-at-a-time sensitivity analysis (±10%). Results showed that SCAM achieved the most coherent monotonic response patterns and the strongest calibration performance, highlighting the value of integrating scientifically justified shape constraints into environmental regression models. Although conclusions remain conditional on the adopted simulation framework, the study demonstrates how constrained nonlinear models can improve interpretability and robustness in data-limited contamination assessments.