## Materials and Methods
Porosity and pore volume were measured using helium and the permeability was measured using nitrogen. The confining pressure for the measurements was 500 psi. The microCT images were acquired in the SkyScan 1173 equipment (Bruker, Kontich, Belgium). Image reconstruction was done using the InstaRecon(r) software (InstaRecon, Champaign, IL, USA), and image processing was done using Avizo software (Thermo Fisher Scientific, Waltham, MA, USA). For the NMR experiments, the Oxford MQC--5 MHz (Oxford Instruments, Abingdon, UK) equipment was used to acquire the transverse relaxation time (T2) using a CPMG sequence from [PERSON] et al. [53]. For that, the sample was saturated with NaCl synthetic brine with a 50,000 ppm concentration. Equilibrium contact angle (CA) and interfacial tension (IFT) were measured in a Drop Shape Analysis (DSA) Hastelloy high-pressure system (Kruss, Hamburg, Germany and Eurotechnica, Bargteheide, Germany) equipped with high-precision pumps, pressure indicators, and a temperature-controlling device. All experiments were carried out at 60 \({}^{\circ}\)C and 6.895 MPa (1000 psi). Rock samples were cleaned using toluene and methanol, and CA was measured in both clean and aged rocks. The latter were aged for one day in formation brine and 15 days in crude oil at 60 \({}^{\circ}\)C. For the tests considering carbonated brine, both fluids (oil and brine) were saturated with CO\({}_{2}\) at test conditions prior to the experiment to reduce mass transfer effects. A complete description of CA and IFT measurement experimental procedures both in the presence and absence of CO\({}_{2}\) can be found in prior works [38; 43]. Primary drainage capillary pressure curves were obtained using a Core Lab Instruments A-200 centrifuge. The coquina plugs were saturated in brine and immersed in oil. The speed of rotation was increased in steps (from 500 to 9000 rpm), recording the volume of brine produced from the rock on each step. The speed of rotation was converted to capillary pressure using the equations presented by [PERSON] and [PERSON] [54], and the inverse method described in [PERSON] et al. [55] was applied to analyze the data. The simulations were performed using PORE software, developed by Engineering Simulation and Scientific Software (ESSS) Corporation, Florianopolis, Brazil. PORE constructs a pore network model from 3D microCT images, which are binarized based on a user-selected threshold. The binary image is then developed into a pore network model of pores and throats based on the Maximum Spheres [56; 57]. The shapes of the pores and throats are then established based on their Form Factor \"G\" which is defined as the cross-sectional area divided by the cross-sectional perimeter squared [58]. Euclidean geometrical shapes such as circles, squares, equilateral triangles, and prisms each have a unique form factor, and the pores and throats are individually transformed into the closest matching geometrical shape [59]. Multi-phase flow is then simulated based on the conservation of momentum with increasing pressure steps sequentially from pore to pore through the connecting throats where the radii of curvature in the Young-Laplace equation is controlled by the geometrical shapes of the pores and throats [60; 61]. The contact angle and the interfacial tension are user-input parameters and were selected based on lab experiments, and are presented in Section 3.2.