diff --git a/ProcessLib/HeatTransportBHE/BHE/BHETypes.h b/ProcessLib/HeatTransportBHE/BHE/BHETypes.h
index d0f9cec124afa2827fdb05733c63d1d08f5ea123..e284f64a3b24e8976b83a2f9c1aba32b3c255ef1 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHETypes.h
+++ b/ProcessLib/HeatTransportBHE/BHE/BHETypes.h
@@ -13,6 +13,8 @@
 
 #include <boost/variant.hpp>
 #include "BHE_1U.h"
+#include "BHE_CXA.h"
+#include "BHE_CXC.h"
 
 namespace ProcessLib
 {
@@ -20,7 +22,7 @@ namespace HeatTransportBHE
 {
 namespace BHE
 {
-using BHETypes = boost::variant<BHE_1U>;
+using BHETypes = boost::variant<BHE_1U, BHE_CXA, BHE_CXC>;
 }  // end of namespace BHE
 }  // end of namespace HeatTransportBHE
 }  // end of namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h
index 4f8c67d13732ec77002971f202d15cec6d7379ed..672115c9cd8fd62036ff45395bb539fd05418b7c 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h
@@ -46,6 +46,7 @@ public:
            PipeConfiguration1U const& pipes);
 
     static constexpr int number_of_unknowns = 4;
