We present new results on robust model reduction for partial differential equations. Our contribution is threefold: 1.) The stabilization is achieved via closure models for reduced order models (ROMs), where we use Lyapunov robust control theory to design a new stabilizing closure model that is robust with respect to parametric uncertainties; 2.) The free parameters in the proposed ROM stabilization method are autotuned using a data-driven multi-parametric extremum seeking (MES) optimization algorithm; and 3.) The challenging 3D Boussinesq equation numerical test-bed is used to demonstrate the advantages of the proposed method.