<p>The Parker–Oldenburg (ParkerO) method is a widely used physics-based framework for inverting density-interfaces from gravity data. However, its conventional assumption–a single interface with uniform density contrast–fails to represent multilayered, vertically heterogeneous geological settings. To address this limitation, we introduce an Enhanced Varying-Density Parker–Oldenburg (EVD-PO) formulation that explicitly incorporates multiple density-interfaces and depth-dependent, vertically varying density contrasts into gravity anomaly modeling. By rederiving the Fourier-domain forward operator and embedding it within a stabilized iterative scheme, EVD-PO enables accurate and computationally efficient recovery of density-interface geometry. We validate the method in two realistic western Pacific case studies: Moho depth recovery and seafloor topography reconstruction, using independent reference data for comparison. In both applications, EVD-PO resolves finer structural details and yields lower errors compared to the reference data. These results demonstrate that EVD-PO effectively captures the nonlinear contributions from interface relief and density heterogeneity, offering a robust, high-accuracy tool for imaging complex subsurface structures from gravity data.</p>