+    static constexpr int number_of_grout_zones = 2;
 
     std::array<double, number_of_unknowns> pipeHeatCapacities() const;
 
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_CXA.cpp b/ProcessLib/HeatTransportBHE/BHE/BHE_CXA.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2a305b3765ad268bc9a62ff038274a6920a015bf
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_CXA.cpp
@@ -0,0 +1,193 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "BHE_CXA.h"
+
+#include <boost/math/constants/constants.hpp>
+#include "FlowAndTemperatureControl.h"
+#include "Physics.h"
+#include "ThermoMechanicalFlowProperties.h"
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+namespace BHE
+{
+BHE_CXA::BHE_CXA(BoreholeGeometry const& borehole,
+                 RefrigerantProperties const& refrigerant,
+                 GroutParameters const& grout,
+                 FlowAndTemperatureControl const& flowAndTemperatureControl,
+                 PipeConfigurationCXA const& pipes)
+    : BHECommon{borehole, refrigerant, grout, flowAndTemperatureControl},
+      _pipes(pipes)
+{
+    // Initialize thermal resistances.
+    auto values = apply_visitor(
+        [&](auto const& control) {
+            return control(refrigerant.reference_temperature,
+                           0. /* initial time */);
+        },
+        flowAndTemperatureControl);
+    updateHeatTransferCoefficients(values.flow_rate);
+}
+
+std::array<double, BHE_CXA::number_of_unknowns> BHE_CXA::pipeHeatCapacities()
+    const
+{
+    double const& rho_r = refrigerant.density;
+    double const& specific_heat_capacity = refrigerant.specific_heat_capacity;
+    double const& rho_g = grout.rho_g;
+    double const& porosity_g = grout.porosity_g;
+    double const& heat_cap_g = grout.heat_cap_g;
+
+    return {{/*i*/ rho_r * specific_heat_capacity,
+             /*o*/ rho_r * specific_heat_capacity,
+             /*g*/ (1.0 - porosity_g) * rho_g * heat_cap_g}};
+}
+
+std::array<double, BHE_CXA::number_of_unknowns> BHE_CXA::pipeHeatConductions()
+    const
+{
+    double const& lambda_r = refrigerant.thermal_conductivity;
+    double const& rho_r = refrigerant.density;
+    double const& Cp_r = refrigerant.specific_heat_capacity;
+    double const& alpha_L = _pipes.longitudinal_dispersion_length;
+    double const& porosity_g = grout.porosity_g;
+    double const& lambda_g = grout.lambda_g;
+
+    double const velocity_norm_inner_pipe = std::abs(_flow_velocity);
+    double const velocity_norm_annulus = std::abs(_flow_velocity_annulus);
+
+    // Here we calculate the laplace coefficients in the governing
+    // equations of BHE. These governing equations can be found in
+    // 1) Diersch (2013) FEFLOW book on page 952, M.120-122, or
+    // 2) Diersch (2011) Comp & Geosci 37:1122-1135, Eq. 23-25.
+    return {{// pipe i, Eq. 23
+             (lambda_r + rho_r * Cp_r * alpha_L * velocity_norm_annulus),
+             // pipe o, Eq. 24
+             (lambda_r + rho_r * Cp_r * alpha_L * velocity_norm_inner_pipe),
+             // pipe g, Eq. 25
+             (1.0 - porosity_g) * lambda_g}};
+}
+
+std::array<Eigen::Vector3d, BHE_CXA::number_of_unknowns>
+BHE_CXA::pipeAdvectionVectors() const
+{
+    double const& rho_r = refrigerant.density;
+    double const& Cp_r = refrigerant.specific_heat_capacity;
+
+    return {{// pipe i, Eq. 23
+             {0, 0, -rho_r * Cp_r * _flow_velocity_annulus},
+             // pipe o, Eq. 24
+             {0, 0, rho_r * Cp_r * _flow_velocity},
+             // grout g, Eq. 25
+             {0, 0, 0}}};
+}
+
+constexpr std::pair<int, int> BHE_CXA::inflow_outflow_bc_component_ids[];
+
+void BHE_CXA::updateHeatTransferCoefficients(double const flow_rate)
+
+{
+    auto const tm_flow_properties_annulus =
+        calculateThermoMechanicalFlowPropertiesAnnulus(
+            _pipes.inner_outflow_pipe,
+            _pipes.outer_pipe,
+            borehole_geometry.length,
+            refrigerant,
+            flow_rate);
+
+    _flow_velocity_annulus = tm_flow_properties_annulus.velocity;
+
+    auto const tm_flow_properties =
+        calculateThermoMechanicalFlowPropertiesPipe(_pipes.inner_outflow_pipe,
+                                                    borehole_geometry.length,
+                                                    refrigerant,
+                                                    flow_rate);
+
+    _flow_velocity = tm_flow_properties.velocity;
+    _thermal_resistances =
+        calcThermalResistances(tm_flow_properties.nusselt_number,
+                               tm_flow_properties_annulus.nusselt_number);
+}
+
+/// Nu_o is the Nusselt number of inner pipe, Nu_i is the Nusselt number of
+/// annulus.
+std::array<double, BHE_CXA::number_of_unknowns> BHE_CXA::calcThermalResistances(
+    double const Nu_o, double const Nu_i)
+{
+    static constexpr double pi = boost::math::constants::pi<double>();
+
+    double const& D = borehole_geometry.diameter;
+    double const& lambda_r = refrigerant.thermal_conductivity;
+    double const& lambda_g = grout.lambda_g;
+    double const& d_o = _pipes.outer_pipe.diameter;
+    double const& d_i = _pipes.inner_outflow_pipe.diameter;
+    double const& b_o = _pipes.outer_pipe.wall_thickness;
+    double const& b_i = _pipes.inner_outflow_pipe.wall_thickness;
+    double const& lambda_p_o = _pipes.outer_pipe.wall_thermal_conductivity;
+    double const& lambda_p_i =
+        _pipes.inner_outflow_pipe.wall_thermal_conductivity;
+
+    // thermal resistances due to advective flow of refrigerant in the _pipes
+    // Eq. 58, 59, 60 in Diersch_2011_CG
+    double const d_h = d_o - d_i - 2 * b_i;
+    double const d_i_o = d_i + 2 * b_i;
+    double const d_o_o = d_o + 2 * b_o;
+    double const R_adv_o = 1.0 / (Nu_o * lambda_r * pi);
+    double const R_adv_a_i = 1.0 / (Nu_i * lambda_r * pi) * (d_h / d_i_o);
+    double const R_adv_b_i = 1.0 / (Nu_i * lambda_r * pi) * (d_h / d_o);
+
+    // thermal resistance due to thermal conductivity of the pipe wall material
+    // Eq. 66
+    double const R_con_o = std::log(d_i_o / d_i) / (2.0 * pi * lambda_p_i);
+    double const R_con_i = std::log(d_o_o / d_o) / (2.0 * pi * lambda_p_o);
+
+    // Eq. 68
+    double const chi =
+        std::log(std::sqrt(D * D + d_o_o * d_o_o) / std::sqrt(2) / d_o_o) /
+        std::log(D / d_o_o);
+    // Eq. 69
+    // thermal resistances of the grout
+    double const R_g = std::log(D / d_o_o) / 2 / pi / lambda_g;
+    // Eq. 67
+    // thermal resistance due to the grout transition.
+    double const R_con_b = chi * R_g;
+    // Eq. 56 and 57
+    double const R_ff = R_adv_o + R_adv_a_i + R_con_o;
+    double const R_fig = R_adv_b_i + R_con_i + R_con_b;
+
+    // thermal resistance due to grout-soil exchange
+    double const R_gs = (1 - chi) * R_g;
+
+    return {{R_fig, R_ff, R_gs}};
+
+    // keep the following lines------------------------------------------------
+    // when debuffing the code, printing the R and phi values are needed--------
+    // std::cout << "Rfig =" << R_fig << " Rff =" << R_ff << " Rgs =" << R_gs
+    // << "\n"; double phi_fig = 1.0 / (R_fig * S_i);
+    // double phi_ff = 1.0 / (R_ff * S_g1); double phi_gs = 1.0 / (R_gs * S_gs);
+    // std::cout << "phi_fig ="
+    // << phi_fig << " phi_ff =" << phi_ff << " phi_gs =" << phi_gs << "\n";
+    // -------------------------------------------------------------------------
+}
+
+double BHE_CXA::updateFlowRateAndTemperature(double const T_out,
+                                             double const current_time)
+{
+    auto values = apply_visitor(
+        [&](auto const& control) { return control(T_out, current_time); },
+        flowAndTemperatureControl);
+    updateHeatTransferCoefficients(values.flow_rate);
+    return values.temperature;
+}
+}  // namespace BHE
+}  // namespace HeatTransportBHE
+}  // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_CXA.h b/ProcessLib/HeatTransportBHE/BHE/BHE_CXA.h
new file mode 100644
index 0000000000000000000000000000000000000000..5cdadeb6c04bdc3f674969125b979ae7d4451c25
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_CXA.h
@@ -0,0 +1,151 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <Eigen/Eigen>
+
+#include "BaseLib/Error.h"
+
+#include "BHECommon.h"
+#include "FlowAndTemperatureControl.h"
+#include "PipeConfigurationCXA.h"
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+namespace BHE
+{
+/**
+ * The BHE_CXA class is the realization of Coaxial pipe with Annular type of the
+ * Borehole Heate Exchanger. In this class, the pipe heat capacity,
+ * pipe heat conduction, pipe advection vectors are intialized according to the
+ * geometry of CXA type of BHE. For CXA type of BHE, 3 primary unknowns are
+ * assigned on the 1D BHE elements. They are the temperature in inflow pipe
+ * T_in, temperature in outflow pipe T_out, temperature of the grout zone
+ * surrounding the inflow pipe T_g. These primary variables are solved
+ * according to heat convection and conduction equations on the pipes and also
+ * in the grout zone. The interaction of the CXA type of BHE and the
+ * surrounding soil is regulated through the thermal resistance values, which
+ * are calculated specifically during the initialization of the class.
+ */
+class BHE_CXA final : public BHECommon
+{
+public:
+    BHE_CXA(BoreholeGeometry const& borehole,
+            RefrigerantProperties const& refrigerant,
+            GroutParameters const& grout,
+            FlowAndTemperatureControl const& flowAndTemperatureControl,
+            PipeConfigurationCXA const& pipes);
+
+    static constexpr int number_of_unknowns = 3;
+    static constexpr int number_of_grout_zones = 1;
+
+    std::array<double, number_of_unknowns> pipeHeatCapacities() const;
+
+    std::array<double, number_of_unknowns> pipeHeatConductions() const;
+
+    std::array<Eigen::Vector3d, number_of_unknowns> pipeAdvectionVectors()
+        const;
+
+    template <int NPoints, typename SingleUnknownMatrixType,
+              typename RMatrixType, typename RPiSMatrixType,
+              typename RSMatrixType>
+    void assembleRMatrices(
+        int const idx_bhe_unknowns,
+        Eigen::MatrixBase<SingleUnknownMatrixType> const& matBHE_loc_R,
+        Eigen::MatrixBase<RMatrixType>& R_matrix,
+        Eigen::MatrixBase<RPiSMatrixType>& R_pi_s_matrix,
+        Eigen::MatrixBase<RSMatrixType>& R_s_matrix) const
+    {
+        switch (idx_bhe_unknowns)
+        {
+            case 0:  // PHI_fig
+                R_matrix.block(0, 2 * NPoints, NPoints, NPoints) +=
+                    -1.0 * matBHE_loc_R;
+                R_matrix.block(2 * NPoints, 0, NPoints, NPoints) +=
+                    -1.0 * matBHE_loc_R;
+
+                R_matrix.block(0, 0, NPoints, NPoints) +=
+                    1.0 * matBHE_loc_R;  // K_i
+                R_matrix.block(2 * NPoints,
+                               2 * NPoints,
+                               NPoints,
+                               NPoints) += 1.0 * matBHE_loc_R;  // K_ig
+                return;
+            case 1:  // PHI_ff
+                R_matrix.block(0, NPoints, NPoints, NPoints) +=
+                    -1.0 * matBHE_loc_R;
+                R_matrix.block(NPoints, 0, NPoints, NPoints) +=
+                    -1.0 * matBHE_loc_R;
+
+                R_matrix.block(0, 0, NPoints,
+                               NPoints) += 1.0 * matBHE_loc_R;  // K_i1
+                R_matrix.block(NPoints, NPoints, NPoints, NPoints) +=
+                    1.0 * matBHE_loc_R;  // K_o
+                return;
+            case 2:  // PHI_gs
+                R_s_matrix += matBHE_loc_R;
+
+                R_pi_s_matrix.block(2 * NPoints, 0, NPoints, NPoints) +=
+                    -1.0 * matBHE_loc_R;
+
+                R_matrix.block(2 * NPoints, 2 * NPoints, NPoints,
+                               NPoints) += matBHE_loc_R;  // K_ig
+                return;
+            default:
+                OGS_FATAL(
+                    "Error!!! In the function BHE_CXA::assembleRMatrices, "
+                    "the index of bhe unknowns is out of range! ");
+        }
+    }
+
+    /// Return the inflow temperature for the boundary condition.
+    double updateFlowRateAndTemperature(double T_out, double current_time);
+
+    double thermalResistance(int const unknown_index) const
+    {
+        return _thermal_resistances[unknown_index];
+    }
+
+    static constexpr std::pair<int, int> inflow_outflow_bc_component_ids[] = {
+        {0, 1}};
+
+private:
+    // Placing it here before using it in the cross_section_areas.
+    PipeConfigurationCXA const _pipes;
+
+public:
+    std::array<double, number_of_unknowns> const cross_section_areas = {
+        {_pipes.outer_pipe.area() - _pipes.inner_outflow_pipe.outsideArea(),
+         _pipes.inner_outflow_pipe.area(),
+         borehole_geometry.area() - _pipes.outer_pipe.outsideArea()}};
+
+private:
+    void updateHeatTransferCoefficients(double const flow_rate);
+
+    std::array<double, number_of_unknowns> calcThermalResistances(
+        double const Nu_o, double const Nu_i);
+
+private:
+    /// PHI_fig, PHI_ff, PHI_gs;
+    /// Here we store the thermal resistances needed for computation of the heat
+    /// exchange coefficients in the governing equations of BHE.
+    /// These governing equations can be found in
+    /// 1) Diersch (2013) FEFLOW book on page 958, M.3, or
+    /// 2) Diersch (2011) Comp & Geosci 37:1122-1135, Eq. 90-97.
+    std::array<double, number_of_unknowns> _thermal_resistances;
+
+    /// Flow velocity inside the pipes and annulus. Depends on the flow_rate.
+    double _flow_velocity, _flow_velocity_annulus;
+};
+}  // namespace BHE
+}  // namespace HeatTransportBHE
+}  // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_CXC.cpp b/ProcessLib/HeatTransportBHE/BHE/BHE_CXC.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b13df15c940b722fb796190e635103ce82b21cc2
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_CXC.cpp
@@ -0,0 +1,192 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "BHE_CXC.h"
+
+#include <boost/math/constants/constants.hpp>
+#include "FlowAndTemperatureControl.h"
+#include "Physics.h"
+#include "ThermoMechanicalFlowProperties.h"
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+namespace BHE
+{
+BHE_CXC::BHE_CXC(BoreholeGeometry const& borehole,
+                 RefrigerantProperties const& refrigerant,
+                 GroutParameters const& grout,
+                 FlowAndTemperatureControl const& flowAndTemperatureControl,
+                 PipeConfigurationCXC const& pipes)
+    : BHECommon{borehole, refrigerant, grout, flowAndTemperatureControl},
+      _pipes(pipes)
+{
+    // Initialize thermal resistances.
+    auto values = apply_visitor(
+        [&](auto const& control) {
+            return control(refrigerant.reference_temperature,
+                           0. /* initial time */);
+        },
+        flowAndTemperatureControl);
+    updateHeatTransferCoefficients(values.flow_rate);
+}
+
+std::array<double, BHE_CXC::number_of_unknowns> BHE_CXC::pipeHeatCapacities()
+    const
+{
+    double const& rho_r = refrigerant.density;
+    double const& specific_heat_capacity = refrigerant.specific_heat_capacity;
+    double const& rho_g = grout.rho_g;
+    double const& porosity_g = grout.porosity_g;
+    double const& heat_cap_g = grout.heat_cap_g;
+
+    return {{/*i*/ rho_r * specific_heat_capacity,
+             /*o*/ rho_r * specific_heat_capacity,
+             /*g*/ (1.0 - porosity_g) * rho_g * heat_cap_g}};
+}
+
+std::array<double, BHE_CXC::number_of_unknowns> BHE_CXC::pipeHeatConductions()
+    const
+{
+    double const& lambda_r = refrigerant.thermal_conductivity;
+    double const& rho_r = refrigerant.density;
+    double const& Cp_r = refrigerant.specific_heat_capacity;
+    double const& alpha_L = _pipes.longitudinal_dispersion_length;
+    double const& porosity_g = grout.porosity_g;
+    double const& lambda_g = grout.lambda_g;
+
+    double const velocity_norm_inner_pipe = std::abs(_flow_velocity);
+    double const velocity_norm_annulus = std::abs(_flow_velocity_annulus);
+
+    // Here we calculate the laplace coefficients in the governing
+    // equations of BHE. These governing equations can be found in
+    // 1) Diersch (2013) FEFLOW book on page 952, M.120-122, or
+    // 2) Diersch (2011) Comp & Geosci 37:1122-1135, Eq. 26-28.
+    return {{// pipe i, Eq. 26
+             (lambda_r + rho_r * Cp_r * alpha_L * velocity_norm_inner_pipe),
+             // pipe o, Eq. 27
+             (lambda_r + rho_r * Cp_r * alpha_L * velocity_norm_annulus),
+             // pipe g, Eq. 28
+             (1.0 - porosity_g) * lambda_g}};
+}
+
+std::array<Eigen::Vector3d, BHE_CXC::number_of_unknowns>
+BHE_CXC::pipeAdvectionVectors() const
+{
+    double const& rho_r = refrigerant.density;
+    double const& Cp_r = refrigerant.specific_heat_capacity;
+
+    return {{// pipe i, Eq. 26
+             {0, 0, -rho_r * Cp_r * _flow_velocity},
+             // pipe o, Eq. 27
+             {0, 0, rho_r * Cp_r * _flow_velocity_annulus},
+             // grout g, Eq. 28
+             {0, 0, 0}}};
+}
+
+constexpr std::pair<int, int> BHE_CXC::inflow_outflow_bc_component_ids[];
+
+void BHE_CXC::updateHeatTransferCoefficients(double const flow_rate)
+
+{
+    auto const tm_flow_properties_annulus =
+        calculateThermoMechanicalFlowPropertiesAnnulus(_pipes.inner_inflow_pipe,
+                                                       _pipes.outer_pipe,
+                                                       borehole_geometry.length,
+                                                       refrigerant,
+                                                       flow_rate);
+
+    _flow_velocity_annulus = tm_flow_properties_annulus.velocity;
+
+    auto const tm_flow_properties =
+        calculateThermoMechanicalFlowPropertiesPipe(_pipes.inner_inflow_pipe,
+                                                    borehole_geometry.length,
+                                                    refrigerant,
+                                                    flow_rate);
+
+    _flow_velocity = tm_flow_properties.velocity;
+    _thermal_resistances =
+        calcThermalResistances(tm_flow_properties_annulus.nusselt_number,
+                               tm_flow_properties.nusselt_number);
+}
+
+/// Nu_o is the Nusselt number of inner pipe, Nu_i is the Nusselt number of
+/// annulus.
+std::array<double, BHE_CXC::number_of_unknowns> BHE_CXC::calcThermalResistances(
+    double const Nu_annulus, double const Nu_inner_inflow_pipe)
+{
+    static constexpr double pi = boost::math::constants::pi<double>();
+
+    double const& D = borehole_geometry.diameter;
+    double const& lambda_r = refrigerant.thermal_conductivity;
+    double const& lambda_g = grout.lambda_g;
+    double const& d_o = _pipes.outer_pipe.diameter;
+    double const& d_i = _pipes.inner_inflow_pipe.diameter;
+    double const& b_o = _pipes.outer_pipe.wall_thickness;
+    double const& b_i = _pipes.inner_inflow_pipe.wall_thickness;
+    double const& lambda_p_o = _pipes.outer_pipe.wall_thermal_conductivity;
+    double const& lambda_p_i =
+        _pipes.inner_inflow_pipe.wall_thermal_conductivity;
+
+    // thermal resistances due to advective flow of refrigerant in the _pipes
+    // Eq. 73, 74, 75 in Diersch_2011_CG
+    double const d_h = d_o - d_i - 2 * b_i;
+    double const d_i_o = d_i + 2 * b_i;
+    double const d_o_o = d_o + 2 * b_o;
+    double const R_adv_i = 1.0 / (Nu_inner_inflow_pipe * lambda_r * pi);
+    double const R_adv_a_o = 1.0 / (Nu_annulus * lambda_r * pi) * (d_h / d_i_o);
+    double const R_adv_b_o = 1.0 / (Nu_annulus * lambda_r * pi) * (d_h / d_o);
+
+    // thermal resistance due to thermal conductivity of the pipe wall material
+    // Eq. 81
+    double const R_con_o = std::log(d_o_o / d_o) / (2.0 * pi * lambda_p_o);
+    double const R_con_i = std::log(d_i_o / d_i) / (2.0 * pi * lambda_p_i);
+
+    // Eq. 83
+    double const chi =
+        std::log(std::sqrt(D * D + d_o_o * d_o_o) / std::sqrt(2) / d_o_o) /
+        std::log(D / d_o_o);
+    // Eq. 84
+    // thermal resistances of the grout
+    double const R_g = std::log(D / d_o_o) / 2 / pi / lambda_g;
+    // Eq. 82
+    // thermal resistance due to the grout transition.
+    double const R_con_b = chi * R_g;
+    // Eq. 71 and 72
+    double const R_ff = R_adv_i + R_adv_a_o + R_con_i;
+    double const R_fog = R_adv_b_o + R_con_o + R_con_b;
+
+    // thermal resistance due to grout-soil exchange
+    double const R_gs = (1 - chi) * R_g;
+
+    return {{R_ff, R_fog, R_gs}};
+
+    // keep the following lines------------------------------------------------
+    // when debuffing the code, printing the R and phi values are needed--------
+    // std::cout << "Rfog =" << R_fog << " Rff =" << R_ff << " Rgs =" << R_gs
+    // << "\n"; double phi_fog = 1.0 / (R_fog * S_i);
+    // double phi_ff = 1.0 / (R_ff * S_g1); double phi_gs = 1.0 / (R_gs * S_gs);
+    // std::cout << "phi_fog ="
+    // << phi_fog << " phi_ff =" << phi_ff << " phi_gs =" << phi_gs << "\n";
+    // -------------------------------------------------------------------------
+}
+
+double BHE_CXC::updateFlowRateAndTemperature(double const T_out,
+                                             double const current_time)
+{
+    auto values = apply_visitor(
+        [&](auto const& control) { return control(T_out, current_time); },
+        flowAndTemperatureControl);
+    updateHeatTransferCoefficients(values.flow_rate);
+    return values.temperature;
+}
+}  // namespace BHE
+}  // namespace HeatTransportBHE
+}  // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_CXC.h b/ProcessLib/HeatTransportBHE/BHE/BHE_CXC.h
new file mode 100644
index 0000000000000000000000000000000000000000..9c40e9b71917519de70545a845d6bf4e7341787e
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_CXC.h
@@ -0,0 +1,151 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <Eigen/Eigen>
+
+#include "BaseLib/Error.h"
+
+#include "BHECommon.h"
+#include "FlowAndTemperatureControl.h"
+#include "PipeConfigurationCXC.h"
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+namespace BHE
+{
+/**
+ * The BHE_CXC class is the realization of Coaxial pipe with Centred type of the
+ * Borehole Heate Exchanger. In this class, the pipe heat capacity,
+ * pipe heat conduction, pipe advection vectors are intialized according to the
+ * geometry of CXC type of BHE. For CXC type of BHE, 3 primary unknowns are
+ * assigned on the 1D BHE elements. They are the temperature in inflow pipe
+ * T_in, temperature in outflow pipe T_out, temperature of the grout zone
+ * surrounding the outflow pipe T_g. These primary variables are solved
+ * according to heat convection and conduction equations on the pipes and also
+ * in the grout zone. The interaction of the CXC type of BHE and the
+ * surrounding soil is regulated through the thermal resistance values, which
+ * are calculated specifically during the initialization of the class.
+ */
+class BHE_CXC final : public BHECommon
+{
+public:
+    BHE_CXC(BoreholeGeometry const& borehole,
+            RefrigerantProperties const& refrigerant,
+            GroutParameters const& grout,
+            FlowAndTemperatureControl const& flowAndTemperatureControl,
+            PipeConfigurationCXC const& pipes);
+
+    static constexpr int number_of_unknowns = 3;
+    static constexpr int number_of_grout_zones = 1;
+
+    std::array<double, number_of_unknowns> pipeHeatCapacities() const;
+
+    std::array<double, number_of_unknowns> pipeHeatConductions() const;
+
+    std::array<Eigen::Vector3d, number_of_unknowns> pipeAdvectionVectors()
+        const;
+
+    template <int NPoints, typename SingleUnknownMatrixType,
+              typename RMatrixType, typename RPiSMatrixType,
+              typename RSMatrixType>
+    void assembleRMatrices(
+        int const idx_bhe_unknowns,
+        Eigen::MatrixBase<SingleUnknownMatrixType> const& matBHE_loc_R,
+        Eigen::MatrixBase<RMatrixType>& R_matrix,
+        Eigen::MatrixBase<RPiSMatrixType>& R_pi_s_matrix,
+        Eigen::MatrixBase<RSMatrixType>& R_s_matrix) const
+    {
+        switch (idx_bhe_unknowns)
+        {
+            case 0:  // PHI_ff
+                R_matrix.block(0, NPoints, NPoints, NPoints) +=
+                    -1.0 * matBHE_loc_R;
+                R_matrix.block(NPoints, 0, NPoints, NPoints) +=
+                    -1.0 * matBHE_loc_R;
+
+                R_matrix.block(0, 0, NPoints,
+                               NPoints) += 1.0 * matBHE_loc_R;  // K_i
+                R_matrix.block(NPoints, NPoints, NPoints, NPoints) +=
+                    1.0 * matBHE_loc_R;  // K_o
+                return;
+            case 1:  // PHI_fog
+                R_matrix.block(NPoints, 2 * NPoints, NPoints, NPoints) +=
+                    -1.0 * matBHE_loc_R;
+                R_matrix.block(2 * NPoints, NPoints, NPoints, NPoints) +=
+                    -1.0 * matBHE_loc_R;
+
+                R_matrix.block(NPoints, NPoints, NPoints, NPoints) +=
+                    1.0 * matBHE_loc_R;  // K_o
+                R_matrix.block(2 * NPoints,
+                               2 * NPoints,
+                               NPoints,
+                               NPoints) += 1.0 * matBHE_loc_R;  // K_og
+                return;
+            case 2:  // PHI_gs
+                R_s_matrix += matBHE_loc_R;
+
+                R_pi_s_matrix.block(2 * NPoints, 0, NPoints, NPoints) +=
+                    -1.0 * matBHE_loc_R;
+
+                R_matrix.block(2 * NPoints, 2 * NPoints, NPoints,
+                               NPoints) += matBHE_loc_R;  // K_og
+                return;
+            default:
+                OGS_FATAL(
+                    "Error!!! In the function BHE_CXC::assembleRMatrices, "
+                    "the index of bhe unknowns is out of range! ");
+        }
+    }
+
+    /// Return the inflow temperature for the boundary condition.
+    double updateFlowRateAndTemperature(double T_out, double current_time);
+
+    double thermalResistance(int const unknown_index) const
+    {
+        return _thermal_resistances[unknown_index];
+    }
+
+    static constexpr std::pair<int, int> inflow_outflow_bc_component_ids[] = {
+        {0, 1}};
+
+private:
+    // Placing it here before using it in the cross_section_areas.
+    PipeConfigurationCXC const _pipes;
+
+public:
+    std::array<double, number_of_unknowns> const cross_section_areas = {
+        {_pipes.inner_inflow_pipe.area(),
+         _pipes.outer_pipe.area() - _pipes.inner_inflow_pipe.outsideArea(),
+         borehole_geometry.area() - _pipes.outer_pipe.outsideArea()}};
+
+private:
+    void updateHeatTransferCoefficients(double const flow_rate);
+
+    std::array<double, number_of_unknowns> calcThermalResistances(
+        double const Nu_o, double const Nu_i);
+
+private:
+    /// PHI_ff, PHI_fog, PHI_gs;
+    /// Here we store the thermal resistances needed for computation of the heat
+    /// exchange coefficients in the governing equations of BHE.
+    /// These governing equations can be found in
+    /// 1) Diersch (2013) FEFLOW book on page 958, M.3, or
+    /// 2) Diersch (2011) Comp & Geosci 37:1122-1135, Eq. 90-97.
+    std::array<double, number_of_unknowns> _thermal_resistances;
+
+    /// Flow velocity inside the pipes and annulus. Depends on the flow_rate.
+    double _flow_velocity, _flow_velocity_annulus;
+};
+}  // namespace BHE
+}  // namespace HeatTransportBHE
+}  // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/CreateBHECXA.cpp b/ProcessLib/HeatTransportBHE/BHE/CreateBHECXA.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..571042f2210d3badded73931dd7d6b8a3a9334f1
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/BHE/CreateBHECXA.cpp
@@ -0,0 +1,101 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "CreateBHECXA.h"
+#include "BaseLib/ConfigTree.h"
+#include "CreateFlowAndTemperatureControl.h"
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+namespace BHE
+{
+BHE::BHE_CXA createBHECXA(
+    BaseLib::ConfigTree const& config,
+    std::map<std::string,
+             std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
+        curves)
+{
+    //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__borehole}
+    auto const& borehole_config = config.getConfigSubtree("borehole");
+    const double borehole_length =
+        //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__borehole__length}
+        borehole_config.getConfigParameter<double>("length");
+    const double borehole_diameter =
+        //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__borehole__diameter}
+        borehole_config.getConfigParameter<double>("diameter");
+    BoreholeGeometry const borehole{borehole_length, borehole_diameter};
+
+    //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__pipes}
+    auto const& pipes_config = config.getConfigSubtree("pipes");
+    //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__pipes__outer}
+    Pipe const outer_pipe = createPipe(pipes_config.getConfigSubtree("outer"));
+    Pipe const inner_pipe =
+        //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__pipes__inner}
+        createPipe(pipes_config.getConfigSubtree("inner"));
+    const double pipe_longitudinal_dispersion_length =
+        //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__pipes__longitudinal_dispersion_length}
+        pipes_config.getConfigParameter<double>(
+            "longitudinal_dispersion_length");
+    PipeConfigurationCXA const pipes{inner_pipe, outer_pipe,
+                                     pipe_longitudinal_dispersion_length};
+
+    //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__grout}
+    auto const& grout_config = config.getConfigSubtree("grout");
+    const double grout_density =
+        //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__grout__density}
+        grout_config.getConfigParameter<double>("density");
+    const double grout_porosity =
+        //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__grout__porosity}
+        grout_config.getConfigParameter<double>("porosity");
+    const double grout_heat_capacity =
+        //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__grout__heat_capacity}
+        grout_config.getConfigParameter<double>("heat_capacity");
+    const double grout_thermal_conductivity =
+        //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__grout__thermal_conductivity}
+        grout_config.getConfigParameter<double>("thermal_conductivity");
+    GroutParameters const grout{grout_density, grout_porosity,
+                                grout_heat_capacity,
+                                grout_thermal_conductivity};
+
+    //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__refrigerant}
+    auto const& refrigerant_config = config.getConfigSubtree("refrigerant");
+    double const refrigerant_density =
+        //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__refrigerant__density}
+        refrigerant_config.getConfigParameter<double>("density");
+    double const refrigerant_viscosity =
+        //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__refrigerant__viscosity}
+        refrigerant_config.getConfigParameter<double>("viscosity");
+    double const refrigerant_heat_capacity =
+        //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__refrigerant__specific_heat_capacity}
+        refrigerant_config.getConfigParameter<double>("specific_heat_capacity");
+    double const refrigerant_thermal_conductivity =
+        //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__refrigerant__thermal_conductivity}
+        refrigerant_config.getConfigParameter<double>("thermal_conductivity");
+    double const refrigerant_reference_temperature =
+        //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__refrigerant__reference_temperature}
+        refrigerant_config.getConfigParameter<double>("reference_temperature");
+    RefrigerantProperties const refrigerant{
+        refrigerant_viscosity, refrigerant_density,
+        refrigerant_thermal_conductivity, refrigerant_heat_capacity,
+        refrigerant_reference_temperature};
+
+    auto const flowAndTemperatureControl = createFlowAndTemperatureControl(
+        //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__flow_and_temperature_control}
+        config.getConfigSubtree("flow_and_temperature_control"),
+        curves,
+        refrigerant);
+
+    return {borehole, refrigerant, grout, flowAndTemperatureControl, pipes};
+}
+}  // namespace BHE
+}  // namespace HeatTransportBHE
+}  // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/CreateBHECXA.h b/ProcessLib/HeatTransportBHE/BHE/CreateBHECXA.h
new file mode 100644
index 0000000000000000000000000000000000000000..3534a8eb8bb2c683afe771f919eeaff655660803
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/BHE/CreateBHECXA.h
@@ -0,0 +1,31 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "BHE_CXA.h"
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+namespace BHE
+{
+BHE::BHE_CXA createBHECXA(
+    BaseLib::ConfigTree const& bhe_conf,
+    std::map<std::string,
+             std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
+        curves);
+}  // namespace BHE
+}  // namespace HeatTransportBHE
+}  // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/CreateBHECXC.cpp b/ProcessLib/HeatTransportBHE/BHE/CreateBHECXC.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5852eea80981e3e4b3b437c9de14b372a9b8559b
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/BHE/CreateBHECXC.cpp
@@ -0,0 +1,101 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "CreateBHECXC.h"
+#include "BaseLib/ConfigTree.h"
+#include "CreateFlowAndTemperatureControl.h"
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+namespace BHE
+{
+BHE::BHE_CXC createBHECXC(
+    BaseLib::ConfigTree const& config,
+    std::map<std::string,
+             std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
+        curves)
+{
+    //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__borehole}
+    auto const& borehole_config = config.getConfigSubtree("borehole");
+    const double borehole_length =
+        //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__borehole__length}
+        borehole_config.getConfigParameter<double>("length");
+    const double borehole_diameter =
+        //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__borehole__diameter}
+        borehole_config.getConfigParameter<double>("diameter");
+    BoreholeGeometry const borehole{borehole_length, borehole_diameter};
+
+    //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__pipes}
+    auto const& pipes_config = config.getConfigSubtree("pipes");
+    //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__pipes__inner}
+    Pipe const inner_pipe = createPipe(pipes_config.getConfigSubtree("inner"));
+    Pipe const outer_pipe =
+        //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__pipes__outer}
+        createPipe(pipes_config.getConfigSubtree("outer"));
+    const double pipe_longitudinal_dispersion_length =
+        //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__pipes__longitudinal_dispersion_length}
+        pipes_config.getConfigParameter<double>(
+            "longitudinal_dispersion_length");
+    PipeConfigurationCXC const pipes{inner_pipe, outer_pipe,
+                                     pipe_longitudinal_dispersion_length};
+
+    //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__grout}
+    auto const& grout_config = config.getConfigSubtree("grout");
+    const double grout_density =
+        //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__grout__density}
+        grout_config.getConfigParameter<double>("density");
+    const double grout_porosity =
+        //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__grout__porosity}
+        grout_config.getConfigParameter<double>("porosity");
+    const double grout_heat_capacity =
+        //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__grout__heat_capacity}
+        grout_config.getConfigParameter<double>("heat_capacity");
+    const double grout_thermal_conductivity =
+        //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__grout__thermal_conductivity}
+        grout_config.getConfigParameter<double>("thermal_conductivity");
+    GroutParameters const grout{grout_density, grout_porosity,
+                                grout_heat_capacity,
+                                grout_thermal_conductivity};
+
+    //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__refrigerant}
+    auto const& refrigerant_config = config.getConfigSubtree("refrigerant");
+    double const refrigerant_density =
+        //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__refrigerant__density}
+        refrigerant_config.getConfigParameter<double>("density");
+    double const refrigerant_viscosity =
+        //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__refrigerant__viscosity}
+        refrigerant_config.getConfigParameter<double>("viscosity");
+    double const refrigerant_heat_capacity =
+        //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__refrigerant__specific_heat_capacity}
+        refrigerant_config.getConfigParameter<double>("specific_heat_capacity");
+    double const refrigerant_thermal_conductivity =
+        //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__refrigerant__thermal_conductivity}
+        refrigerant_config.getConfigParameter<double>("thermal_conductivity");
+    double const refrigerant_reference_temperature =
+        //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__refrigerant__reference_temperature}
+        refrigerant_config.getConfigParameter<double>("reference_temperature");
+    RefrigerantProperties const refrigerant{
+        refrigerant_viscosity, refrigerant_density,
+        refrigerant_thermal_conductivity, refrigerant_heat_capacity,
+        refrigerant_reference_temperature};
+
+    auto const flowAndTemperatureControl = createFlowAndTemperatureControl(
+        //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__flow_and_temperature_control}
+        config.getConfigSubtree("flow_and_temperature_control"),
+        curves,
+        refrigerant);
+
+    return {borehole, refrigerant, grout, flowAndTemperatureControl, pipes};
+}
+}  // namespace BHE
+}  // namespace HeatTransportBHE
+}  // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/CreateBHECXC.h b/ProcessLib/HeatTransportBHE/BHE/CreateBHECXC.h
new file mode 100644
index 0000000000000000000000000000000000000000..0aabfc4194420b31925dcfb50e8c8339acd61964
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/BHE/CreateBHECXC.h
@@ -0,0 +1,31 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include "BHE_CXC.h"
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+namespace BHE
+{
+BHE::BHE_CXC createBHECXC(
+    BaseLib::ConfigTree const& bhe_conf,
+    std::map<std::string,
+             std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
+        curves);
+}  // namespace BHE
+}  // namespace HeatTransportBHE
+}  // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/CreateFlowAndTemperatureControl.cpp b/ProcessLib/HeatTransportBHE/BHE/CreateFlowAndTemperatureControl.cpp
index 2015ea78dd431b0c76353b7cbdd6033423d63cc6..80d1795e44ab59122b5bb435c3277e980aea15c2 100644
--- a/ProcessLib/HeatTransportBHE/BHE/CreateFlowAndTemperatureControl.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/CreateFlowAndTemperatureControl.cpp
@@ -9,9 +9,9 @@
  *
  */
 
