## 4 Numerical Results and Discussion
### Flow Behaviours in DFNs  
Nonlinear flow behavior exists not only in a single fracture or a small number of fractures but also in DFNs. Experimental and numerical research on a single fracture can also be applied to the complex fracture network. A fracture network is generated according to the method described in Section 3.2, as shown in Figure 14. The side length of the DFN model is 50 m. As the permeability of the rock matrix is too low, it is considered that the fluid only flows in the fractures. The horizontal in-situ stress \(\sigma_{\rm x}\) and vertical in-situ stress \(\sigma_{\rm y}\) are applied to the model. The in-situ stress values are 1 MPa, 3 MPa, 5 MPa, 7 MPa, and 10 MPa, respectively, and the lateral pressure coefficient is 1. The upper boundary of the model is the hydraulic boundary. The water pressure at the lower boundary is 0, and the left and right boundaries are impermeable boundaries. The range of the hydraulic gradient is 10\({}^{-5}\)\(\sim\)10\({}^{0}\). The initial fracture aperture values considered are 1 mm, 1.5 mm, and 2 mm. The statistical distribution parameters of the two groups of fractures and the material parameter values of the matrix and fracture are shown in Tables 5 and 6, respectively. The material parameters are assigned to CPE6 element for the rock matrix and CPE6 MP element for the fracture, respectively, before solving by ABAQUS. \begin{table}
\begin{tabular}{c c c c} \hline \hline
**Cases** & **Inlet (s)** & **Outlet (s)** & **Remarks** \\ \hline Case 1 & 1 & 2 & One inlet and one outlet \\ Case 2 & 1 & 1,2 & One inlet and two outlets \\ Case 3 & 1,2 & 2 & Two inlets and one outlet \\ Case 4 & 1,2 & 1,2 & Two inlets and two outlets \\ \hline \hline \end{tabular}
\end{table} 
Table 4: Simulated cases. Figure 13: Comparison of experimental results with the numerical results of Model 2. Figure 15 shows that there is a nonlinear relationship between the hydraulic gradient and discharge. With an increase in the hydraulic gradient, the nonlinear phenomenon becomes more obvious, and this nonlinear relationship can be described well by the Forchheimer law. The numerical results show that the curvature of the curve increases with the hydraulic gradient, which indicates that the nonlinearity of Q and J becomes stronger under a higher flow velocity and the correlation coefficient R\({}^{2}\) > 0.99. Figure 15 shows that a larger in-situ stress induces smaller discharge under the same hydraulic gradient. Since the fracture would close under the influence of the normal stress, the larger the normal stress is, the smaller the equivalent hydraulic aperture will be, which induces the lower permeability. [PERSON] et al. [51] uses [PERSON]'s simulation method to simulate the nonlinear flow of fracture surfaces with different roughnesses. This nonlinear characteristic is consistent with the simulation results of this study; this view was also confirmed by [PERSON] et al. [50]. In addition, the influence of stress on aperture was considered in this paper. The relationship between the transmissivity (the ratio of the discharge to hydraulic gradient) and hydraulic gradient of the model is shown in Figure 16. The curve can be roughly divided into three stages. In the first stage, when the hydraulic gradient is less than \(10^{-3}\), the flow velocity is low in a fracture, and the effect of inertia can be ignored. The transmissivity is constant, which indicates that the fluid flow is linear and that the cubic law is applicable. In the second stage, when the hydraulic gradient is less than \(10^{-2}\) but greater than \(10^{-3}\), the transmissivity changes slowly with the hydraulic gradient, indicating that there is a weak inertial region in this stage. In the third stage, when the hydraulic gradient is greater than \(10^{-2}\), the flow velocity is high, and the effect of inertia plays a leading role,