In this paper we propose a (suboptimal) H2 model reduction method for a class of linear network systems that describe diffusively coupled networks. To preserve a network structure, we form a reduced-order model by using the characteristic matrix of a graph clustering so that the reduced-order model has less number of vertices. Then, we formulate the model reduction problem as a nonconvex optimization problem with binary variables, aiming for a graph clustering that minimizes the H2-norm of the approximation error. Based on the controllability and the observability Gramians of the error system we derive an optimization problem with mixed-binary variables and then we propose a convex relaxation of the binary variables, leading to a smooth opt...