The diffusion of atomic hydrogen through metals is governed by the hops of such atoms among interstitial lattice sites and their chances of remaining trapped in dislocations, impurities, grain boundaries, etc. These complex mechanisms need to be understood and modeled in order to study their effects on the ductility and fracture resistance of metals and motivates the development of numerical methods. When focusing on hydrogen transport in metals, the main difficulties are the treatment of the boundary conditions and distinguishing the evolution of the free and trapped hydrogen atoms. In this talk, we will present a new formulation that models transient solutions to this problem by employing an incremental variational principle whose stationarity conditions include the thermodynamic equilibrium between the two types of hydrogen atoms. This formulation ensures, among other things, the symmetry of the tangent operator leading to large computational savings. Numerical examples of stationary and transient solutions obtained with this method will be presented.