Accurate reduced-order models of turbulent flows have been traditionally constructed with the proper orthogonal decomposition (POD), however the method has been limited to prototypical flows over a narrow parameter range. An orthogonal complement subspace method is developed here to treat inhomogeneous boundary conditions while implicitly coupling the velocity and temperature fields of turbulent convection. A new flux matching procedure is formulated as a state space residual series expansion to efficiently model parameter dependent convection, greatly extending the utility of the reduced-order modeling framework. An illustrative test case of turbulent channel flow over heated blocks shows flow and thermal models with 95% accuracy over the domain can be produced, while simultaneously reducing the number of degrees of freedom by a factor of 10(3). Error bounds are formulated and provide a posteriori error estimates for the reduced-order model. (c) 2006 Published by Elsevier Ltd.