This paper presents a time-domain modeling and shape optimization framework for microwave structures including metasurfaces, based on a nodal Discontinuous Galerkin Time-Domain (DGTD) method. In particular, we employ an unstructured mesh and the inherent mesh refinement ability of DGTD to model multi-scale geometries via adaptive mesh refinement. More importantly, we integrate a multi-tier local time-stepping technique into the time integration of DGTD, which significantly alleviates the cost introduced by mesh refinement. We further present a flexible parameterization approach by defining the contours of metasurface unit cells by B-spline curves, which allows us to explore a wide range of smooth shapes by only a few design variables. Besides, we adopt a Polynomial Chaos-Kriging (PCK) surrogate method to approximate the full-wave DGTD model, reducing the computational cost for optimization problems that require time-consuming simulations by three orders of magnitude. The B-spline parameterization and the PCK surrogate model are combined with a pattern search-based Pareto front algorithm to optimize metasurfaces for multiple design objectives. The proposed optimization framework is also applicable to other microwave structures. We demonstrate the effectiveness of our approach through the optimization of an omega-bianisotropic Huygens’ metasurface unit cell and an E-plane microwave T-junction.