Numerical and physical modeling are the two main tools available for predicting the influence of water waves on coastal or offshore structures. Both models have their strengths and weaknesses. An integrated use of numerical and physical modeling which exploits their advantages can provide an optimal description of full-scale, realistic engineering problems. In this series of two papers, we report on a generalized three-dimensional (3D) deterministic coupling theory for multidirectional nonlinear wave propagation from a numerical model to a physical model up to second order. In this work, the second-order coupling theory developed by [Z. Yang, S. Liu, H.B. Bingham and J. Li, 2014a. Second-order coupling of numerical and physical wave tanks for 2D irregular waves. Part I Formulation, implementation and numerical properties. Coast. Eng. 92, 48–60] and the ad hoc unified wave generation theory developed by [H. Zhang, H.A., Schäffer, K.P., Jakobsen, 2007. Deterministic combination of numerical and physical coastal wave models. Coast. Eng. 54, 171–186.] have been extended to include the coupling between multidirectional nonlinear waves. In part I of this article series, the full second-order 3D coupling theory for multidirectional nonlinear waves is been described in detail. A novel generalized fully-nonlinear motion boundary equation has been derived, which allows the interface between the numerical and physical wave domains to be a 3D arbitrarily shaped wavemaker system. The corresponding 3D wave coupling equation is given for I-, L-, and O-shaped wavemaker layouts. The new formulation is presented in a unified form for an arbitrarily shaped layout of planar wavemakers, and covers both hinged and piston wavemaker types. For the second-order dispersive correction of the paddle stroke, the super-harmonic and sub-harmonic wave-wave interactions have been taken into account. For practical implementation, a typical discretization method for the 3D coupling equations is derived by considering the case of the I-shaped piston wavemaker. The new 3D coupling equations are solved by a combined five-point Lagrange interpolation and the fourth-order Runge–Kutta scheme. Numerical evaluations of the implementation have been conducted by considering a theoretical second-order unidirectional wave field over a range of spatial and time resolutions. The results thus obtained indicate that the proposed discretization scheme is accurate and effective. In addition, the precision of the discrete scheme is observed to be closely related to the spatial and temporal resolution. A separate experimental validation of the theory is presented in Part II.