Abstract: We present a workflow based on physics-informed neural networks (PINNs) to model multiphase flow in fractured porous media. After validating the method on a synthetic problem, it was applied to an experimental dataset where brine was injected into a CO$_2$-saturated fractured shale core, with the fractures and fluid distribution imaged using CT scans. A domain decomposition approach was used for matrix and fractures, with a multi-network architecture to separately compute the state variables of each domain. A novel hybrid pre-training strategy was applied to enable the model to capture the system's complex multiscale attributes. By history matching of multi-fidelity observations, flow parameters were determined, with multiple initializations performed to assess uncertainty and uniqueness. The workflow showed high precision in retrieving key flow characteristics and accounted for multiscale effects. Additionally, the accuracy and computational efficiency of the proposed approach were compared with an existing method, numerical simulation, demonstrating superiority by order(s) of magnitude in both aspects. To our knowledge, this is the first workflow to efficiently solve inverse modeling of multiphase flow in fractured porous media using noisy, multifidelity real-world data.