+#include "BaseLib/Algorithm.h"
 #include "BaseLib/ConfigTree.h"
 #include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"
-#include "BaseLib/Algorithm.h"
 
 #include "CreateFlowAndTemperatureControl.h"
 #include "RefrigerantProperties.h"
@@ -37,10 +37,10 @@ FlowAndTemperatureControl createFlowAndTemperatureControl(
         double const flow_rate = config.getConfigParameter<double>("flow_rate");
 
         auto const& temperature_curve = *BaseLib::getOrError(
-                    curves,
-                    //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__flow_and_temperature_control__TemperatureCurveConstantFlow__temperature_curve}
-                    config.getConfigParameter<std::string>("temperature_curve"),
-                    "Required temperature curve not found.");
+            curves,
+            //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__flow_and_temperature_control__TemperatureCurveConstantFlow__temperature_curve}
+            config.getConfigParameter<std::string>("temperature_curve"),
+            "Required temperature curve not found.");
 
         return TemperatureCurveConstantFlow{flow_rate, temperature_curve};
     }
@@ -59,10 +59,10 @@ FlowAndTemperatureControl createFlowAndTemperatureControl(
     if (type == "FixedPowerFlowCurve")
     {
         auto const& flow_rate_curve = *BaseLib::getOrError(
-                    curves,
-                    //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__flow_and_temperature_control__FixedPowerFlowCurve__flow_rate_curve}
-                    config.getConfigParameter<std::string>("flow_rate_curve"),
-                    "Required flow rate curve not found.");
+            curves,
+            //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__flow_and_temperature_control__FixedPowerFlowCurve__flow_rate_curve}
+            config.getConfigParameter<std::string>("flow_rate_curve"),
+            "Required flow rate curve not found.");
 
         //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__flow_and_temperature_control__FixedPowerFlowCurve__power}
         double const power = config.getConfigParameter<double>("power");
@@ -71,6 +71,22 @@ FlowAndTemperatureControl createFlowAndTemperatureControl(
                                    refrigerant.specific_heat_capacity,
                                    refrigerant.density};
     }
