In this work we formulate an efficient method for the description of fully many-body localized systems in weak contact with thermal environments at temperature T. The key idea is to exploit the representation of the system in terms of quasi-local integrals of motion (l-bits) to efficiently derive the generator for the quantum master equation in Born-Markov approximation. We, moreover, show how to compute the steady state of this equation efficiently by using quantum-jump Monte-Carlo techniques as well as by deriving approximate kinetic equations of motion. As an example, we consider a one-dimensional disordered extended Hubbard model for spinless fermions, for which we derive the l-bit representation approximately by employing a recently proposed method valid in the limit of strong disorder and weak interactions. Coupling the system to a global thermal bath, we study the transport between two leads with different chemical potentials at both of its ends. We find that the temperature-dependent current is captured by an interaction-dependent version of Mott's law for variable range hopping, where transport is enhanced/lowered depending on whether the interactions are attractive or repulsive, respectively. We interpret these results in terms of spatio-energetic correlations between the l-bits.