Numerical simulations presented in this work were performed on a single-processor workstation (2.5GHz, 1Gb RAM). We performed a direct numerical simulation, integrating the dynamical equations of motion (4) and (5) using pseudo-spectral method with resolution of wavenumbers. Numerical integrator used for advancing in time was RK7(8) presented in [31]. Time step was where is the period of the shortest wave on the axis. Approximate processor time for this work was 4.5 weeks.

In our numerical experiment, we force the system in the -space ring with and . This ring is located at the low wavenumber part of the -space in order to generate energy cascade toward large 's, but we deliberately avoid forcing even longer waves () because our experience shows that this would lead to undesirable strong anisotropic effects. In the ring, we fix the shape to coincide with the ZF spectrum, , and hence we set , . These fixed amplitudes were then multiplied by random phase factors. Thus, surface and velocity potential were set to , where were uniformly distributed in and

where for the surface and for velocity potential. Damping was applied in both small and large wavenumber regions. At small wavenumbers inside the forcing ring, we applied an adaptive damping to prevent formation of undesirable ``condensate'' which could spoil isotropy and locality of scale interactions. At large wavenumbers, damping is needed to absorb the energy cascade and, therefore, to avoid ``bottleneck'' spectrum accumulation near the cutoff wavenumber. In our simulations, we implemented the damping as a low-pass filter applied to the -space variables at each time step. The damping function had the form

Nonlinearity parameter was set to , which is a sufficient value to produce a resonance broadening for supporting energy cascade.