A new forward modeling approach to simulate the extension/closure and orientation statistics of evolving fracture networks is presented. The model is fully dynamical and couples fracturing to other processes. Thus, fracturing affects hydrology and, in turn, its development is affected by fluid pressure. In this way, highly pressured fluids can enhance their own migration while low-pressured ones may become trapped. Fracturing affects the stress tensor through rock volumetric changes and stress affects fracture dynamics. Thus, to predict fracture dynamics, we coevolve a set of fracture variables simultaneously with mass, momentum, and energy conservation equations for the solid and multiple fluid phases. In our computational approach, a representative set of putative fractures of a range of orientations is introduced. The time-dependent properties of each realized fracture is calculated by using the stress component normal to its fracture plane, pressure, and rock properties. The volumetric strain caused by fracturing allows for a self-limiting feedback that we account for using a nonlinear incremental stress rheology. The anisotropic fracture permeability is obtained using the predicted fracture network statistics. Thus, the coupling between rock deformation, notably fracturing, and hydrology is accounted for as is that between fracturing and stress. The statistical aspect of fracture network dynamics is described by assuming a probability distribution characterizing variations in rock strength within a nominally uniform lithology. The dependence of fracture density and length on the rate of fluid pressure or stress variation is thereby captured. Embedding the model in a three-dimensional basin finite element simulator, we illustrate the dynamical nature of the location and character of fracture zones in a sedimentary basin.