Summary: Machine learning methods have in various ways emerged as a useful tool for modeling the dynamics of physical systems in the context of partial differential equations (PDEs). Nonlinear conservation laws (NCLs) of the form \(u_t + f(u)_x = 0\) play a vital role within the family of PDEs. A main challenge with NCLs is that solutions contain discontinuities. That is, one or several jumps of the form \((u_L(t), u_R(t))\) with \(u_L \neq u_R\) may move in space and time such that information about \(f(u)\) in the interval associated with this jump is not present in the observation data. Moreover, the lack of regularity in the solution \(u(x, t)\) prevents use of physics informed neural network (PINN) and similar methods. The purpose of this work is to propose a method to identify the nonlinear flux function \(f(u, \beta)\) with variable parameters of an unknown scalar conservation law \[ u_t + f(u, \beta)_x = 0\tag{\(\ast\)} \] with \(u\) as the dependent variable and \(\beta\) as the parameter. In a recent work we introduced a framework coined ConsLaw-Net that combines a symbolic multi-layer neural network and an entropy-satisfying discrete scheme to learn the nonlinear, unknown flux function \(f(u; \beta)\) for various fixed \(\beta\). Learning the flux function with variable parameters, marked as \(f(u, \beta)\), is more challenging since it requires an understanding of the relationship between the variable \(u\) and parameter \(\beta\) as well. In this work we demonstrate how to couple ConsLaw-Net to the Linear Regression Neural Network (LRNN) to learn the functional form of the two variable function \(f(u, \beta)\). In addition, ConsLaw-Net is here further developed and made more generic by using a refined discrete scheme combined with a more general symbolic neural network (S-Net). We experimentally demonstrate that the upgraded ConsLaw-Net integrated with LRNN is well-suited for learning tasks, achieving more accurate identification than existing learning approaches when applied to problems that involve general classes of flux functions \(f(u, \beta)\). The investigations of this work are restricted to the class of polynomial, rational functions.