  The paper develops an unfitted finite element method for solving the Darcy
system of equations posed in a network of fractures embedded in a porous
matrix. The approach builds on the Hughes--Masud stabilized formulation of the
Darcy problem and the trace finite element method. The system of fractures is
allowed to cut through the background mesh in an arbitrary way. Moreover, the
fractures are not triangulated in the common sense and the junctions of
fractures are not fitted by the mesh. To couple the flow variables at multiple
fracture junctions, we extend the Hughes--Masud formulation by including
penalty terms to handle interface conditions. One observation made here is that
by over-penalizing the pressure continuity interface condition one can avoid
including additional jump terms along the fracture junctions. This simplifies
the formulation while ensuring the optimal convergence order of the method. The
application of the trace finite element allows to treat both planar and
curvilinear fractures with the same ease. The paper presents convergence
analysis and assesses the performance of the method in a series of numerical
experiments. For the background mesh we use an octree grid with cubic cells.
The flow in the fracture can be easily coupled with the flow in matrix, but we
do not pursue the topic of discretizing such coupled system here.
