Fractures and faults play an important role in controlling the flow and transport properties in a reservoir and that is the main reason for which their characterization is very important in hydrocarbon exploration. In addition to this, the possibility of characterizing fractures can represent a great advantage in other fields, like geothermal exploration, hydrofracturing applications or volcanic risk evaluation. Seismic simulation by finite-difference modeling has been implemented as a tool to characterize fracture networks by testing their seismic signatures. The medium in which fractures are placed is represented by a squared model with the sides measuring 5000m; two receivers arrays, each including 1000 receivers, have been implemented: the first at the top of the model (to record the reflected wavefield), the second at a depth of 4000m (to record the transmitted wavefield). Fractures are randomly positioned in the center of the model, between 1250m and 3750m. In the case of a layered medium, only the middle layer, 1500m thick, contains random fractures. The sensitivity of seismic wave propagation to fractures has been analyzed by testing the influence of different fracture features: length, orientation, density (number of fractures). Even a single fracture affects the incident seismic signal, producing diffracted/scattered and transmitted waves. Moreover, its orientation significantly affects the reflected wavefield in the time domain in terms of amplitude and complexity of the response and in the frequency domain as well, where peak amplitude and peak frequency change depending on the fractures orientation. More than the orientation, the fracture length strongly affects the seismic signal in time and frequency domains. Considering a network of fractures, the imprint of the fracture orientation on the reflected wavefield is significant only in the frequency-wavenumber domain; on the other hand, it is much stronger for the transmitted wavefield in the frequency domain, where the peak amplitude and the peak frequency undergo high variation: in particular, the horizontal fractures produce the strongest frequencies attenuation and the lowest peak frequency. The fracture length variation produces the most interesting signature in the time domain for the reflected wavefield (an increase in the fracture length produces longer coda waves) and in the frequency domain for the transmitted field: in particular, the longest fractures produce the strongest frequencies attenuation and the lowest peak frequency. Significant is the signature of the fracture density, but it is particularly strong on the transmitted wavefield in the frequency domain: in particular, the highest number of fractures produce the same effect of the horizontal and longest fractures (strongest frequencies attenuation and lowest peak frequency). Eventually, the introduction of a fractured layer yields strong change in the incident signal, which is highly attenuated and disrupted. The results of the present work can have implications in all the fields where the detection and the characterization of fractures in the subsurface is vital, such as the geothermal exploration or the hydraulic fracturing applications. Despite of the complexity of a fractured systems, some recurrent trends and characteristic responses have been determined. However, it should kept in mind that fracture features are strictly related (they influence each other) and that different features can yield similar responses: this means that fracture characteristics can not be straightforward inferred from their seismic signatures.