We describe an alternative procedure for obtaining approximate Brueckner orbitals in ab initio electronic structure theory. Whereas approximate Brueckner orbitals have traditionally been obtained by mixing the orbitals until the coefficients of singly substituted determinants in the many-electron wavefunction become zero, we remove singly substituted determinants at the outset and obtain orbitals which minimize the total electronic energy. Such orbitals may be described as variational Brueckner orbitals. These two procedures yield the same set of exact Brueckner orbitals in the full configuration interaction limit but differ for truncated wavefunctions. We consider the simplest variant of this approach in the context of coupled-cluster theory, optimizing orbitals for the coupled-cluster doubles (CCD) model. An efficient new method is presented for solving the coupled equations defining the energy, doubles amplitudes, and orbital mixing parameters. Results for several small molecules indicate nearly identical performance between the traditional Brueckner CCD method and the variational Brueckner orbital CCD approach. However, variational Brueckner orbitals offer certain advantages: they simplify analytic gradients by removing the need to solve the coupled-perturbed Brueckner coupled-cluster equations for the orbital response, and their straightforward extensions for inactive orbitals suggests possible uses in size-extensive models of nondynamical electron correlation. Application to O4+ demonstrates the utility of variational Brueckner orbitals in symmetry breaking cases.