A method for simulating forward wavefront propagation in heterogeneous tissue is discussed. The intended application of this method is for the study of aberration produced when performing ultrasound imaging through a layer of soft tissue. A one-way wave equation that permits smooth variation in all acoustically important variables is derived. This equation also describes tissue exhibiting nonlinear elasticity and arbitrary frequency-dependent relaxation. A numerical solution to this equation is found by means of operator splitting and propagation along the spatial depth coordinate. The numerical solution is accurate when compared to analytical solutions for special cases, and when compared to numerical solutions of the full wave equation by other methods. The presented implementation provides a fast numerical method for studying the impact of aberration in medical ultrasound imaging through soft tissue--both on the transmitted beam and the nonlinearly generated harmonic beam.