AbstractFor theoretical studies, reaction norm evolution in a changing environment can be modeled by means of the multivariate breeder's equation, with the reaction norm parameters treated as traits in their own right. This is, however, not a feasible approach for the use of field data, where the intercept and slope values are not available. An alternative approach is to use infinite‐dimensional characters and smooth covariance function estimates found by, e.g., random regression. This is difficult because of the need to find, for example, polynomial basis functions that fit the data reasonably well over time, and because reaction norms in multivariate cases are correlated, such that they cannot be modeled independently. Here, I present an alternative approach based on a multivariate linear mixed model of any order, with dynamical incidence and residual covariance matrices that reflect the changing environment. From such a mixed model follows a dynamical BLUP model for the estimation of the individual reaction norm parameter values at any given parent generation, and for updating of the mean reaction norm parameter values from generation to generation by means of Robertson's secondary theorem of natural selection. This will, for example, make it possible to disentangle the microevolutionary and plasticity components in climate change responses. The BLUP model incorporates the additive genetic relationship matrix in the usual way, and overlapping generations can easily be accommodated. Additive genetic and environmental model parameters are assumed to be known and constant, but it is discussed how they can be estimated by means of a prediction error method. The identifiability by the use of field or laboratory data containing environmental, phenotypic, fitness, and additive genetic relationship data is an important feature of the proposed model.