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 condition. Unless explicitly accounted for during ISM initialization, the surface mass balance forcing from the ESM is far from equilibrated with the ISM flux divergence. Upon coupling to ESM forcing, the result is a shock and unphysical and un- desirable transients in ice geometry and other state variables. Under the assumption of thermomechanical quasi-equilibrium, we present an optimization approach for obtaining optimal ISM initial conditions when coupling to ESMs. This approach targets finding an ISM initial condition characterized by the basal sliding coefficient and the basal topography fields that balances a best fit to surface velocity observations against the minimization of undesirable and unphysical transients when coupling to ESM forcing for forward model simulations. A quasi-Newton method is used to solve the resulting large-scale, PDE-constrained optimization problem, where the 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 to the initialization 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.