Networks for monitoring hazardous environmental spatio-temporal processes can evolve over time and eventually yield unrepresentative data. Estimates of model parameters for those processes may become biased. The paper presents a method for jointly modelling and fitting a (Gaussian) environmental process and (binary) selection process where the number of sites is exceptionally large, ruling out use of standard Bayesian computational methods for model fitting. Instead the paper relies on the integrated nested Laplace approximation. Within that framework, the selection model includes a number of important features: retention (once selected, sites tend to remain in); repulsion (new sites will not placed near established ones); noncompliance detection (sites tend to be placed where environmental hazards are large). A model residual term picks up unknown selection factors that may vary temporally, and independently from the environmental process. A case study shows that closures of a large number of monitoring stations in the United Kingdom while ambient levels of black smoke declined, were made in a preferential way; ignoring this leads to biased estimation of model parameters.