Smoothed-Heaviside

Smoothed Heaviside surface boundary condition for coupling surface and subsurface water flow, from Wang et al. (2020) (Advances in Water Resources 137, 103499), Eqs. 3 and 5. The component switches a single surface node continuously between a ponded (saturated) boundary and an unsaturated boundary as the pressure head h crosses zero, replacing a hard step with the smooth approximations eta_omega(h) = (1/pi)*atan(h/omega) + 1/2 (smoothed Heaviside, Eq. 5a) and delta_omega(h) = (1/pi)*omega/(omega^2 + h^2) (smoothed Dirac delta, Eq. 5b; delta_omega = d eta_omega/dh). The head obeys the boundary condition d/dt[eta_omega(h)*h] = P - I_infil (Eq. 3), stated as the explicit ODE D(h) = (P - I_infil)/(delta_omega*h + eta_omega). h is the single integrated state; eta_omega and delta_omega are algebraic (observed). Migrated from EnvironmentalTransport.jl@a48559c4 src/surface_runoff.jl:232 (`HeavisideBoundaryCondition` @component). Because the head equation conserves eta_omega(h)*h - (P - I_infil)*t, every ODE-integration test below is checked against the exact conservation law, not merely against a prior solver run.