The particle finite element method, or PFEM, has shown great capabilities to simulate free surface flows with strong topological changes. Indeed, thanks to its inherent capability to track the evolution of moving boundaries and deforming domains, natural phenomena such as dam breaks, landslides, sloshing and many others can be simulated. An important challenge of this method remains to adapt the mesh as it is deformed. Indeed, due to the continuous displacement of the particles, a frequent remeshing step is required. In the case of extremely strong topological changes, one may even argue that remeshing is necessary at each time step. In the PFEM, remeshing means inserting and removing nodes, i.e., the particles of the simulation, to maintain high quality elements throughout the domain, without altering the shape of the fluid. In [1], we presented a method to perform quality-based Delaunay refinement in the PFEM for two-dimensional simulations. We showed that this not only greatly improves the smoothness of the free surface boundaries but also reduces mass conservation errors caused by remeshing. In the present work, we extend the approach to three dimensions. First, we use a modified version of the well-known alpha-shape algorithm to detect the domain, by using the advected mesh of the previous time step as a predicate for detecting the new domain. Next, we refine the closed surface with a longest-edge splitting algorithm. Finally, the volume mesh is refined with a quality-based node insertion at element circumcenters. The use of a size field allows us to perform adaptive refinement during the simulation. We will present the mesh refinement technique and show its capabilities on simulations of complex free surface flows.