Over-exploitation of groundwater in coastal aquifers is one of the major reasons for saltwater intrusion. Saltwater intrusion contaminates the freshwater resource in coastal aquifers by increasing the salinity levels of the groundwater. A small quantity of saltwater is enough to contaminate the large quantity of freshwater in a coastal aquifer. Therefore, saltwater intrusion should be considered seriously. Numerical models are needed to enhance the understanding of saltwater intrusion and its related phenomena extensively. The Motooka region in the Fukuoka prefecture in Japan is a coastal area, where groundwater is the main water resource for green house agriculture and domestic use. Over-exploitation of groundwater in this coastal region has resulted in saltwater intrusion and thus in the contamination of the freshwater aquifer. In addition to the more obvious effects of saltwater intrusion, fluctuations in salinity caused by such intrusion is a crucial problem to address, since even slight changes in salinity of the water use for agricultural purposes significantly affect the crop's growth and yield. So far, a research on the salinity fluctuations with groundwater pumping and their effects on seasonal recharge of groundwater in the Motooka region has not been conducted. Therefore, in this study a three-dimensional density-dependent solute transport flow model is developed to simulate the salinity fluctuations due to groundwater pumping. In the present numerical study, the emphasis is on the development of conceptual, mathematical and numerical model of variable density flow and solute transport and its application to simulate the salinity fluctuation due to groundwater pumping at different rates. The model is based on the “transition zone” approach, which considers the interface as a miscible zone where freshwater and saltwater is mixing while maintaining a density gradient across the freshwater/saltwater interface. The transition zone approach requires simultaneous solutions of the governing water flow and solute transport equations. To this end, the model incorporates three fundamental equations in flow and solute transport, namely Darcy's law, general groundwater flow equation and the advection dispersion solute transport equation. The groundwater flow equation and solute transport equation are coupled by the equation of state to produce the salt concentration at each time step for whole flow domain. The finite difference method is used as the numerical technique to solve the partial differential equations of flow and transport under an implicit scheme. The method of characteristics is applied to solve the advection term in the solute transport equation. A non-uniform discretized grid system is adopted in the flow domain allocating relatively small grid sizes to pumping well locations. To achieve reliable results, relevant and important hydro-geological parameters are assigned to the numerical model after considering the hydrological situation of the Motooka region. Different boundary conditions are assigned considering the dominant hydrological processes those are believed to be in effect in the selected area. The numerical results obtained from the model demonstrate the salinity variation due to groundwater pumping and seasonal recharge rates from year 2001 to 2007 under the influence of saltwater intrusion. The results also reveal that model is capable of correctly simulating the physical processes. A comparison of the measured and modeled electric conductivities shows reasonable agreement.