+
+    if (type == "PowerCurveConstantFlow")
+    {
+        auto const& power_curve = *BaseLib::getOrError(
+            curves,
+            //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__flow_and_temperature_control__PowerCurveConstantFlow__power_curve}
+            config.getConfigParameter<std::string>("power_curve"),
+            "Required power curve not found.");
+
+        //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__flow_and_temperature_control__PowerCurveConstantFlow__flow_rate}
+        double const flow_rate = config.getConfigParameter<double>("flow_rate");
+
+        return PowerCurveConstantFlow{power_curve, flow_rate,
+                                      refrigerant.specific_heat_capacity,
+                                      refrigerant.density};
+    }
     OGS_FATAL("FlowAndTemperatureControl type '%s' is not implemented.",
               type.c_str());
 }
diff --git a/ProcessLib/HeatTransportBHE/BHE/FlowAndTemperatureControl.h b/ProcessLib/HeatTransportBHE/BHE/FlowAndTemperatureControl.h
index f4de023261b9a608fff38c8663120f7c5f316845..8be20bde1b8fa55a5023955800c143dc4266e4b9 100644
--- a/ProcessLib/HeatTransportBHE/BHE/FlowAndTemperatureControl.h
+++ b/ProcessLib/HeatTransportBHE/BHE/FlowAndTemperatureControl.h
@@ -64,9 +64,24 @@ struct FixedPowerFlowCurve
     double density;
 };
 
