We demonstrate the potential of a recently developed grid-based eikonal solver for tracking phases comprising reflection branches, transmission branches, or a combination of these, in 3D heterogeneous layered media. The scheme is based on a multi-stage fast marching approach that reinitialises the wavefront from each interface it encounters as either a reflection or transmission. The use of spherical coordinates allows wavefronts and traveltimes to be computed at local, regional, and semi-global scales. Traveltime datasets for a large variety of seismic experiments can be predicted, including reflection, wide-angle reflection and refraction, local earthquake, and teleseismic. A series of examples are presented to demonstrate potential applications of the method. These include: (1) tracking active and passive source wavefronts in the presence of a complex subduction zone; (2) earthquake hypocentre relocation in a laterally heterogeneous 3D medium; (3) joint inversion of wide-angle and teleseismic datasets for P-wave velocity structure in the crust and upper mantle. Results from these numerical experiments show that the new scheme is highly flexible, robust and efficient, a combination seldom found in either grid- or ray-based traveltime solvers. The ability to track arrivals for multiple data classes such as wide-angle and teleseismic is of particular importance, given the recent momentum in the seismic imaging community towards combining active and passive source datasets in a single tomographic inversion.