A turbidity time series, recorded every 15 min in a first-order Scottish stream for more than a year, along with two covariates (stage height and rainfall), is considered. Turbidity time series have complex dynamics because they are nonlinear, nonnormal, nonstationary, with a long memory, and present missing values.Given these issues the turbidity process is analyzed by Markov switching autoregressive models under the Bayesian paradigm. Since the multimodality of the posterior density novel evolutionary Monte Carlo (EMC) algorithms incorporating a few original features are developed to better traverse the posterior surface and escape from local basins of attraction. This because a population of chains are processed in parallel exchanging information one another and with different temperatures attached to each chain. These advanced EMC algorithms allow performing both Bayesian inference and model choice. Hence, it is possible to efficiently fit the actual data, reconstruct the sequence of hidden states, restore the missing values, and classify the observations into a few regimes, providing new insight on turbidity dynamics. A comparison with different nonlinear time series models is also proposed. Finally, a simulation study on the selection of the tuning factors of the EMC algorithm is presented.