Wind turbine controllers are used to optimize the performance of wind turbines such as to reduce power variations and fatigue and extreme loads on wind turbine components. Accurate tuning and design of modern controllers must be done using low-order models that accurately captures the aeroelastic response of the wind turbine. The purpose of this thesis is to investigate the necessary model complexity required in aeroelastic models used for controller design and to analyze and propose methods to design low-order aeroelastic wind turbine models that are suited for model-based control design. The thesis contains a characterization of the dynamics that influence the open-loop aeroelastic frequency response of a modern wind turbine, based on a high-order aeroelastic wind turbine model. One main finding is that the transfer function from collective pitch to generator speed is affected by two low-frequency non-minimum phase zeros. To correctly predict the non-minimum phase zeros, it is shown to be essential to include lateral tower and blade flap degrees of freedom. The thesis describes and analyzes various methods to design low-order aeroelastic models of wind turbines. Low-order models are designed by modal truncation by using the aeroelastic mode shapes of a fully flexible wind turbine. To capture the effect of shed vorticity and dynamic stall, a relatively large number of aerodynamically dominated modes are required, due to the assumption of independent annular flow tubes in the Blade Element Momentum theory (BEM). A set of accurate reduced-order models are subsequently designed assuming quasi-steady aerodynamics, by truncation with a set of low-frequency mode shapes. In a comparison, the balanced truncation method is found to be able to capture the effect of the shed vorticity and dynamic stall using only few states. A set of reduced-order models obtained at various operating points are shown to be easily connected by interpolation and are thereby suited for gain-scheduling control design. A new method is proposed to reduce separately the number of structural and aerodynamic states in aeroelastic models by using a set of structural and aerodynamic basis functions. Accurate approximation of the low-frequency blade response is obtained using a set of purely structural blade mode shapes and by static residualization of the high-frequency blade modes. The effect of shed vorticity and dynamic stall on the blade response can be captured using a set of aerodynamic slave modes of the low-frequency structural flap modes. Future work is to test the proposed method on others models of the unsteady aerodynamic forces on wind turbines and to use the reduced models for controller design.