Various strategies are proposed for enforcing stress-free boundary conditions on deformed domains using a weak-form-based discretization in a Cartesian frame of reference. Due to the irregularity of the computational domain and the particular type of boundary condition, a coupling between the Cartesian velocity components is introduced. As such, a different computational kernel is required for the numerical solution of the associated vector Helmholtz equation, as contrasted to what is used for the scalar unknowns of the equations. Three approaches are presented, aimed towards the exact or approximate implementation of zero tangential stress (traction) at the deformed boundary while ensuring a decoupling of the velocity components in the solution of the vector Helmholtz equation. Two of these strategies, those which approximate the free slip boundary condition, are applied to the propagation of an internal solitary wave (ISW) of depression over a deformed bathymetry. The spatial structure and amplitude of the resultant pseudo-traction, which is accurately predicted by a simple scaling estimate, are explored as a function of the ISW-based Reynolds number, Re. For the Re values considered, the pseudo-traction is negligible with respect to the corresponding no-slip tangential shear stress. The pseudo-traction-induced, time-integrated loss of ISW energy is found to be significantly weaker than the associated viscous dissipation in the interior water column. Although, at laboratory-scale or oceanic Re, an approximate free-slip boundary condition is found to yield negligible pseudo-traction, this might not be the case when an elevated eddy viscosity is used in this context as a surrogate for no-slip turbulent bottom boundary layer dynamics.