A Boussinesq-type model is derived which is accurate to O(kh)4 and which retains the full representation of the fluid kinematics in nonlinear surface boundary condition terms, by not assuming weak nonlinearity. The model is derived for a horizontal bottom, and is based explicitly on a fourth-order polynomial representation of the vertical dependence of the velocity potential. In order to achieve a (4,4) Padé representation of the dispersion relationship, a new dependent variable is defined as a weighted average of the velocity potential at two distinct water depths. The representation of internal kinematics is greatly improved over existing O(kh)2 approximations, especially in the intermediate to deep water range. The model equations are first examined for their ability to represent weakly nonlinear wave evolution in intermediate depth. Using a Stokes-like expansion in powers of wave amplitude over water depth, we examine the bound second harmonics in a random sea as well as nonlinear dispersion and stability effects in the nonlinear Schrödinger equation for a narrow-banded sea state. We then examine numerical properties of solitary wave solutions in shallow water, and compare model performance to the full solution of Tanaka (1986) as well as the level 1, 2 and 3 solutions of Shields & Webster (1988).