+struct PowerCurveConstantFlow
+{
+    FlowAndTemperature operator()(double const T_out, double const time) const
+    {
+        double const power = power_curve.getValue(time);
+        return {flow_rate, power / flow_rate / heat_capacity / density + T_out};
+    }
+    MathLib::PiecewiseLinearInterpolation const& power_curve;
+
+    double flow_rate;
+    double heat_capacity;
+    double density;
+};
+
 using FlowAndTemperatureControl = boost::variant<TemperatureCurveConstantFlow,
                                                  FixedPowerConstantFlow,
-                                                 FixedPowerFlowCurve>;
+                                                 FixedPowerFlowCurve,
+                                                 PowerCurveConstantFlow>;
 
 }  // namespace BHE
 }  // namespace HeatTransportBHE
diff --git a/ProcessLib/HeatTransportBHE/BHE/PipeConfigurationCXA.h b/ProcessLib/HeatTransportBHE/BHE/PipeConfigurationCXA.h
new file mode 100644
index 0000000000000000000000000000000000000000..bf94568419a103d6659108184f903e91c10df387
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/BHE/PipeConfigurationCXA.h
@@ -0,0 +1,34 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <boost/math/constants/constants.hpp>
+
+#include "BaseLib/ConfigTree.h"
+#include "Pipe.h"
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+namespace BHE
+{
+struct PipeConfigurationCXA
+{
+    Pipe const inner_outflow_pipe;
+    Pipe const outer_pipe;
+
+    double const longitudinal_dispersion_length;
+};
+}  // namespace BHE
+}  // namespace HeatTransportBHE
+}  // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/PipeConfigurationCXC.h b/ProcessLib/HeatTransportBHE/BHE/PipeConfigurationCXC.h
new file mode 100644
index 0000000000000000000000000000000000000000..3344d9e4c690453c7fd757699fa0d21501837d82
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/BHE/PipeConfigurationCXC.h
@@ -0,0 +1,34 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <boost/math/constants/constants.hpp>
+
+#include "BaseLib/ConfigTree.h"
+#include "Pipe.h"
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+namespace BHE
+{
+struct PipeConfigurationCXC
+{
+    Pipe const inner_inflow_pipe;
+    Pipe const outer_pipe;
+
+    double const longitudinal_dispersion_length;
+};
+}  // namespace BHE
+}  // namespace HeatTransportBHE
+}  // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/CreateHeatTransportBHEProcess.cpp b/ProcessLib/HeatTransportBHE/CreateHeatTransportBHEProcess.cpp
index ce0730be389245c894f9cf9763201c304b5f45cd..e26fff25b166b284f7463e27c5b112b23beff633 100644
--- a/ProcessLib/HeatTransportBHE/CreateHeatTransportBHEProcess.cpp
+++ b/ProcessLib/HeatTransportBHE/CreateHeatTransportBHEProcess.cpp
@@ -16,6 +16,8 @@
 
 #include "BHE/BHETypes.h"
 #include "BHE/CreateBHE1U.h"
