Hydraulic fracturing or hydraulic stimulation is one of the most effective methods of enhancing hot dry rock (HDR) geothermal system productivity. In the present study, network models of "fractal geometry" approximate a 3D structure of fractured rocks. The fracture network models are generated by distributing fractures randomly in space and assuming the fractal equation correlating the number of fractures and fracture lengths. Based on this approach, a mathematical model of hydraulic rock fracturing is proposed. The model incorporates approximations of the fracture mechanical behavior drawn from the rock mechanics literature, a very simplified analysis of the operative physical processes, and mapping of the connectivity of fracture network to a cubic regular grid. It is assumed that the flow properties of a stochastic fracture network depend on the fluid pressure. Along with the fractal-type distribution of the fracture lengths, the fracture surfaces are also assumed to follow fractal geometry. The latter allows numerical simulation of the natural rock fracture dilation caused by fracture shear offset. On of the problems that can be resolved by fracture surface modeling is the apparent limitation on the number of fractures that can be analyzed experimentally. In this respect, the suggested mathematical model can be used to simulate fractal surfaces identical to fractures found in the natural rocks. Taken together, these approaches permit the approximate engineering resolution of the multi-parametric, highly complex mechanical problem.