This paper proposes an aggregation-based model reduction method for nonlinear models of multi-zone building thermal dynamics. The full-order model, which is already a lumped-parameter approximation, quickly grows in state space dimension as the number of zones increases. An advantage of the proposed method, apart from being applicable to the nonlinear thermal models, is that the reduced model obtained has the same structure and physical intuition as the original model. The key to the methodology is an analogy between a continuous-time Markov chain and the linear part of the thermal dynamics. A recently developed aggregation-based method of Markov chains is employed to aggregate the large state space of the full-order model into a smaller one. Simulations are provided to illustrate tradeoffs between modeling error and computation time. (C) 2014 Elsevier Ltd. All rights reserved.