+#include "BHE/CreateBHECXA.h"
+#include "BHE/CreateBHECXC.h"
 #include "HeatTransportBHEProcess.h"
 #include "HeatTransportBHEProcessData.h"
 
@@ -186,6 +188,18 @@ std::unique_ptr<Process> createHeatTransportBHEProcess(
             bhes.push_back(BHE::createBHE1U(bhe_config, curves));
             continue;
         }
+
+        if (bhe_type == "CXA")
+        {
+            bhes.push_back(BHE::createBHECXA(bhe_config, curves));
+            continue;
+        }
+
+        if (bhe_type == "CXC")
+        {
+            bhes.push_back(BHE::createBHECXC(bhe_config, curves));
+            continue;
+        }
         OGS_FATAL("Unknown BHE type '%s'.", bhe_type.c_str());
     }
     // end of reading BHE parameters -------------------------------------------
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE-impl.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE-impl.h
index 9d4a201947ab69f2b848accd923dd84ef6372e66..e30497703e01df48a8226bf7f09494338aa76e25 100644
--- a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE-impl.h
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE-impl.h
@@ -28,7 +28,8 @@ HeatTransportBHELocalAssemblerBHE<ShapeFunction, IntegrationMethod, BHEType>::
     : _process_data(process_data),
       _integration_method(integration_order),
       _bhe(bhe),
