Multiphase flow and transport phenomena within fractures are important because fractures often represent primary flow conduits in otherwise low permeability rock. Flows within the fracture, between the fracture and the adjacent matrix, and through the pore space within the matrix typically happen on different length and time scales. Capturing these scales experimentally is difficult. It is therefore useful to have a computational tool that establishes the exact position and shape of fluid/fluid interfaces in realistic fracture geometries. The level set method is such a tool. Our progressive quasistatic (PQS) algorithm based on the level set method finds detailed, pore-level fluid configurations satisfying the Young-Laplace equation at a series of prescribed capillary pressures. The fluid volumes, contact areas and interface curvatures are readily extracted from the configurations. The method automatically handles topological changes of the fluid volumes as capillary pressure varies. It also accommodates arbitrarily complicated shapes of solid confining surfaces. Here we apply the PQS method to analytically defined fracture faces and aperture distributions, to geometries of fractures obtained from high-resolution images of real rocks, and to idealized fractures connected to a porous matrix. We also explicitly model a fracture filled with proppant, using a cooperative rearrangement algorithm to construct the proppant bed and the surrounding matrix. We focus on interface movement between matrix and fracture, and snap-off of non-wetting phase into the fracture during imbibition in particular. The configuration of fluids is strongly affected by asperities in unpropped fractures and by the locally open regions at the proppant/matrix interface. The area of phase in contact with the matrix is nonlinear with phase saturation and strongly hysteretic, and thus transfer functions based on saturations should be used with caution. The effect of coupling fracture capillarity and matrix capillarity on multiphase flow properties depends on the relative sizes of typical pore throats in the matrix and typical aperture in the fracture. The simulations agree with direct obsevations of fluid configurations in fractures.