Summary: We study the problem of determining through a variational data assimilation approach the initial condition for a coupled set of nonlinear partial differential equations from which a model trajectory emerges in agreement with a given set of time-distributed observations. The partial differential equations describe an idealized coupled atmospheric-ocean model on a rotating torus. The model consist of the viscous shallow-water equations in geophysical scaling that represents the large-scale atmospheric dynamics coupled via a simplified but physically plausible coupling to a model that represents the large-scale ocean dynamics and consists of the incompressible two-dimensional Navier-Stokes equations and an advection-diffusion equation. We propose a variational algorithm (4D-Var) of the coupled data assimilation problem that is solvable and computable. This algorithm relies on the use of a variational cost functional that is tailored to the regularity of the coupled model as well as to the regularity of the observations by means of derivative-based norms. We support this proposal by developing regularity results for an idealized coupled atmospheric-ocean model using the concept of classical solutions. Based on these results we formulate a suitable cost functional. For this cost functional we prove the existence of optimal initial conditions in the sense of minimizers of the cost functional and characterize the minimizers by a first-order necessary condition involving the coupled adjoint equations. We prove local convergence of a gradient-based descent algorithm to optimal initial conditions using second-order adjoint information. Instrumental for our results is the use of suitable Sobolev norms instead of the standard Lebesgue norms in the cost functional. The index of the actual Sobolev space provides an additional scale selective mechanisms in the variational algorithm.