Material models that depend on multiple independent variables are often necessary for accurate continuum numerical simulations, particularly for applications that involve large stresses and deformations. It is rare that purely physics-based models are used in continuumm simulations because of the attendant computational cost. Instead, experimental and microscale simulation data are expressed as phenomenological models and fed into continuum simulations. As the number of independent variables increases (e.g., internal state variables (ISVs)), such models are not only difficult to design but also need exponentially larger amounts of data to parameterize accurately. In this talk we examine an alternate procedure for developing multi-variable phenomenological constitutive models via multi-layer perceptron neural networks. The method is applied to experimental crush-curves and bulk modulus data for dry concrete sand. We find that reasonable network topologies and selected optimization parameters can be discovered by a factorial design of experiments procedure. Networks containing rectified linear unit (ReLU) activation functions are observed to produce excellent fits to the input data as well as acceptable generalization, provided derivatives of the neural network models are not required in the continuum simulations.