Dynamic stress evolution during earthquake rupture contains information of fault frictional behavior that governs dynamic rupture propagation. Most of earthquake stress drop and evolution studies are based on kinematic slip inversions. Several dynamic inversion methods in the literature require dynamic rupture modeling that makes them cumbersome with limited applicability. In this study, we develop a fault-stress model of earthquake sources in the framework of the representation theorem. We then propose a dynamic stress inversion method based on the fault-stress model to directly invert for dynamic stress evolution process on the fault plane by fitting seismic data. In this inversion method, we calculate numerical Green’s function once only, using an explicit finite element method EQdyna with a unit change of shear or normal stress on each subfault patch. A linear least-squares procedure is used to invert for stress evolution history on the fault. To stabilize the inversion process, we apply several constraints including zero normal slip (no separation or penetration of the fault), non-negative shear slip, and moment constraint. The method performs well and reliably on a synthetic model, a checkerboard model and the 2016 M_w 5.0 Cushing (Oklahoma) earthquake. The proposed fault-stress model of earthquake sources with inversion techniques such as one presented in this study provides a new paradigm for earthquake source studies using seismic data, with a potential of deciphering more physics from seismic recordings of earthquakes.