One of the most important causes of the freshwater shortage in estuarine area is the increasing seawater intrusion into the river. To simulate seawater intrusion properly, two important factors should be considered. One is the bidirectional and time-dependent coupling effects between river discharges and tidal forces. The other is the three-dimensional and stratified structure of dynamic processes involved. However, these two factors have rarely been investigated simultaneously, or they were often simplified in previous researches, especially for the estuary connected with an upstream river network through multiple outlets such as the Pearl River Estuary (PRE). In order to consider these two factors, a numerical modeling system, which couples a one-dimensional river network model with a three-dimensional unstructured-grid Finite-Volume Coastal Ocean Model (FVCOM), has been developed and successfully applied to the simulation of seawater intrusion into rivers emptying into the PRE. By treating the river network with a one-dimensional model, computational efficiency has been improved. With coupling 1D and 3D models, the specification of up-stream boundary conditions becomes more convenient. Simulated results are compared with field measured data. Good agreement indicates that the modeling system may correctly capture the physical processes of seawater intrusion into rivers.