This paper presents an incremental variational method to assimilate the observed tidal harmonic constants using a frequency domain linearized shallow water equation. A cost function was constructed with tidal boundary conditions and tidal forcing as its control (independent) variables. To minimize the cost function, optimal boundary conditions and tidal forcing were derived using a conventional dual 4-Dimensional Variational (4D-Var) Physical-space Statistical Analysis System. The tangent linear and adjoint model were solved by using a finite element method. By adapting the incremental form, the variational method streamlines the workflow to provide the incremental correction to the boundary conditions and tidal forcing of a hydrodynamic forward model. The method was tested for semi-diurnal M2 tides in a regional sea with a complex tidal system. The results demonstrate a 65–72% reduction of tidal harmonic constant vector error by assimilating the observed M2 tidal harmonic constants. In addition to improving the tides of a hydrodynamic model by optimizing boundary conditions and tidal forcing, the method computes a spatially varying uncertainty of individual tidal constituents in the model. The method provides a versatile tool for mapping the spatially continuous tides and currents in coastal and estuarine waters by assimilating the harmonic constants of individual tidal constituents of observed tides and currents.