In this paper, we consider the problem of finding surrogate models for large-scale second-order linear time-invariant systems with inhomogeneous initial conditions. For this class of systems, the superposition principle allows us to decompose the system behavior into three independent components. The first behavior corresponds to the transfer between the input and output having zero initial conditions. In contrast, the other two correspond to the transfer between the initial position or the initial velocity and the output when no input is applied. Based on this superposition of systems, our goal is to propose model reduction schemes that allow to preserve the second-order structure in the surrogate models. To this aim, we introduce tailored second-order Gramians for each system component and compute them numerically, solving Lyapunov equations. As a consequence, two methodologies are proposed. The first one consists in reducing each of the components independently using a suitable balanced truncation procedure. The sum of these reduced systems provides an approximation of the original system. This methodology allows flexibility on the order of the reduced-order model. The second proposed methodology consists in extracting the dominant subspaces from the sum of Gramians to build the projection matrices leading to a surrogate model. Additionally, we discuss error bounds for the overall output approximation. Finally, the proposed methods are illustrated by means of benchmark problems. [ABSTRACT FROM AUTHOR]