This paper presents the mathematical formulation of an integrated surface-subsurface-overland model for flow and transport of thermal energy and salinity. The model contains three major components: a 3D surface water module, a 2D overland module, and a 3D subsurface module. The surface water module simulates flow and transport in the main channel of rivers/estuaries based on the 3D Navier-Stokes equations with or without hydrostatic assumptions. The moving free surface is explicitly handled by solving the kinematic boundary condition equation and a moving grid method with a node-repositioning algorithm is used to track the deformation of the water surface. The surface water transport module solves the energy equation for spatial-temporal distributions of temperature and the mass transport equations for the salinity field. The numerical solution is obtained using the finite element method or the mixed Lagrangian-Eulerian (particle tracking) and finite element method. The overland module solves the 2D simplified diffusive wave equations, and the transport equations for temperature and salinity using the finite element method. The groundwater module adopts the modified Richard's equation, the Darcy's law, and the transport equations to model the flow and transport of temperature and salinity through saturated-unsaturated porous media using the finite element method. The integrated model also enables the dynamic surface-subsurface, surface-overland, and overland-subsurface interactions through the couplings of flows and transports at the interface. The integrated model has been tested and applied to Loxahatchee Estuary for the investigation of floodplain habitat hydroperiod in the Northwest Fork of the river.