This paper presents the development and application of an integrated surface water and groundwater model to simulate hydrodynamics and thermal and salinity in tidal water bodies and subsurface media. The hydrodynamic module for tidal waters solves three-dimensional Navier-Stokes equations with or without the hydrostatic assumptions. The Richards equation is used to simulate the subsurface flow in both vadose and saturated zones. The Boussinesq approximation is employed to deal with the buoyancy force due to temperature and salinity variations. The moving free surface is explicitly handled by solving the kinematic boundary condition equation using a node-repositioning algorithm. The transport module solves the energy equation for temperature distribution and the mass transport equations for the salinity fields. The Arbitrary Lagrangian-Eulerian (ALE) representation is adopted for all transport equations including momentum transport. The solution is obtained with finite element methods or a combination of finite element and Semi-Lagrangian (particle tracking) methods. The model has been calibrated with tidal and salinity measurements at five locations in Loxahatchee River Estuaries. The calibration of groundwater flow and subsurface salinity transport is not as satisfactory. Further study to downscale the MODFLOW hydraulic conductivity is needed to better reflect the fine grid used to discretize the groundwater module in the present study. A better conceptual model to deal with the interface between surface water and subsurface water and to handle flood plain is warranted.