This thesis presents an energy-conservative data-driven approach in modelling the closure terms of the Navier-Stokes equations casted through the Variational Multiscale (VMS) framework. For context, the VMS framework is applied in designing stabilised finite element methods for multiscale phenomena in which stability is not guaranteed. Under standard applications of this framework with appropriate cubature rules and iteration routines, the only error source is model related in the assumptions for the fine scales, thereby influencing the closure terms of the coarse-scale equations. This issue can be rectified by predicting the perfect closure terms with the application of neural networks. To maintain energy conservation (as a generality, not as a guarantee), a constrained objective function was formulated, which was then translated into an unconstrained problem suitable for present day neural network training tools via the Augmented Lagrangian method. The features of the eight neural networks used were extracted from two spatiotemporal stencils. A priori evaluation of the fine-scale model demonstrated a tendency to overpredict closure terms for the Re_tau = 180 turbulent channel flow test case, yielding desired energy dissipative behaviour. This was true for both regimes of the flow in which the closure terms yield energy gain and energy loss via backscatter and cascade respectively. Unverified a posteriori simulation of the test case with the model displayed excessively high energy dissipation, in addition to high frequency oscillations, which were undampened by the time marching method used.