Reliable estimation of geomechanical properties (i.e., Young's modulus and Poisson's ratio) of shale can provide important constraints to guide production by identifying sweet spots and locations for effective hydraulic fracturing. Amplitude Variation with Offset (AVO) inversion may make important contributions to such tasks, since it is an effective method to estimate elastic parameters of target layers. However, to estimate geomechanical and even elastic properties of shale from AVO inversion is not a trivial task, since organic-rich shale formations are often anisotropic. Our objective is to apply a nonlinear AVO inversion using the exact Zoeppritz solutions instead of its linear approximation. It allows estimation of seismic anisotropy and furthermore estimates of anisotropy of geomechanical properties. We first reformulate the exact Zoeppritz equations for reflection coefficients in terms of four parameters (one ratio of background P-wave and S-wave velocities, and three contrasts of P-wave and S-wave velocities, and density). An adjoint state technique is applied to compute the gradient of reflection amplitudes modeled by the parameters. This allows the nonlinear AVO inversion possible. We then propose a workflow to estimate seismic anisotropy and geomechanical properties of organic-rich shale. It is based on analyses of results from the AVO inversion. The anisotropy of the model shale is related to the kerogen volume fraction values using measured well logs and laboratory data for various shale formations. By applying inversion tests, we determine behaviors of the AVO inversion solutions developed for isotropic media when the target shale formation instead has seismic anisotropy related to organic content. These tests show that the inversion accurately determines horizontal P-wave and S-wave velocities and underestimates density when a far angle range is applied with input data. When the angle range is small, the inversion can obtain reliable vertical velocities, and correct density. Therefore, seismic anisotropy of the model can be estimated by comparing these inverted horizontal and vertical velocities. In addition, geomechanical properties of the model are also reliably determined in both horizontal and vertical directions. In contrast to most conventional AVO inversions based on linear approximations of the Zoeppritz equations, the proposed Zoeppritz AVO inversion is not limited by assumptions of weak contrasts and seismic isotropy. This allows better estimations of elastic and geomechanical properties and their anisotropy for unconventional shale play which are highly anisotropic and often surrounded by hard layers to generate strong contrasts. Given reliable inference of geomechanical properties from the AVO inversion, the results can directly impact to quantify fracability of unconventional play. Consequently, the workflow for nonlinear AVO inversion contributes to optimization of well placement, stimulated reservoir volume (SRV), and completion design of unconventional reservoir development.