We address complications in the coupling of a dynamic ice sheet model (ISM) and forcing from an Earth System Model (ESM), which arise because of the unknown ISM initial conditions. Unless explicitly accounted for during ISM initialization, the ice sheet is far from thermomechanical equilibrium with the surface mass balance forcing from the ESM. Upon coupling to ESM forcing, the result is a shock and unphysical and undesirable transients in ice geometry and other state variables. Under the assumption of thermomechanical equilibrium, we present an approach for finding ISM initial conditions - characterized by optimization of the basal sliding coefficient and basal topography fields - that balance a best fit to surface velocity and basal topography observations against the minimization of unphysical transients when coupling to surface mass balance forcing. A quasi-Newton method isused to solve the resulting large-scale, PDE-constrained optimization problem, where the cost function gradients with respect to the parameter fields are computed using adjoints. After studying properties of our approach on a synthetic test problem, we apply the method towards obtaining optimal initial conditions for a model of the Greenland ice sheet.  Our results show that, in the presence of uncertainties in the basal topography, ice thickness should also be treated as an optimization variable. While the focus here is on the coupling between an ISM and ESM-derived surface mass balance, the method is easily extended to include optimal coupling to forcing from an ocean model through submarine melt rates.