SUMMARYMacro-hybrid penalized ÿnite element approximations are studied for steady ÿltration problems with seawater intrusion. On the basis of nonoverlapping domain decompositions with vertical interfaces, sections of coastal aquifers are decomposed into subsystems with simpler geometries and small scales, interconnected via transmission conditions of pressure and ux continuity. Corresponding local penalized formulations are derived from the global penalized variational formulation of the two-free boundary ow problem, with continuity transmission conditions modelled variationally in a dual sense. Then, macrohybrid ÿnite element approximations are derived for the system, deÿned on independent subdomain grids. Parallel relaxation penalty-duality algorithms are proposed from ÿxed-point problem characterizations. Numerical experiments exemplify the macro-hybrid penalized theory, showing a good agreement with previous primal conforming penalized ÿnite element approximations (Comput. Methods Appl. Mech. Engng. 2000; 190:609-624).