We propose a novel approach to model viscoelasticity materials using neural networks, which capture rate-dependent and nonlinear constitutive relations. However, inputs and outputs of the neural networks are not directly observable, and therefore common training techniques with input-output pairs for the neural networks are inapplicable... To that end, we develop a novel computational approach to both calibrate parametric and learn neural-network-based constitutive relations of viscoelasticity materials from indirect displacement data in the context of multi-physics interactions. We show that limited displacement data hold sufficient information to quantify the viscoelasticity behavior. We formulate the inverse computation---modeling viscoelasticity properties from observed displacement data---as a PDE-constrained optimization problem and minimize the error functional using a gradient-based optimization method. The gradients are computed by a combination of automatic differentiation and physics constrained learning. The effectiveness of our method is demonstrated through numerous benchmark problems in geomechanics and porous media transport. (read more)