A variable-fidelity method can remarkably improve the efficiency of a design optimiza-tion based on a high-fidelity and expensive numerical simulation, with assistance of lower-fidelity and cheaper simulation(s). However, most existing works only incorporatetwo"levels of fidelity, and thus efficiency improvement is very limited. In order to reduce the number of high-fidelity sim-ulations as many as possible, there is a strong need to extend it to three or more fidelities. This arti-cle proposes a novel variable-fidelity optimization approach with application to aerodynamic design. Its key ingredient is the theory and algorithm of a Multi-level Hierarchical Kriging (MHK), which is referred to as a surrogate model that can incorporate simulation data with arbi-trary levels of fidelity. The high-fidelity model is defined as a CFD simulation using a fine grid and the lower-fidelity models are defined as the same CFD model but with coarser grids, which are determined through a grid convergence study. First, sampling shapes are selected for each level of fidelity via technique of Design of Experiments (DoE). Then, CFD simulations are conducted and the output data of varying fidelity is used to build initial MHK models for objective (e.g. CD) and constraint (e.g. CL, Cm) functions. Next, new samples are selected through infill-sampling criteria and the surrogate models are repetitively updated until a global optimum is found. The proposed method is validated by analytical test cases and applied to aerodynamic shape opti-mization of a NACA0012 airfoil and an ONERA M6 wing in transonic flows. The results confirm that the proposed method can significantly improve the optimization efficiency and apparently out-performs the existing single-fidelity or two-level-fidelity method.