The scale of large finite element (FE) models may nowadays easily exceed 10 7 or even 10 8 degrees of freedom (DOF), leading to excessive calculation times when performing transient simulations. Such long simulation times hinder effective structural or thermal design and optimization and make any engineering insight into a problem difficult. The Krylov subspace-based model order reduction (MOR) is a reduction method based on projection of a discretized model onto a lower dimension subspace. The paper presents a methodology based on this method in the context of thermal transient simulation of a large scale subsea equipment FE model. The finite element model mesh size exceeds 30 × 10 6 DOFs. The problem has nonzero initial conditions (ICs) and has to be transformed into a problem with zero ICs in order to apply the Krylov based MOR. Coupling the Krylov based MOR models employs a novel technique involving coupling through their surface interfaces. The approach is compared with the solution obtained using a full FE simulation which takes about 7 days to solve with a fine time step. The results are compared using an error norm which computes maximum absolute difference of temperature fields over time taking the full FE simulation with the fine time step as a reference. The study shows that applying the proposed method using Krylov MOR for performing thermal transient simulations is valid and leads to substantial reduction of the computational time.