This paper presents a robust and scalable surface tracking algorithm for fluid dynamics simulations on non-body-fitted meshes, such as those used in embedded boundary methods. The algorithm is specifically designed to address scenarios where multiple solid surfaces come into contact or interpenetrate — situations that are difficult to handle with conventional body-fitted approaches. It computes three types of geometric quantities: (1) subdomain classification of background mesh nodes, (2) intersections between edges in the background mesh and the embedded boundaries, and (3) shortest distances from background mesh nodes to embedded boundaries. The algorithm builds on fundamental tools from computational geometry, including flood-fill, ray-tracing, and continuous collision detection [1]. When surface penetration is detected, overlapped regions in the background mesh are identified using the edge-surface intersections from the involved surfaces. Then, embedded surface elements that bound these regions are selectively excluded to reflect the physical boundary evolution according to user-specified modeling options. To improve computational efficiency and scalability, the algorithm utilizes space-partitioning trees (e.g., k-d trees), bounding volume hierarchies, and localized subdomain queries. The performance of the interface tracking algorithm is demonstrated through several numerical examples involving surface contact and penetration. REFERENCES [1] J.D. Foley, A. van Dam, S.K. Feinder, and J.F. Hughes, Computer Graphics: Principles and Practice, 2nd ed., Addison Wesley, (1990).