An improved method for computing the three-dimensional (3D) first-order Lagrangian residual velocity ( u L ) is established. The method computes tidal body force using the harmonic constants of the zeroth-order tidal current. Compared with using the tidal-averaging method to compute the tidal body force, the proposed method filters out the clutter other than the single-frequency tidal input from the open boundary and obtains u L that is more consistent with the analytic solution. Based on the new method, u L is calculated for a wide bay with a longitudinal topography. The strength and pattern of u L are mostly determined by the parts of the tidal body force related to the vertical mixing of the Stokes’ drift and the Coriolis effect, with a minor contribution from the advection effect. The geometrical shape of the bay can influence u L through the topographic gradient. The magnitude of u L increases with the increases in tidal energy input and vertical eddy viscosity and decreases in terms of the bottom friction coefficient.