A two-dimensional volume of fluid (VOF) method was used to simulate liquid movement within unsaturated fracture intersections and fracture networks. Fragmentation and merging of liquid interfaces are automatically handled, without resorting to adaptive mesh refinement or interface repairing algorithms. Wetting effects are modeled by imposing different contact angles depending on whether the liquid interface is advancing, receding, or essentially stationary. The numerical examples and experimental validation presented in this study clearly demonstrate the potential value of the VOF method for the modeling of liquid motion in unsaturated fractures, fracture junctions and fracture networks.