Theexactwithinpotentialflowintegral equation approach of Evans and Ford [Proc. R. Soc. London Ser. A452, 373 (1996)] for thenormalsolitary wave, is here generalized to 2‐layer, ‘‘internal’’ solitary waves. This differs in its mathematical form from other exact integral equation methods based on the complex velocity potential. For both ‘‘rigidlid’’ (i.e., flat toplayer surface) and ‘‘free‐surface’’ boundary conditions, a set of coupled non‐linear integral equations are derived by an application of Green’s theorem. For each point on the layer interface(s), these describe functional constraints on the profiles and interface fluid velocity moduli; theexactprofiles and velocities being those forms that satisfy these constraints atallsuch interface points. Using suitable parametric representations of the profiles and interface velocity moduli as functions of horizontal distance,x, and utilizingtailoredquadraturemethods [Int. J. Comput. Math. B6, 219 (1977)], numerical solutions were obtained by the Newton–Raphson method that are highly accurate even atlargeamplitudes. For ‘‘rigid lid’’ boundary conditions, internal wave solutions are presented for layer density and depth ratios typical ofoceanicinternal wave phenomena as found in the Earth’s marginal seas. Their various properties, i.e., mass, momentum, energy, circulation, phase and fluid velocities, streamline profiles, internal pressures, etc., are evaluated and compared, where possible, with observed properties of such phenomena as reported, for example, from the Andaman Sea. The nature of the limiting (or ‘‘maximum’’) internal wave is investigated asymptotically and argued to be consistent with two ‘‘surge’’ regions separating the outskirts flow from a wide mid‐section region of uniform ‘‘conjugateflow’’ as advocated by Turner and Vanden‐Broeck [Phys. Fluids31, 286 (1988)]. ©1996 American Institute of Physics.