Phosphorus (P) enrichment in soils has been documented in the Santa Fe River watershed (SFRW, 3585 km) in north-central Florida. Yet the environmental factors that control P distribution in soils across the landscape, with potential contribution to water quality impairment, are not well understood. The main goal of this study was to develop soil-landscape P models to support a "precision soil conservation" approach combining fine-scale (i.e., site-specific) and coarse-scale (i.e., watershed-extent) assessment of soil P. The specific objectives were to: (i) identify those environmental properties that impart the most control on the spatial distribution of soil Mehlich-1 extracted P (MP) in the SFRW; (ii) model the spatial patterns of soil MP using geostatistical methods; and (iii) assess model quality using independent validation samples. Soil MP data at 137 sites were fused with spatially explicit environmental covariates to develop soil MP prediction models using univariate (lognormal kriging, LNK) and multivariate methods (regression kriging, RK, and cokriging, CK). Incorporation of exhaustive environmental data into multivariate models (RK and CK) improved the prediction of soil MP in the SFRW compared with the univariate model (LNK), which relies solely on soil measurements. Among all tested environmental covariates, land use and vegetation related properties (topsoil) and geologic data (subsoil) showed the largest predictive power to build inferential models for soil MP. Findings from this study contribute to a better understanding of spatially explicit interactions between soil P and other environmental variables, facilitating improved land resource management while minimizing adverse risks to the environment.