<p>The accuracy of GNSS-Acoustic (GNSS-A) seafloor positioning is much affected by spatiotemporal variability in the ocean Sound Speed Field (SSF). The conventional depth-invariant horizontal sound speed gradient models represent the water-column sound speed structure using a single effective depth-averaged gradient, which cannot adequately describe the vertical heterogeneity of horizontal gradients and therefore introduce model errors and positioning biases. To address this issue, we propose a layered sound speed gradient model and a corresponding joint inversion framework for precise GNSS-A seafloor positioning. The proposed method parameterizes horizontal sound speed gradients in multiple depth layers and couples adjacent layers through depth-weighted interlayer continuity constraints, which jointly estimates seafloor transponder coordinates, the depth-invariant temporal perturbation, and layer-wise horizontal gradients. Simulation results under depth-dependent horizontal gradient scenarios show that the conventional single-layer model introduces systematic positioning biases, whereas the proposed layered model significantly improves positioning accuracy, accurately estimates the temporal perturbation and layer-wise horizontal gradients, and is robust for a broad range of layering configurations and constraint parameters. Field experiments using real GNSS-A observations further demonstrate the practical value of the proposed method. The results show that the proposed approach improves the fitting of acoustic observations, maintains short-term repeatability, and yields consistent multi-epoch coordinate time series together with reasonable site-velocity estimates. These findings indicate that incorporating a layered horizontal gradient structure can improve GNSS-A seafloor positioning and provides an interpretable modeling and inversion framework for long-term seafloor deformation monitoring and characterization of ocean environmental variability.</p>