-      _element_id(e.getID())
+      _element_id(e.getID()),
+      _number_of_grout_zones(bhe.number_of_grout_zones)
 {
     // need to make sure that the BHE elements are one-dimensional
     assert(e.getDimension() == 1);
@@ -185,7 +186,7 @@ void HeatTransportBHELocalAssemblerBHE<ShapeFunction, IntegrationMethod,
     local_K
         .template block<temperature_size, temperature_size>(temperature_index,
                                                             temperature_index)
-        .noalias() += 2.0 * _R_s_matrix;
+        .noalias() += _number_of_grout_zones * _R_s_matrix;
 
     // debugging
     // std::string sep =
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE.h
index 55741410a1e9b5e4dfa9dd30243ed244c77d940c..0336f5fc542b018f00024919c5bbf66aff85c35b 100644
--- a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE.h
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE.h
@@ -86,6 +86,8 @@ private:
 
     std::size_t const _element_id;
 
+    const int _number_of_grout_zones;
+
     SecondaryData<typename ShapeMatrices::ShapeType> _secondary_data;
 
     typename ShapeMatricesType::template MatrixType<bhe_unknowns_size,
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/LocalDataInitializer.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/LocalDataInitializer.h
index f1ece8767289f2564232e89d091afeed09bed818..a7f19a74b031a1dd9b26e5622b2940f29618478a 100644
--- a/ProcessLib/HeatTransportBHE/LocalAssemblers/LocalDataInitializer.h
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/LocalDataInitializer.h
@@ -261,6 +261,20 @@ private:
                     e, boost::get<BHE::BHE_1U>(bhe),
                     std::forward<ConstructorArgs>(args)...}};
             }
+
+            if (bhe.type() == typeid(BHE::BHE_CXA))
+            {
+                return LADataIntfPtr{new LADataBHE<ShapeFunction, BHE::BHE_CXA>{
+                    e, boost::get<BHE::BHE_CXA>(bhe),
+                    std::forward<ConstructorArgs>(args)...}};
+            }
+
+            if (bhe.type() == typeid(BHE::BHE_CXC))
+            {
+                return LADataIntfPtr{new LADataBHE<ShapeFunction, BHE::BHE_CXC>{
+                    e, boost::get<BHE::BHE_CXC>(bhe),
+                    std::forward<ConstructorArgs>(args)...}};
+            }
             OGS_FATAL(
                 "Trying to create local assembler for an unknown BHE type.");
         };