A 2-D fully nonlinear numerical wave tank (NWT) is developed to investigate the Bragg resonance scattered by submerged multi-array breakwaters for random waves. This model is based on a boundary integral equation method with linear element scheme. The fully nonlinear free surface boundary condition is treated using the Mixed Eulerian-Lagrangian method and the 4th-order Runge-Kutta method. The incident random waves are generated by JONSWAP spectrum at one end of the wave tank. Two damping zones are deployed at both ends of the NWT to absorb the energy of the reflected and transmitted waves. In the regular wave cases, the results of Bragg reflection calculated are in good agreement with that of other experiments and numerical models. In addition, the simulated spectrum of random waves is also verified by the original input spectrum. The results of the random waves have the same trend as those of the Regular waves. The reflection coefficient for random waves at the first peak of resonance is about 70 percent of that of the regular wave, but the frequency of band width of Bragg effect has become wider and this advantage may compensate the peak reduction. Finally, we may conclude that the present model is adequate to use as a tool for coastal protection. Systematic studies for random waves propagating over series submerged breakwaters are conducted. The Bragg reflection will be enhanced with the increase of relative height, the length of bars, the number of breakwaters, and the toe angle of submerged breakwaters. In this study, it also reveals that the frequency of peak reflection for higher breakwaters has down shift phenomenon.