A density‐dependent, variably saturated groundwater flow and solute transport model was used to investigate the influence of swash motions on subsurface flow and moisture dynamics in beach aquifers with heterogeneous distributions of hydraulic conductivity (K) and capillarity. The numerical simulations were performed within a Monte Carlo framework using field measurements conducted in the swash zone of a sandy beach. Our results show that heterogeneous capillarity causes spatially variable capillary rise above the groundwater table. In response to swash motions, heterogeneity creates capillary barriers that result in pockets of elevated moisture content beneath the swash zone. These moisture hotspots persist within the unsaturated zone even at ebb tide when the swash motions recede seaward. Heterogeneous capillarity also results in highly tortuous preferential flow paths and alters the flow rates from the sand surface to the water table. Heterogeneous K greatly enhances the seawater infiltration into the swash zone and modulates its spatial distribution along the beach surface. Due to heterogeneous K and capillarity, complex mixing patterns emerge. Both strain‐dominated and vorticity‐dominated flow regions develop and dissipate as tides and waves move across the beach surface. Complex mixing patterns of seawater percolating from the swash zone surface to the water table, with localized areas of high and low mixing intensities, are further demonstrated by analysis of dilution index. Our findings reveal the influence of geologic heterogeneity on swash zone moisture and flow dynamics, which may have important implications for sediment transport and chemical processing in beach aquifers.