In this study we investigate the dynamics of the region that includes Greece, the Aegean Sea and Asia Minor. In a least-squares inversion we solve for a continuous strain rate field, and corresponding velocity field, that satisfies 872 GPS data. The estimate of the geodetic strain rate field provides constraints for our dynamic analysis. Next, we separately solve the depth integrated 3-D force balance equations for depth-integrated deviatoric stresses within the lithosphere, in which body force input comes from differences in vertically integrated vertical stress, or differences in gravitational potential energy per unit area (GPE). These GPE estimates calibrate the absolute magnitudes of deviatoric stresses that are acting within the lithosphere. Further, we investigate the sensitivity of our stress field solutions by using two different crustal structure models: one from compiled crustal structure estimates obtained primarily from relatively recent seismic observations and the other from the Crust 2.0 model. In an iterative least-squares inversion we then solve for stress field boundary conditions that, when added to the contribution of deviatoric stresses associated with GPE differences, provides a best fit to the directions of principal axes and relative magnitudes of the principal axes of the rates of strain obtained in the kinematic analysis. Robust features that arise from the boundary condition solution are NNE forcing along the southern boundary east of about 33 degrees E (0.5-1.2 x 10(12) N m(-1)), with a rapid anticlockwise rotation of forces to the west of this, along with an outward pulling force (similar to 0.4 x 10(12) N m(-1)) directed SSW along the entire Hellenic Arc segment. This force system along the Hellenic Arc can be interpreted as a result of slab rollback. The total depth integrated 3-D deviatoric stresses in the final dynamic solution provides an excellent match to the deformation indicators throughout the region, with vertically integrated stress magnitudes of order 0.5-2.5 x 10(12) Nm(-1). We use constraints from derived stress magnitudes, together with GPS-defined scalar values of strain rate magnitude, to define bulk effective viscosities of the lithosphere. Depth-averaged effective viscosities for the entire lithosphere are high within the Black Sea, of order 0.7-3 x 10(23) Pa-s, relative to surrounding continental lithosphere. North Anatolian shear zone, northern Aegean Sea and Gulf of Corinth are characterized by low depth averaged viscosities of order 1-5 x 10(21) Pa-s. Deviatoric stresses from GPE differences and boundary condition effects combine in surprising ways in some regions, resulting in near total stress cancellation in areas such as the southern Aegean Sea and portions of the central Anatolian block. GPE differences combine with boundary condition effects along the eastern segment of the North Anatolian Fault (NAF) in a way that is compatible with the hypothesis that motion on the NAF was facilitated by slab detachment beneath East Anatolia and dynamic uplift of East Anatolian Plateau. In general, GPE differences play a nearly equal role as boundary condition influences in their contribution to the total deviatoric stress field. The low depth integrated deviatoric stress magnitudes throughout the region suggest that zones of active deformation are facilitated by dramatic weakening mechanisms throughout the lithosphere.