The permeability of fractures, including natural and hydraulic, are essential parameters for the modeling of fluid flow in conventional and unconventional fractured reservoirs. However, traditional analytical cubic law (CL-based) models used to estimate fracture permeability show unsatisfactory performance when dealing with different dynamic complexities of fractures. This work presents a data-driven, physics-included model based on machine learning as an alternative to traditional methods. The workflow for the development of the data-driven model includes four steps. Step 1: Identify uncertain parameters and perform Latin Hypercube Sampling (LHS). We first identify the uncertain parameters which affect the fracture permeability. We then generate training samples using LHS. Step 2: Perform training simulations and collect inputs and outputs. In this step, high-resolution simulations with parallel computing for the Navier-Stokes equations (NSEs) are run for each of the training samples. We then collect the inputs and outputs from the simulations. Step 3: Construct an optimized data-driven surrogate model. A data-driven model based on machine learning is then built to model the nonlinear mapping between the inputs and outputs collected from Step 2. Herein, Artificial Neural Network (ANN) coupling with Bayesian optimization algorithm is implemented to obtain the optimized surrogate model. Step 4: Validate the proposed data-driven model. In this step, we conduct blind validation on the proposed model with high-fidelity simulations. We further test the developed surrogate model with newly generated fracture cases with a broad range of roughness and tortuosity under different Reynolds numbers. We then compare its performance to the reference NSEs solutions. Results show that the developed data-driven model delivers good accuracy exceeding 90% for all training, validation, and test samples. This work introduces an integrated workflow for developing a data-driven, physics-included model using machine learning to estimate fracture permeability under complex physics (e.g., inertial effect). To our knowledge, this technique is introduced for the first time for the upscaling of rock fractures. The proposed model offers an efficient and accurate alternative to the traditional upscaling methods that can be readily implemented in reservoir characterization and modeling workflows.