Realistic lunar tides of the Great Bay Estuarine System have been simulated using a fixed boundary finite element numerical model as described by Ip et al. <1998>. It is a two-dimensional, nonlinear, time-stepping model with a groundwater component coupled to a kinematic force balance to facilitate the realistic drainage and filling of elements during ebb and flood, respectively. The numerical model reproduces the observed M2 tides as described by Swift and Brown <1983>, and it successfully captures qualitatively correct residual currents and transports, realistic mass-conserving flooding and dewatering of the tidal flats, and current asymmetry between flood and ebb. The simulation results were sensitive to local bathymetry and the implemented friction law. The accuracy of the model is demonstrated by comparison with the 1975 Great Bay study <Swift and Brown, 1983> in terms of tidal elevations at 14 tidal stations and 4 cross-sectionally averaged current measurements in the estuary. Quantitatively, the model results show good agreement with observations, displaying correlation coefficients of ≥0.96 in surface elevation and ≥0.95 in averaged current, with average RMS errors of 0.12 m and 0.26 m s-1, respectively. In addition, tidal flat hydrodynamics, characteristic distributions of residual current, sediment bed load transport, and influence of topography on the overall circulation in the region are also discussed.