A statistical meso-damage thermo-hydro-mechanical (T-H-M-D) coupling model is developed to investigate the trans-scale progressive cracking progress of the shale reservoir stimulated by liquid nitrogen (LN2). The theoretical framework encompasses: (1) local thermal equilibrium theory to formulate the heat exchange between the rock matrix and the cold fluid; (2) compressible fluid seepage theory to simulate the mass transport in pore media; (3) classical thermo-poroelatic theory to analyze the matrix deformation characteristics; (4) elastic-brittle-damage theory with double damage functions to identify the shale failure process at low temperatures; (5) two independent statistical distribution functions to describe the shale meso-scale heterogeneity in both thermos-mechanical and hydraulic properties; (6) Helmholtz free energy functions to calculate temperature- and pressure-dependent fluid thermodynamic properties; (7) experimental and empirical formulas to characterize temperature- and damage-related rock properties. The numerical model is implemented into the finite element method to analyze the underlying response mechanisms of the shale reservoir interacting with LN2. Simulation results show that LN2-shale interaction can induce multiple initiation positions, resulting in multiple asymmetric macro-cracks. The failure mechanism of the shale during cryogenic fracturing is dominated by elastic-brittle tensile-induced damage. Cryogenic shock to enhance shale elastic modulus and induce thermal tensile stress is beneficial to forming complex fracture networks. It is also found both fracture initiation pressure and induced fracture morphology are more sensitive to fluid injection temperature relative to reservoir temperature.