Summary: In this work, an initial boundary value problem for a system of viscoelastic wave equations with nonlinear boundary source term of the form \[ \begin{gathered} (u_i)_{t t} - \Delta u_i - \Delta(u_i)_{t t} + \mathop{\int}\limits_0^t g_i(t - s) \Delta u_i(s) d s - \Delta(u_i)_t = 0,\quad \text{in} \quad \Omega \times(0, T), \\ u_i(x, 0) = \varphi_i(x),\quad (u_i)_t(x, 0) = \psi_i(x),\quad \text{in} \quad \Omega, \\ u_i(x, t) = 0,\quad \text{on} \quad \Gamma_0 \times(0, T), \\ \partial_\nu(u_i)_{t t} + \partial_\nu u_i - \mathop{\int}\limits_0^t g_i(t - s) \partial_\nu u_i(s) d s + \partial_\nu(u_i)_t + f_i(u) = 0,\quad \text{on} \quad \Gamma_1 \times(0, T), \end{gathered} \] where \(i = 1,\ldots, l\) (\(l \geq 2\)) is considered in a bounded domain \(\Omega\) in \(\mathbb{R}^N\) (\(N \geq 1\)). By the Faedo-Galerkin approximation method we obtain existence and uniqueness of weak solutions. Under appropriate assumptions on initial data and the relaxation functions, we establish general decay and blow up results associated to solution energy. Estimates for lifespan of solutions are also given.