Different turbulence closure schemes significantly influence the structure of Lagrangian residual velocity (LRV), yet the underlying mechanisms remain inadequately understood. Under constant eddy viscosity coefficient conditions, the different parameter β (the ratio of the eddy viscosity term to the local acceleration term) predominantly governs the LRV by modulating the Lagrangian mean barotropic pressure gradient, while the Lagrangian mean eddy viscosity term exerts a negative feedback effect. Under varying eddy viscosity conditions, dominant components of total Lagrangian mean eddy viscosity term vary across turbulence closure schemes and the governing mechanisms of LRV become increasingly complex. The pioneering research meticulously tracks particle motion from the zero-velocity initial phase, establishing an equivalence between LRV and the two-time Lagrangian integrals of the acceleration and other dynamic terms. The two-time Lagrangian integrals of local acceleration and horizontal advection terms play the dominant roles in the LRV, while the vertical advection terms provide supplementary effects. Under non-stratified conditions, the two-time Lagrangian integrals of barotropic component generally counterbalance the two-time Lagrangian integrals of eddy viscosity component. In stratified contexts, upper-layer LRV in the along-estuary direction is influenced primarily by the two-time Lagrangian integrals of barotropic component, contributing up to half of the total LRV, while the lower-layer inflow is significantly shaped by the combined interaction of two-time Lagrangian integrals of eddy viscosity and baroclinic components. In the cross-estuary direction, unlike the along-estuary direction, the two-time Lagrangian integrals of eddy viscosity component contribute negatively to LRV.