def add_phase_interconnections(net, snow_partitioning_n, voxel_size=1,
                               marching_cubes_area=False,
                               alias=None):
    r"""
    This function connects networks of two or more phases together by
    interconnecting neighbouring nodes inside different phases.

    The resulting network can be used for the study of transport and kinetics
    at interphase of two phases.

    Parameters
    ----------
    network : 2D or 3D network
        A dictoionary containing structural information of two or more
        phases networks. The dictonary format must be same as porespy
        region_to_network function.

    snow_partitioning_n : tuple
        The output generated by snow_partitioning_n function. The tuple should
        have phases_max_labels and original image of material.

    voxel_size : scalar
        The resolution of the image, expressed as the length of one side of a
        voxel, so the volume of a voxel would be **voxel_size**-cubed.  The
        default is 1, which is useful when overlaying the PNM on the original
        image since the scale of the image is alway 1 unit lenth per voxel.

    marching_cubes_area : bool
        If ``True`` then the surface area and interfacial area between regions
        will be causing the marching cube algorithm. This is a more accurate
        representation of area in extracted network, but is quite slow, so
        it is ``False`` by default.  The default method simply counts voxels
        so does not correctly account for the voxelated nature of the images.

    alias : dict (Optional)
        A dictionary that assigns unique image label to specific phase.
        For example {1: 'Solid'} will show all structural properties associated
        with label 1 as Solid phase properties.
        If ``None`` then default labelling will be used i.e {1: 'Phase1',..}.

    Returns
    -------
    A dictionary containing network information of individual and connected
    networks. The dictionary names use the OpenPNM convention so it may be
    converted directly to an OpenPNM network object using the ``update``
    command.

    """
    # -------------------------------------------------------------------------
    # Get alias if provided by user
    im = snow_partitioning_n.im
    al = _create_alias_map(im, alias=alias)
    # -------------------------------------------------------------------------
    # Find interconnection and interfacial area between ith and jth phases
    conns1 = net['throat.conns'][:, 0]
    conns2 = net['throat.conns'][:, 1]
    label = net['pore.label'] - 1

    num = snow_partitioning_n.phase_max_label
    num = [0, *num]
    phases_num = sp.unique(im * 1)
    phases_num = sp.trim_zeros(phases_num)
    for i in phases_num:
        loc1 = sp.logical_and(conns1 >= num[i - 1], conns1 < num[i])
        loc2 = sp.logical_and(conns2 >= num[i - 1], conns2 < num[i])
        loc3 = sp.logical_and(label >= num[i - 1], label < num[i])
        net['throat.{}'.format(al[i])] = loc1 * loc2
        net['pore.{}'.format(al[i])] = loc3
        if i == phases_num[-1]:
            loc4 = sp.logical_and(conns1 < num[-1], conns2 >= num[-1])
            loc5 = label >= num[-1]
            net['throat.boundary'] = loc4
            net['pore.boundary'] = loc5
        for j in phases_num:
            if j > i:
                pi_pj_sa = sp.zeros_like(label)
                loc6 = sp.logical_and(conns2 >= num[j - 1], conns2 < num[j])
                pi_pj_conns = loc1 * loc6
                net['throat.{}_{}'.format(al[i], al[j])] = pi_pj_conns
                if any(pi_pj_conns):
                    # ---------------------------------------------------------
                    # Calculates phase[i] interfacial area that connects with
                    # phase[j] and vice versa
                    p_conns = net['throat.conns'][:, 0][pi_pj_conns]
                    s_conns = net['throat.conns'][:, 1][pi_pj_conns]
                    ps = net['throat.area'][pi_pj_conns]
                    p_sa = sp.bincount(p_conns, ps)
                    # trim zeros at head/tail position to avoid extra bins
                    p_sa = sp.trim_zeros(p_sa)
                    i_index = sp.arange(min(p_conns), max(p_conns) + 1)
                    j_index = sp.arange(min(s_conns), max(s_conns) + 1)
                    s_pa = sp.bincount(s_conns, ps)
                    s_pa = sp.trim_zeros(s_pa)
                    pi_pj_sa[i_index] = p_sa
                    pi_pj_sa[j_index] = s_pa
                    # ---------------------------------------------------------
                    # Calculates interfacial area using marching cube method
                    if marching_cubes_area:
                        ps_c = net['throat.area'][pi_pj_conns]
                        p_sa_c = sp.bincount(p_conns, ps_c)
                        p_sa_c = sp.trim_zeros(p_sa_c)
                        s_pa_c = sp.bincount(s_conns, ps_c)
                        s_pa_c = sp.trim_zeros(s_pa_c)
                        pi_pj_sa[i_index] = p_sa_c
                        pi_pj_sa[j_index] = s_pa_c
                    net['pore.{}_{}_area'.format(al[i],
                                                 al[j])] = (pi_pj_sa *
                                                            voxel_size ** 2)
    return net