For surface seismic experiments we require Green's functions that are the response to an impulse at a single point on the surface. However the phase shift algorithm requires initial conditions at all locations on the surface,. Since horizontal waves, reflected waves, and overturned waves can all contribute to the surface wavefield it appears that you must have already modeled the data in order to find the initial condition. To avoid this problem I calculate a dip limited Green's function so that only events traveling at less than some angle to the vertical (say 80 degrees) are propagated. I assume that none of this energy reaches the surface. The initial condition is then just an impulse at the source location that is dip filtered to an appropriate dip range. While this does not calculate the Green's function correctly at high dips, it should be adequate for the dip ranges typical in seismic data.
Figure shows the amplitude and phase plots for a mono-frequency Green's function calculated for a constant velocity medium. The amplitude should be radially symmetric but it decays near the surface because of the dip limitation. Figure shows one time slice through the full Green's function in the (x,z,t) domain. This was calculated by inverse transforming a complete set of mono-frequency Green's functions. The dip limitation in the initial conditions is clearly visible as an amplitude decay of the wavefront at steeper angles.
Because the simple one-way equation only propagates energy downwards, it cannot correctly model overturned waves. Claerbout 1985 explained how to overcome this problem by saving the upgoing energy at each depth step, and then use an upwards extrapolation phase to pick up the upgoing energy and propagate it to the surface. This method of modeling overturned waves clearly requires twice as much computation as a simple downward extrapolation.