AbstractThe present paper is concerned with the numerical approximation of a three‐parameter family of Boussinesq systems. The systems have been proposed as models of the propagation of long internal waves along the interface of a two‐layer system of fluids with rigid‐lid condition for the upper layer and under a Boussinesq regime for the flow in both layers. We first present some theoretical properties of the systems on well‐posedness, conservation laws, Hamiltonian structure, and solitary‐wave solutions, using the results for analogous models for surface wave propagation. Then the corresponding periodic initial‐value problem is discretized in space by the spectral Fourier Galerkin method and for each system, error estimates for the semidiscrete approximation are proved. The spectral semidiscretizations are numerically integrated in time by a fourth‐order Runge–Kutta‐composition method based on the implicit midpoint rule. Numerical experiments illustrate the accuracy of the fully discrete scheme, in particular its ability to simulate accurately solitary‐wave solutions of the systems.