Although variations in the density,ρ, are naturally incorporated into finite-difference parabolic equation (PE) solvers, split-step PE algorithms traditionally account for density changes by adding to the refractive index terms containing derivatives ofρ. As a consequence, geoacoustic density profiles that contain step discontinuities at layer interfaces must be smoothed appropriately before these extra terms can be evaluated. In this paper, a new hybrid method is proposed for treating density inhomogeneities in the split-step PE. This approach involves splitting the differential operator into density-independent and density-dependent components. While the former terms may be evaluated with the split-step Fourier technique, the influence of density changes is handled through a finite-difference procedure. Such an algorithm can be easily implemented in recently proposed hybrid split-step/finite-difference and split-step/Lanczos PE solvers. [J. Acoust. Soc. Am.96, 396–405 (1994)].