River water pollution is one of the most important environmental issues. Advection–dispersion equations are used to study the temporal changes in pollutant concentration along the study river reach. The use of advection–dispersion equations in investigating how the concentration of pollution is transformed requires a lot of data including river cross-section characteristics, dispersion coefficient, and upstream and downstream boundary conditions, etc. therefore, the corresponding calculations are very costly, difficult and time-consuming. In the present study, instead of using the mentioned equations, the linear Muskingum method (used in previous studies for flood routing) and the particle swarm optimization (PSO) algorithm was used for the first time to calculate the temporal changes in pollution concentration at different stream locations. The presented solution in the presented study is very accurate and only requires the temporal changes in concentration in the upstream and downstream of the study river reach and for this reason, it is very low-cost and easy to use and requires less time to collect data and perform calculations. In the proposed method, the parameters (X, K, ∆t) of the linear Muskingum method were optimized using the PSO algorithm, and by dividing the temporal changes in the input concentration into three areas of the beginning (the input concentration is greater than the output concentration), the peak (the maximum input and output concentrations) and the end (the output concentration is greater than the input concentration) areas, the accuracy of the calculations increased. The mentioned method was studied for different lengths (first case of x = 50 m (up) and x = 75 m (down), second case of x = 50 m (up) and x = 100 m (down), third case of x = 75 m (up) and x = 100 m (down)) and the mean relative error (MRE) of the total, peak area and the relative error of the maximum concentration using constant parameters for the first case were calculated as 7.08, 1.02, and 2.34 percent, for the second case as 7.41, 11.06 and 6.69 percent, and for the third case as 6.75, 3.59 and 5.42 percent, respectively. If three parameters of (X, K, ∆t) are used, the mentioned values improved by 31.3, 63.7 and 65.5 percent, respectively compared to the case of using constant parameters.